Late last week I was chasing down some numerical precision questions I had about certain color-space computations. I found myself trying to determine, once and for all, what tristimulus values to use for the CIE standard illuminant called D65.

CIE defines D65 using a table of spectral values, 1 nm spacing, over the range 300 nm to 830 nm.

So how do you get the corresponding tristimulus values for the CIE 1931 standard color observer? The CIE technical report on colorimetry, CIE 15:2004, suggests that the “rigorous” way to do it is using a dot product of the 1 nm spectral values with the 1 nm xbar, ybar, and zbar tables for the standard observer. But the technical report contains written tristimulus values that do not quite agree with that computation. And there are other values in use out there as well.

This table shows everything I was able to track down:

If you round all these values to five fractional digits, there are 10 unique triples among the 13 rows shown above. Even if you round to just four fractional digits, there are 6 unique triples.

Sigh.

]]>I got the idea for these test matrices from Gene Golub years ago, but the mathematical foundation comes from a paper by Divakar Viswanath and Nick Trefethen.

This post is about two MATLAB programs included in the collection from *Numerical Computing with MATLAB*. The program `golub` generates badly conditioned integer test matrices. The program `lugui` is an interactive graphical interface that allows you to experiment with pivot strategies in Gaussian elimination.

Set the random number generator to the state it has when MATLAB starts and generate a 6-by-6 example.

```
rng('default')
A = golub(6)
```

A = 1 3 11 0 -11 -15 18 55 209 15 -198 -277 -23 -33 144 532 259 82 9 55 405 437 -100 -285 3 -4 -111 -180 39 219 -13 -9 202 346 401 253

It's not obvious that this matrix is badly conditioned, but it is.

cond(A)

ans = 3.0115e+12

In case you needed another demonstration that the determinant does not indicate near singularity.

```
format short
det(A)
```

ans = 1.0000

The elements of the inverse vary over nine orders of magnitude and show perverse scaling that is a direct consequence of the way this matrix was created.

format short e X = inv(A)

X = 2.8508e+09 -1.5513e+08 3.1110e+06 1.4854e+06 -1.6559e+05 -2.0380e+04 -1.3362e+09 7.2711e+07 -1.4582e+06 -6.9624e+05 7.7615e+04 9.5520e+03 1.0422e+08 -5.6713e+06 1.1373e+05 5.4306e+04 -6.0540e+03 -7.4500e+02 1.2586e+07 -6.8489e+05 1.3735e+04 6.5580e+03 -7.3100e+02 -9.0000e+01 -8.4225e+05 4.5833e+04 -9.1900e+02 -4.3900e+02 4.9000e+01 6.0000e+00 -1.3830e+05 7.5260e+03 -1.5100e+02 -7.2000e+01 8.0000e+00 1.0000e+00

In a 1998 paper Divakar Viswanath and Nick Trefethen demonstrate that the 2-norm condition number of an *n* -by- *n* triangular matrix with normally distributed entries below the diagonal and ones on the diagonal grows exponentially almost surely at a rate *c^n* with *c* = 1.30568...

I want my matrices to have integer entries so I use randomly distributed integers, but the exponential behavior must be nearly the same.

Here is the code for `golub`.

```
type golub
```

function A = golub(n) %GOLUB Badly conditioned integer test matrices. % GOLUB(n) is the product of two random integer n-by-n matrices, % one of them unit lower triangular and one unit upper triangular. % LU factorization without pivoting fails to reveal that such % matrices are badly conditioned. % See also LUGUI. % Copyright 2014 Cleve Moler % Copyright 2014 The MathWorks, Inc. s = 10; L = tril(round(s*randn(n)),-1)+eye(n); U = triu(round(s*randn(n)),1)+eye(n); A = L*U;

`Lugui` lets you use the mouse to select the pivot elements for Gaussian elimination. Or, you can choose diagonal pivoting, partial pivoting or complete pivoting. Here are animated gifs of these three choices.

Diagonal pivoting actually defies the conventional wisdom and picks the *smallest* possible element at each stage as the pivot, namely 1.0. But this choice means that division by the pivot involves no roundoff. Moreover the choice leads to the reconstruction of the original triangular factors.

```
[L,U] = lugui_noshow(A,'diagonal')
cond_L = cond(L)
cond_U = cond(U)
```

L = 1 0 0 0 0 0 18 1 0 0 0 0 -23 36 1 0 0 0 9 28 -2 1 0 0 3 -13 -1 7 1 0 -13 30 15 16 -8 1 U = 1 3 11 0 -11 -15 0 1 11 15 0 -7 0 0 1 -8 6 -11 0 0 0 1 11 24 0 0 0 0 1 -6 0 0 0 0 0 1 cond_L = 8.4806e+06 cond_U = 7.9889e+05

Partial pivoting gives only a weak indication that this matrix is badly conditioned.

```
[L,U,p] = lugui_noshow(A,'partial');
condition_estimate_partial = max(abs(diag(U)))/min(abs(diag(U)))
```

condition_estimate_partial = 5.8208e+06

Complete pivoting does a much better job.

```
[L,U,p,q] = lugui_noshow(A,'complete');
condition_estimate_complete = max(abs(diag(U)))/min(abs(diag(U)))
```

condition_estimate_complete = 1.5166e+12

The QR decomposition needs column pivoting to reveal the ill conditioning.

[Q,R] = qr(A); condition_estimate_without = max(abs(diag(R)))/min(abs(diag(R))) [Q,R,E] = qr(A); condition_estimate_with = max(abs(diag(R)))/min(abs(diag(R)))

condition_estimate_without = 6.6062e+06 condition_estimate_with = 2.2595e+12

Divakar Viswanath and Nick Trefethen, "Condition Numbers of Random Triangular Matrices," *SIAM J. Matrix Anal. Appl.* 19, 564-581, 1998, <http://epubs.siam.org/doi/pdf/10.1137/S0895479896312869>. Also available at https://people.maths.ox.ac.uk/trefethen/publication/PDF/1998_77.pdf

Cleve Moler, *Numerical Computing with MATLAB*, 2004. <http://www.mathworks.com/moler/ncmfilelist.html>.

Get
the MATLAB code

Published with MATLAB® R2015a

Brett‘s Pick this week is struct2str, by Marco Cococcioni.

Back in May, Sean wrote about different File Exchange submissions that facilitate writing cell arrays to text. That was a useful post–and one that triggered an interesting follow-up discussion between Sean and me: namely, how would one write a struct to string? I write a lot of apps, and frequently want to put the contents of a struct into a listbox, for example. That can be challenging and frustrating.

When I suggested that Sean follow up with a post on writing structs to text, he quickly agreed that that would be interesting. But I think he initially missed my point: he suggested that I could easily write the fieldnames of a structure to a listbox using, for instance,

lb = uicontrol('units','normalized',...'position',[0.1 0.1 0.8 0.8],... 'Style','listbox'); d = dir('.\*'); lb.String = fieldnames(d);

“True,” I replied. “But I want the contents of the struct, not just the fieldnames.” His reply: “How about”

names = {d(:).name}; lb.String = names;

Okay, that gets me the values of a single fieldname (in this example, ‘name’). However, assume that I get the file information for an image, and want to drop it all in a textbox or listbox:

```
structInfo = imfinfo('street2.jpg')
```

structInfo = Filename: 'C:\Program Files\MATLAB\R2015a\toolbox\matlab\demos\...' FileModDate: '22-Mar-2004 19:54:40' FileSize: 39920 Format: 'jpg' FormatVersion: '' Width: 640 Height: 480 BitDepth: 24 ColorType: 'truecolor' FormatSignature: '' NumberOfSamples: 3 CodingMethod: 'Huffman' CodingProcess: 'Sequential' Comment: {}

How could I easily put all of that information into a listbox?

That’s where Marco’s file comes in:

lb.String = struct2str(structInfo);

I had initially done this with a particularly ugly `for` loop, but `struct2str` makes it very easy. (Marco used some loops, too, but his code alleviates for me the pain of writing that loop. Thanks, Marco!)

Note that there are a few empty fields in the `imfinfo` output struct in the usage above. Kudos and swag to the first person who replies with an elegant approach to automatically removing from the struct the fields with empty values. Use any approach you’d like, including other File Exchange submissions. What I’d like to see is this:

structInfo = rmfield(structInfo,{'FormatVersion','FormatSignature','Comment'})

structInfo = Filename: 'C:\Program Files\MATLAB\R2015a\toolbox\matlab\demos\...' FileModDate: '22-Mar-2004 19:54:40' FileSize: 39920 Format: 'jpg' Width: 640 Height: 480 BitDepth: 24 ColorType: 'truecolor' NumberOfSamples: 3 CodingMethod: 'Huffman' CodingProcess: 'Sequential'

…but without the manual specification of the empty fields. Let me know how you would approach this, and I will select (and reward) the one I like best. It’s good to be the blogger. ;)

As always, comments are welcome! Let me know what you think (or respond to my challenge) here, or leave feedback for Marco here.

Get
the MATLAB code

Published with MATLAB® R2015a

I just got asked a question about a good way to find the closest value in a vector that was less than a threshold. My solution is fairly short, and demonstrates some of my favorite MATLAB techniques. I will compare also show you an "obvious" solution.

First let's set up the data for our problem.

thresh = 75; nvals = 10^6; data = 100*rand(1,nvals);

We could solve this by brute force, just looping over the values. Let's try that. I'm going to set the index value to empty (`[]`). That way, if we end up with an array that doesn't meet the criterion, we can tell. Also, I am setting the current minimum value to `-Inf` so any finite value that we find as a valid candidate will be closer to `thresh`, assuming we can find one.

loopindex = []; candidate = -inf; for ind = 1:numel(data) dval = data(ind); if dval < thresh && dval > candidate candidate = dval; loopindex = ind; end end

Next I show you my non-loop solution.

First collect the list of possible data entries - the ones that are less than the threshold value `thresh`. This list is a logical variable, essentially true and false values for each entry in `data`, selecting the candidate values that are less than the `thresh` value.

possibles = data < thresh;

Let's find the actual best value, plus its index into the reduced set from `possibles`. The index we find will not be the index into `data` but rather into a smaller array which is the subset meeting the threshold criteria.

[posmax, posind] = max(data(possibles));

Convert the answer into the correct index in the *original* array, `data`.

inddatapos = find(possibles); % possible indices inddata = inddatapos(posind); % find the index we care about

If `inddata` is empty, then there were no possible values meeting the criterion. So we set the final values accordingly.

if isempty(inddata) posmax = -Inf; end

sameSolution = isequal([inddata posmax],[loopindex candidate])

sameSolution = 1

You might guess correctly which solution is more natural to me at this point :-). I am wondering which solution you prefer, and why? Let me know here.

Get
the MATLAB code

Published with MATLAB® R2015a

A customer contacted Tech Support recently to ask the following question about structuring elements and dilation:

*If I created a structuring element from the matrix [1 1; 1 1], then the origin of the structuring element is the (1,1) or upper-left element. How do I create a structuring element so that the origin is at (1,2) or (2,1)?*

There are two equivalent ways to accomplish this.

In the first way, you add a row or column (or both) to the matrix that you pass to `strel` or `imdilate` so that the center of the new matrix is the desired structuring element origin. For example:

se = strel([0 1 1; 0 1 1; 0 0 0])

se = Flat STREL object containing 4 neighbors. Neighborhood: 0 1 1 0 1 1 0 0 0

In the second way, you pass the original 2-by-2 matrix to `strel` and then call the function `translate` to "move" the origin to the desired location.

se = strel([1 1; 1 1])

se = Flat STREL object containing 4 neighbors. Neighborhood: 1 1 1 1

se2 = translate(se,[-1 0])

se2 = Flat STREL object containing 4 neighbors. Neighborhood: 0 1 1 0 1 1 0 0 0

Get
the MATLAB code

Published with MATLAB® R2015a

**The Problem: Getting the Value of Nested Mask Parameters**

I want to write a script that will look in a model for blocks of a particular type, a subsystem I made, and analyze the values passed to it as dialog parameters. I do not want this script to be an initialization callback, I just want the users to be able to run it whenever they want, like a Model Advisor check.

When my block is at the top of a model, it is very easy. Let's say I have variables "a" and "b" defined in the base workspace with respective values of 2 and 3. If the user specifies "a+b" as value of parameter "c" in my block, I can get the final value of "c" using the code below:

But what if the user puts my block inside another subsystem:

In this case, "a" and "b" are not in the base workspace anymore, they are in the mask workspace of subsystem SS_Level1.

**The Solution**

If you look carefully at the bottom of the reference page for Mask Parameters, you will find one named `maskWSVariables`

**Now it's your turn**

Let us know what you think of this technique, or if you have preferred tricks to access values and properties in a model by leaving a comment here.

]]>A blog reader asked me recently how to visualize the ellipse-based measurements computed by `regionprops`. The reader wanted to superimpose the estimated ellipses on the image.

To refresh your memory: The function `regionprops`, which computes geometrical measurements of image regions, offers several measurements that are based on fitting an ellipse to the region. (If you must know, the ellipse fit is computing by matching second-order moments.) These measurements are:

- MajorAxisLength
- MinorAxisLength
- Orientation

I would plot an ellipse using a parametric equation. Conveniently, the Wikipedia article on ellipses has a parametric equation in a form that's just right:

$$X(t) = X_c + a \cos t \cos \phi - b \sin t \sin \phi$$

$$Y(t) = Y_c + a \cos t \sin \phi + b \sin t \cos \phi$$

where $(X_c,Y_c)$ is the ellipse center, $a$ and $b$ are the major and minor axis lengths, and $\phi$ is the angle between the x-axis and the major axis.

Let's try it with this image, which contains a bunch of ellipse-like objects:

```
url = 'http://blogs.mathworks.com/steve/files/rice-bw.png';
bw = imread(url);
imshow(bw)
```

Ask `regionprops` to compute all the ellipse-related measurements:

s = regionprops(bw,{... 'Centroid',... 'MajorAxisLength',... 'MinorAxisLength',... 'Orientation'})

s = 69x1 struct array with fields: Centroid MajorAxisLength MinorAxisLength Orientation

Here's what the measurements look like for the seventh object:

s(7)

ans = Centroid: [37.9286 226.9732] MajorAxisLength: 30.9655 MinorAxisLength: 9.7466 Orientation: -64.8058

Now write a loop that computes the ellipse curves one at a time, for each object, and superimposes the curve on the image.

figure imshow(bw,'InitialMagnification','fit') t = linspace(0,2*pi,50); hold on for k = 1:length(s) a = s(k).MajorAxisLength/2; b = s(k).MinorAxisLength/2; Xc = s(k).Centroid(1); Yc = s(k).Centroid(2); phi = deg2rad(-s(k).Orientation); x = Xc + a*cos(t)*cos(phi) - b*sin(t)*sin(phi); y = Yc + a*cos(t)*sin(phi) + b*sin(t)*cos(phi); plot(x,y,'r','Linewidth',5) end hold off

That's it!

Enjoy.

Get
the MATLAB code

Published with MATLAB® R2015a

Sean's pick this week is geom2d by David Legland.

Working at MathWorks, I get to use MATLAB every day, learn a lot about it, and get to try everything. But one of the things I forget is the experience of someone who doesn't use it every day but occasionally needs it for a specific task. Thus, when I get the opportunity to solve a problem like this, it reminds me of my early days and reminds me what a privilege it is to have access to professional quality software.

Yesterday, one of our program managers came to my office and wanted me to create a map of the United States for her with different regions color coded. Seems straight forward right? Here's the story of the steps I took to create the map for her.

I knew off the top of my head, that `usamap` in the Mapping Toolbox was the function that I needed to build the map. So I went to the doc page linked above by running:

`doc usamap`

The third example at the bottom of the doc page was exactly what I was looking for. I almost always start with examples because it's easier to modify or delete than to add. As an added bonus, it already even excluded Alaska and Hawaii which she didn't want on her map!

I copied this example in, and ran it. I quickly realized that the line I really needed to modify was the following one; it even had a nice comment noting the current randomness.

faceColors = makesymbolspec('Polygon',... {'INDEX', [1 numel(states)], 'FaceColor', ... polcmap(numel(states))}); %NOTE - colors are random

I figured `polcmap` probably worked like any of the other colormap functions and takes the number of unique colors to make. My task involved
six regions and so I needed six colors.

figure polcmap6 = polcmap(6); colormap(polcmap6) colorbar

Minor victory!

Now I needed to index into this colormap each state's region. I asked for an Excel sheet containing the states and the region indices.

This can be easily read in with `readtable`.

`Tstateregions = readtable('StateRegionsLookup.xlsx');`

Warning: Variable names were modified to make them valid MATLAB identifiers.

With the region indices, I should be able to plot the map color coded by indexing into the colormap:

figure ax = usamap('conus'); states = shaperead('usastatelo', 'UseGeoCoords', true,... 'Selector',... {@(name) ~any(strcmp(name,{'Alaska','Hawaii'})), 'Name'}); faceColors = makesymbolspec('Polygon',... {'INDEX', [1 numel(states)], 'FaceColor', ... polcmap6(Tstateregions.Region,:)}); geoshow(ax, states, 'DisplayType', 'polygon', ... 'SymbolSpec', faceColors) framem off gridm off mlabel off plabel off

But that doesn't look right! So what went wrong, I took a look at the states and compared them:

table({states(:).Name}.', Tstateregions.State,'VariableNames',{'MATLAB','Excel'})

ans = MATLAB Excel ______________________ ________________________ 'Alabama' ''Alabama'' 'Arizona' ''Arizona'' 'Arkansas' ''Arkansas'' 'California' ''California'' 'Colorado' ''Colorado'' 'Connecticut' ''Connecticut'' 'Delaware' ''Delaware'' 'Florida' ''District of Columbia'' 'Georgia' ''Florida'' 'Idaho' ''Georgia'' 'Illinois' ''Idaho'' 'Indiana' ''Illinois'' 'Iowa' ''Indiana'' 'Kansas' ''Iowa'' 'Kentucky' ''Kansas'' 'Louisiana' ''Kentucky'' 'Maine' ''Louisiana'' 'Maryland' ''Maine'' 'Massachusetts' ''Maryland'' 'Michigan' ''Massachusetts'' 'Minnesota' ''Michigan'' 'Mississippi' ''Minnesota'' 'Missouri' ''Mississippi'' 'Montana' ''Missouri'' 'Nebraska' ''Montana'' 'Nevada' ''Nebraska'' 'New Hampshire' ''Nevada'' 'New Jersey' ''New Hampshire'' 'New Mexico' ''New Jersey'' 'New York' ''New Mexico'' 'North Carolina' ''New York'' 'North Dakota' ''North Carolina'' 'Ohio' ''North Dakota'' 'Oklahoma' ''Ohio'' 'Oregon' ''Oklahoma'' 'Pennsylvania' ''Oregon'' 'Rhode Island' ''Pennsylvania'' 'South Carolina' ''Rhode Island'' 'South Dakota' ''South Carolina'' 'Tennessee' ''South Dakota'' 'Texas' ''Tennessee'' 'Utah' ''Texas'' 'Vermont' ''Utah'' 'Virginia' ''Vermont'' 'Washington' ''Virginia'' 'West Virginia' ''Washington'' 'Wisconsin' ''West Virginia'' 'Wyoming' ''Wisconsin'' 'District of Columbia' ''Wyoming''

Ahh ha! Washington DC is appended to the end of MATLAB's shipping state border file where as it's alphabetical in the Excel sheet causing an offset. Simple enough manual fix:

Tstateregions = [Tstateregions([1:7 9:end],:); Tstateregions(8,:)];

figure ax = usamap('conus'); faceColors = makesymbolspec('Polygon',... {'INDEX', [1 numel(states)], 'FaceColor', ... polcmap6(Tstateregions.Region,:)}); geoshow(ax, states, 'DisplayType', 'polygon', ... 'SymbolSpec', faceColors) framem off gridm off mlabel off plabel off

By now you're probably wondering what this has to do with `geom2d`. Well I still had to put the two letter code on each state. To do this I needed a way to identify the centroid of the polygon.
I went to the File Exchange and searched for "polygon centroid", and sure enough `geom2d` has a function named exactly that!

for ii = 1:numel(states) latii = states(ii).Lat; % Remove nans indicating end of border lonii = states(ii).Lon; cent = polygonCentroid(latii(~isnan(latii)),lonii(~isnan(lonii))); textm(cent(1),cent(2),Tstateregions.x2letter(ii)) end

And there you have it.

The polygon centroid function is just one of dozens in this well documented and organized toolbox. There are tools for all
sorts of math, aggregation and plotting of geometric primitives, e.g. polygons, lines, rays, and polynomial curves. Each
group of functions has its own *Contents.m* file so one can quickly see all of the functions and a quick description of what they do. The functions have a full help
description, many with an example. There are also a bunch of example files that have been published so you can see the example
without having to run it. Overall this is a great package to have around and now it can join its older brother, `geom3d` with a Pick of the Week banner!

Do you ever get a random one-off task or activity that is quickly solvable in MATLAB and that makes you feel good?

Give it a try and let us know what you think here or leave a comment for David.

Get
the MATLAB code

Published with MATLAB® R2015b

- What are the units of this signal? Radians or Degrees?
- What is the meaning of that 0.15 value I put there
- What is that signal? A velocity or a position?

The answer to that is a good naming convention.

To describe the important things to consider when defining a naming convention, I am happy to welcome guest blogger Michael Burke. Being part of the MathWorks Consulting Services team, Michael has worked with multiple customers on multiple projects and has seen tons of different naming conventions.

**What’s in a Naming Convention?**

Naming conventions are a method for providing information about variables in your model and generated code. Done well they are easy to understand and work with. However we have all run into projects using Elder Gods naming conventions — unpronounceable, indecipherable and impossible to remember.

So what makes for a good naming convention? Like flavors of ice cream this is somewhat to taste, but like good ice cream there are known good ingredients. In the case of naming conventions, information is the good ingredient.

**Ingredients**

**Who:**

*Explanation:*With the “who” you are telling what the variable is used for; it should quickly let a user/reader know if this is something he needs to review.*Examples:*HorseSpeed, HorseSp (using abbreviations)*Ice cream equivalent:*This is the most important part, this is the high fat cream.

**What (Units):**

*Explanation:*Working with physical units (kg, m, watts) or software constructs (flags, counters, states), units explain the role your data plays.*Examples:*HorseSpeedFurlongsPerFortnight, HorseSpFpf*Ice cream equivalent:*Salt, you forget that it’s in there but it makes everything so much better.

**Where (do I come from):**

*Explanation:*The “where” part of the data tells where the data is written to or originally defined. This can aid in the partitioning of the data.*Examples:*HorseSpeedRaceTrackFurlongsPerFortnight, HorseSpFpfRacet*Ice cream equivalent:*Makes debugging easier so it is like sugar, sweet.

**Role:**

*Explanation:*The role defines the data as a signal, constant or tunable parameter.*Examples:*ConstHorseSpeedRaceTrackFurlongsPerFortnight, conHorseSpFpfRacet*Ice cream equivalent:*Your base flavor such as chocolate or vanilla

**Data type:**

*Explanation:*It is not uncommon for people to include the data type (e.g. double, uint8) in the name of a variable. With the exception of logical variables, this is generally a poor design decision and should be avoided as it limits the portability of the code.*Examples:*DoubleConstHorseSpeedRaceTrackFurlongsPerFortnight, DbConHorseSpFpfRacet*Ice cream equivalent:*Jimmies, sure some 6 year old likes them, but not for long.

**Storage class:**

*Explanation:*Like data type above, including information about storage class limits the long-time portability of the data.*Examples:*LocalDoubleConstHorseSpeedRaceTrackFurlongsPerFortnight, LocDbConHorseSpFpfRacet*Ice cream equivalent:*Month-old diced peanuts sprinkled on top.

At this stage you are no doubt looking at these names with recognition of the Elder God problem mentioned above. An engineer who had been with the company for 15 years would know all the abbreviations and could decipher the meaning but for the rest of us it is a gibbering chant. The good news is that there are guidelines on how to avoid this problem.

**Mixing it up**

When creating a naming convention there are 2 primary guidelines.

**Make it human readable:** Human readability is achieved through 3 things:

- Delimitation - either using CamelCase or under_scores.
- Consistency - always have the same type of information in the same portion of the variable. (e.g. don’t have HorseSpeedMPH and MPH_DonkeySpeed.)
- Limit information - as tempting as it is to place as much information as you can into a variable, generally 2 or 3 bits of information are all it takes before it becomes overwhelming.

**Brevity is the soul of naming conventions:** Abbreviations for naming conventions, especially for units, is key to usability:

- Use standard abbreviations where possible (e.g. MPH is known to be miles per hour)
- Try to use between 2 and 4 characters for the abbreviation
- If you can’t, consider using underscores to delineate

Consider these examples of a battery current variable where “I” is the common abbreviation for current:

`BatteryI`: It is easy to miss the “I” character at the end due to the short length`Battery_I`: I now sticks out due to the “_” characters`BatterCur`: You could also use the abbreviation “Cur” for current.

**Now it's your turn**

I hope the stars are in alignment for you and that you enjoy your naming conventions!

Let us know what you think of those guidelines and if you have other guidelines not covered in this post by leaving a comment here.

]]>The set of colors that can be represented using a particular color space is called the *gamut*. Some L*a*b* color values may be *out of gamut* when converted to RGB. You know that a converted RGB color is out of gamut when any of its component values is less than 0 or greater than 1.

lab = [80 -130 85]; lab2rgb(lab)

ans = -0.6210 0.9537 -0.1926

The negative values demonstrate that the L*a*b* color [80 -130 85] is not in the gamut of the sRGB color space, which is the default RGB color space used by `lab2rgb`. A different RGB color space, Adobe RGB (1998), has a larger gamut than sRGB. Use the 'ColorSpace' parameter to convert a L*a*b* color to the Adobe color space.

lab2rgb(lab,'ColorSpace','adobe-rgb-1998')

ans = 0.1234 0.9522 0.1073

Because all the output values are between 0.0 and 1.0 (inclusive), you can conclude that the L*a*b* color [80 -130 85] is inside the Adobe RGB (1998) gamut.

Get
the MATLAB code

Published with MATLAB® R2015a

A rectangular box, such as a book or a cell phone, thrown in the air can tumble stably about its longest axis, or about its shortest axis, but not about its middle axis.

In his latest book, *Differential Equations and Linear Algebra*, my colleague Gilbert Strang describes an interesting model problem made famous among MIT students by another MIT professor, Alar Toomre. If you throw a rectangular box in the air with a twist, you can make it tumble stably about its longest axis, or about its shortest axis. But, if the lengths of the three sides of the box are different, you cannot make it tumble about its middle-sized axis.

To describe the dynamics of the box, discount the effects of gravity by using a coordinate system centered in the box and that moves with it. Let $x(t)$, $y(t)$, and $z(t)$ be the angular momenta about the three principal axes, and let $I_1$, $I_2$, and $I_3$ be the moments of inertia about these axes. Then Euler's equations are

$$ { d \over dt} \left( \begin{array}{c} x \\ y \\ z \end{array} \right) = \left( \begin{array}{r} ayz \\ bxz\\ cxy \end{array} \right) $$

where

$$ a = 1/I_3 - 1/I_2, \ \ b = 1/I_1 - 1/I_3, \ \ c = 1/I_2 - 1/I_1 $$

If we have

$$ I_1 = 1, \ \ I_2 = 1/2, \ \ I_3 = 1/3 $$

the equations become simply

$$ { d \over dt} \left( \begin{array}{c} x \\ y \\ z \end{array} \right) = \left( \begin{array}{r} yz \\ -2xz\\ xy \end{array} \right) $$

There are six critical points. Any solution to the Euler equations that starts exactly at one of these points remains there.

$$ X = \pm \left( \begin{array}{c} 1 \\ 0 \\ 0 \end{array} \right), \ \ Y = \pm \left( \begin{array}{c} 0 \\ 1 \\ 0 \end{array} \right), \ \ Z = \pm \left( \begin{array}{c} 0 \\ 0 \\ 1 \end{array} \right) $$.

It turns out that any solution that starts out with $x^2+y^2+z^2=1$ retains that property, so solutions travel on the surface of the unit sphere. If you think of the sphere as the earth with $+Z$ at the North Pole, then $+X$ is the point where the Greenwich meridian crosses the equator. This is in the eastern Atlantic, off the coast of West Africa. $+Y$ is the point where the $90^\circ$ meridian east crosses the equator. This is in Indian Ocean, west of Sumatra.

To see what happens near a critical point, we need the Jacobian.

$$ J = \left[ \begin{array}{rrr} 0 & z & y \\ -2z & 0 & -2x \\ y & x & 0 \end{array} \right] $$

At $+X$

$$ J = \left[ \begin{array}{rrr} 0 & 0 & 0 \\ 0 & 0 & -2 \\ 0 & 1 & 0 \end{array} \right] $$

The eigenvalues of $J$ are pure imaginary, namely $\lambda = \pm \sqrt{2}i$, along with $\lambda = 0$. So near $+X$ the solutions behave locally like $\cos{(\sqrt{2}t)}$ and $\sin{(\sqrt{2}t)}$. This critical point, which corresponds to rotation about the long axis, is a stable center. The critical point $-X$, which corresponds to rotation about the long axis in the opposite direction, is also a stable center.

At $+Z$

$$ J = \left[ \begin{array}{rrr} 0 & 1 & 0 \\ -2 & 0 & 0 \\ 0 & 0 & 0 \end{array} \right] $$

We have the same eigenvalues and so rotation about the short axis is also a stable center.

The middle axis, $+Y$, is a different story. Here

$$ J = \left[ \begin{array}{rrr} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \end{array} \right] $$

The eigenvalues of $J$ are real, namely $\lambda = 1, 0,$ and $-1$. Near $+Y$ and $-Y$ the solutions have components that behave locally like $e^t$. These are unstable saddle points.

The surface of the sphere can be divided into four quadrants, with one of the stable critical points, $\pm X$ or $\pm Z$, at the center of each quadrant. Any solution to the differential equation is periodic and describes an orbit that remains within one of the quadrants, circulating around its center.

If the initial point is near the center of a quadrant, the orbit is nearly elliptical and the orbit has period close to $\sqrt{2}\pi$. If the initial point is near the boundary of a quadrant, the orbit sticks near that boundary, making sharp turns at the two saddle points $\pm Y$. The period increases as the initial point gets further from the center.

I think of these as giant atmospheric currents that circulate around one-fourth of a planet.

Here is a link to a program, `tumbling_box.m`, that uses `ode45` to solve these equations and produce the following graphic. I invite you to download the program and experiment yourself. You are presented with a sphere that is empty except for dots at three critical points. You can click repeatedly to provide the initial conditions for the orbits.

You can also click on icons in the figure window toolbar to zoom in and out, and rotate and pan over the sphere. If you zoom in to the stable centers points, you can see the elliptical orbits associated with stability. If you zoom in to the saddle points, you can see the hyperbolic orbits associated with instability. If you rotate the sphere, you can follow one of the unstable orbits around the sphere, back to its starting point.

`Tumbling_box` includes `ode45` event handling code to determine the length of a period of the solution. It also has mouse button down code to interpret the coordinates of a mouse click as a point on the surface of the unit sphere.

Gilbert Strang, *Differential Equations and Linear Algebra*, Wellesley-Cambridge Press, 2014, 502pp. <http://www.wellesleycambridge.com>

Get
the MATLAB code

Published with MATLAB® R2015a

I just downloaded the MATLAB Support Package for Apple iOS Sensors and the latest MATLAB Mobile app for my iPhone.

So the first thing I did was throw and catch my phone while logging the phone's acceleration sensors.

Here's the result:

You can see the toss at 1 second, the catch at 2 seconds, and the free fall in between.

If you have an Android phone, try the MATLAB Support Package for Android Sensors.

Be sure to look at the system requirements for these support packages. You'll need the latest release of MATLAB, R2015a, to use the iOS package. You need R2014a or later to use the Android package, or you can try this earlier package to use with releases R2013a and R2013b.

Just don't blame me if you drop your phone!

]]>Some years ago we added the function `typecast` to MATLAB. I've been wanting to write about this useful little function ever since we significantly speeded it up in the R2014b release.

The function `typecast` converts numeric values from one type to another based on the underlying bytes in the numeric representation.

Clear as mud, right? Let me try a simple example to illustrate. Here is a vector of two numbers, each of which is represented using one byte:

x = uint8([128 6])

x = 128 6

The call `typecast(x,'uint16')` takes those two one-byte numbers and converts them to one 16-bit number by jamming together the original bytes.

```
y = typecast(x,'uint16')
```

y = 1664

Here's how the number 1664 is related to the original pair of numbers:

128 + 6*2^8

ans = 1664

If you give `typecast` four one-byte numbers, it will happily convert them into a pair of two-byte numbers.

p = uint8([128 6 200 3])

p = 128 6 200 3

```
typecast(p,'uint16')
```

ans = 1664 968

Or, if you prefer, it will convert them into one four-byte number:

```
typecast(p,'uint32')
```

ans = 63440512

The function `typecast` performs these conversions between all of the integer and floating-point numeric types. These days I tend to use `typecast` (together with its cousin `swapbytes`) when writing low-level binary file reading code. For files that fit easily into memory, I might just read all the file bytes at once using `fileread`, and then I'll use indexing, `swapbytes`, and `typecast` into meaningful MATLAB data.

Let's examine how long it takes to perform this conversion. First, make a vector containing 256 one-byte values.

bytes = uint8(0:255);

Now use `timeit` to estimate how long it takes to convert those values to 128 two-byte values.

```
f = @() typecast(bytes,'uint16');
t_microseconds = timeit(f) * 1e6
```

t_microseconds = 19.8762

How would you expect the computation time to change as the input vector gets longer?

Here's a plot of computation times that I measured using MATLAB R2014a.

r2014a = load('typecast speedup R2014a ah-eddins-maci'); plot(r2014a.n*256,r2014a.times*1e6) xlabel('Number of input bytes') ylabel('Conversion time (microseconds)') ax = gca; ax.YLim(1) = 0; title('typecast performance (R2014a)')

The computation time appears to be proportional to the input vector length, which is not too surprising.

And here are the results measured using the latest MATLAB (R2015a).

r2015a = load('typecast speedup R2015a ah-eddins-maci'); plot(r2015a.n*256,r2015a.times*1e6) xlabel('Number of input bytes') ylabel('Conversion time (microseconds)') ax = gca; ax.YLim(1) = 0; title('typecast performance (R2015a)')

To make the comparison easier, show the results on the same plot.

plot(r2014a.n*256,r2014a.times*1e6) hold on plot(r2015a.n*256,r2015a.times*1e6) xlabel('Number of input bytes') ylabel('Conversion time (microseconds)') ax = gca; ax.YLim(1) = 0; title('typecast performance') legend({'R2014a','R2015a'})

Now that is a little surprising! Not only is the R2015a computation time lower, it appears to be independent of the input vector length. That's because the implementation has been optimized to avoid making a copy in memory of the input vector's data. (For more about this kind of optimization, see Loren's 2006 blog post about MATLAB memory use for functions and variables.)

Have you had reason to use `typecast` in your own work? If so, please leave a comment. I'd be very interested to hear about your application.

Get
the MATLAB code

Published with MATLAB® R2015a

Jiro's pick this week is `bubbleplot` by Ameya Deoras.

Ameya was a fellow MathWorker and also helped us with this Pick of the Week blog. He's now working at a financial firm.

I'm always impressed with his MATLAB programs. They are extremely useful and very well written. One of my all-time favorites is Intelligent Dynamic Date Ticks, which I also highlighted as a Pick. This was a must-have tool, especially before R2014b, when working with data represented as dates. In R2014b, we introduced a new data type `datetime`, which makes date handling (including plotting) much easier.

Back to this post. `bubbleplot` is a visualization function for high-dimensional data. If you have ever used `scatter` or `scatter3`, you probably know that you can already plot four or five dimensional data. With `bubbleplot`, you'll be able to visualize up to 7 dimensions. Here's an example, borrowing from the examples Ameya provides with his entry.

We'll start with a shipping sample data from the Statistics and Machine Learning Toolbox.

load carsmall Origin = cellstr(Origin); % Create a table for easy overview t = table(Acceleration,Horsepower,MPG,Weight,Model_Year,Cylinders,Origin); t(1:10,:)

ans = Acceleration Horsepower MPG Weight Model_Year Cylinders Origin ____________ __________ ___ ______ __________ _________ ______ 12 130 18 3504 70 8 'USA' 11.5 165 15 3693 70 8 'USA' 11 150 18 3436 70 8 'USA' 12 150 16 3433 70 8 'USA' 10.5 140 17 3449 70 8 'USA' 10 198 15 4341 70 8 'USA' 9 220 14 4354 70 8 'USA' 8.5 215 14 4312 70 8 'USA' 10 225 14 4425 70 8 'USA' 8.5 190 15 3850 70 8 'USA'

Then create a bubbleplot using these variables. For Color and Shape, the function is able to take categorical data as well. See the help for more information.

bubbleplot(t.Acceleration, t.Horsepower, t.MPG,... % X, Y, Z t.Weight, t.Model_Year, t.Cylinders,... % Size, Color, Shape Origin,... % Text 'FontSize', 6); grid on xlabel('Acceleration') ylabel('Horsepower') zlabel('MPG')

The function is well-written, with good use of `inputParser` and plenty of error checking. People wanting to learn good programming practices should have a look. If I were to add one suggestion, I would make it so that the function output (handle objects) gets returned only when it is asked for. When I call the function without an output argument, I don't want "ans" created. Just a minor suggestion.

Thanks for a cool visualization tool, Ameya. Hope you're doing well!

**Comments**

Give this a try, and let us know what you think here or leave a comment for Ameya.

Get
the MATLAB code

Published with MATLAB® R2015a

Finally, this day has come!

You heard me right... We decided to take it to the next level with the Modeling and Simulation Challenge group where you can solve time-based simulation problems using a block-diagram editor.

To find those problems, click on the Modeling and Simulation Challenge group in the main Cody page:

When you click on the link to solve the problem, instead of the typical text box where you would write MATLAB Code, you will see a library of blocks and a canvas where you will drop and connect those blocks.

Here is what my solution for the Mass-Spring-Damper problem looks like:

As you can expect, when you click **Simulate** or **Submit**, the power of Simulink is unleashed to solve the equations you put together.

**Now it's your turn**

Put your gaming shoes on and let's solve some Cody Modeling and Simulation Challenges!

]]>Deep Learning is a very hot topic these days especially in computer vision applications and you probably see it in the news and get curious. Now the question is, how do you get started with it? Today's guest blogger, Toshi Takeuchi, gives us a quick tutorial on artificial neural networks as a starting point for your study of deep learning.

Many of us tend to learn better with a concrete example. Let me give you a quick step-by-step tutorial to get intuition using a popular MNIST handwritten digit dataset. Kaggle happens to use this very dataset in the Digit Recognizer tutorial competition. Let's use it in this example. You can download the competition dataset from "Get the Data" page:

- train.csv - training data
- test.csv - test data for submission

Load the training and test data into MATLAB, which I assume was downloaded into the current folder. The test data is used to generate your submissions.

tr = csvread('train.csv', 1, 0); % read train.csv sub = csvread('test.csv', 1, 0); % read test.csv

The first column is the label that shows the correct digit for each sample in the dataset, and each row is a sample. In the remaining columns, a row represents a 28 x 28 image of a handwritten digit, but all pixels are placed in a single row, rather than in the original rectangular form. To visualize the digits, we need to reshape the rows into 28 x 28 matrices. You can use reshape for that, except that we need to transpose the data, because `reshape` operates by column-wise rather than row-wise.

figure % plot images colormap(gray) % set to grayscale for i = 1:25 % preview first 25 samples subplot(5,5,i) % plot them in 6 x 6 grid digit = reshape(tr(i, 2:end), [28,28])'; % row = 28 x 28 image imagesc(digit) % show the image title(num2str(tr(i, 1))) % show the label end

You will be using the nprtool pattern recognition app from Neural Network Toolbox. The app expects two sets of data:

- inputs - a numeric matrix, each column representing the samples and rows the features. This is the scanned images of handwritten digits.
- targets - a numeric matrix of 0 and 1 that maps to specific labels that images represent. This is also known as a dummy variable. Neural Network Toolbox also expects labels stored in columns, rather than in rows.

The labels range from 0 to 9, but we will use '10' to represent '0' because MATLAB is indexing is 1-based.

1 --> [1; 0; 0; 0; 0; 0; 0; 0; 0; 0] 2 --> [0; 1; 0; 0; 0; 0; 0; 0; 0; 0] 3 --> [0; 0; 1; 0; 0; 0; 0; 0; 0; 0] : 0 --> [0; 0; 0; 0; 0; 0; 0; 0; 0; 1]

The dataset stores samples in rows rather than in columns, so you need to transpose it. Then you will partition the data so that you hold out 1/3 of the data for model evaluation, and you will only use 2/3 for training our artificial neural network model.

n = size(tr, 1); % number of samples in the dataset targets = tr(:,1); % 1st column is |label| targets(targets == 0) = 10; % use '10' to present '0' targetsd = dummyvar(targets); % convert label into a dummy variable inputs = tr(:,2:end); % the rest of columns are predictors inputs = inputs'; % transpose input targets = targets'; % transpose target targetsd = targetsd'; % transpose dummy variable rng(1); % for reproducibility c = cvpartition(n,'Holdout',n/3); % hold out 1/3 of the dataset Xtrain = inputs(:, training(c)); % 2/3 of the input for training Ytrain = targetsd(:, training(c)); % 2/3 of the target for training Xtest = inputs(:, test(c)); % 1/3 of the input for testing Ytest = targets(test(c)); % 1/3 of the target for testing Ytestd = targetsd(:, test(c)); % 1/3 of the dummy variable for testing

- You can start the Neural Network Start GUI by typing the command nnstart.
- You then click the Pattern Recognition Tool to open the Neural Network Pattern Recognition Tool. You can also use the command nprtool to open it directly.
- Click "Next" in the welcome screen and go to "Select Data".
- For
`inputs`, select`Xtrain`and for`targets`, select`Ytrain`. - Click "Next" and go to "Validation and Test Data". Accept the default settings and click "Next" again. This will split the data into 70-15-15 for the training, validation and testing sets.
- In the "Network Architecture", change the value for the number of hidden neurons, 100, and click "Next" again.
- In the "Train Network", click the "Train" button to start the training. When finished, click "Next". Skip "Evaluate Network" and click next.
- In "Deploy Solution", select "MATLAB Matrix-Only Function" and save the generated code. I save it as myNNfunc.m.
- If you click "Next" and go to "Save Results", you can also save the script as well as the model you just created. I saved the simple script as myNNscript.m

Here is the diagram of this artificial neural network model you created with the Pattern Recognition Tool. It has 784 input neurons, 100 hidden layer neurons, and 10 output layer neurons.

Your model learns through training the weights to produce the correct output.

`W` in the diagram stands for *weights* and `b` for *bias units*, which are part of individual neurons. Individual neurons in the hidden layer look like this - 784 inputs and corresponding weights, 1 bias unit, and 10 activation outputs.

If you look inside `myNNfunc.m`, you see variables like `IW1_1` and `x1_step1_keep` that represent the weights your artificial neural network model learned through training. Because we have 784 inputs and 100 neurons, the full layer 1 weights will be a 100 x 784 matrix. Let's visualize them. This is what our neurons are learning!

load myWeights % load the learned weights W1 =zeros(100, 28*28); % pre-allocation W1(:, x1_step1_keep) = IW1_1; % reconstruct the full matrix figure % plot images colormap(gray) % set to grayscale for i = 1:25 % preview first 25 samples subplot(5,5,i) % plot them in 6 x 6 grid digit = reshape(W1(i,:), [28,28])'; % row = 28 x 28 image imagesc(digit) % show the image end

Now you are ready to use `myNNfunc.m` to predict labels for the heldout data in `Xtest` and compare them to the actual labels in `Ytest`. That gives you a realistic predictive performance against unseen data. This is also the metric Kaggle uses to score submissions.

First, you see the actual output from the network, which shows the probability for each possible label. You simply choose the most probable label as your prediction and then compare it to the actual label. You should see 95% categorization accuracy.

Ypred = myNNfun(Xtest); % predicts probability for each label Ypred(:, 1:5) % display the first 5 columns [~, Ypred] = max(Ypred); % find the indices of max probabilities sum(Ytest == Ypred) / length(Ytest) % compare the predicted vs. actual

ans = 1.3988e-09 6.1336e-05 1.4421e-07 1.5035e-07 2.6808e-08 1.9521e-05 0.018117 3.5323e-09 2.9139e-06 0.0017353 2.2202e-07 0.00054599 0.012391 0.00049678 0.00024934 1.5338e-09 0.46156 0.00058973 4.5171e-07 0.00025153 4.5265e-08 0.11546 0.91769 2.1261e-05 0.00031076 1.1247e-08 0.25335 1.9205e-06 1.1014e-06 0.99325 2.1627e-08 0.0045572 1.733e-08 3.7744e-07 1.7282e-07 2.2329e-09 7.6692e-05 0.00011479 0.98698 1.7328e-06 1.9634e-05 0.0011708 0.069215 0.01249 0.00084255 0.99996 0.14511 1.0106e-07 2.9687e-06 0.0033565 ans = 0.95293

You probably noticed that the artificial neural network model generated from the Pattern Recognition Tool has only one hidden layer. You can build a custom model with more layers if you would like, but this simple architecture is sufficient for most common problems.

The next question you may ask is how I picked 100 for the number of hidden neurons. The general rule of thumb is to pick a number between the number of input neurons, 784 and the number of output neurons, 10, and I just picked 100 arbitrarily. That means you might do better if you try other values. Let's do this programmatically this time. `myNNscript.m` will be handy for this - you can simply adapt the script to do a parameter sweep.

sweep = [10,50:50:300]; % parameter values to test scores = zeros(length(sweep), 1); % pre-allocation models = cell(length(sweep), 1); % pre-allocation x = Xtrain; % inputs t = Ytrain; % targets trainFcn = 'trainscg'; % scaled conjugate gradient for i = 1:length(sweep) hiddenLayerSize = sweep(i); % number of hidden layer neurons net = patternnet(hiddenLayerSize); % pattern recognition network net.divideParam.trainRatio = 70/100;% 70% of data for training net.divideParam.valRatio = 15/100; % 15% of data for validation net.divideParam.testRatio = 15/100; % 15% of data for testing net = train(net, x, t); % train the network models{i} = net; % store the trained network p = net(Xtest); % predictions [~, p] = max(p); % predicted labels scores(i) = sum(Ytest == p) /... % categorization accuracy length(Ytest); end

Let's now plot how the categorization accuracy changes versus number of neurons in the hidden layer.

figure plot(sweep, scores, '.-') xlabel('number of hidden neurons') ylabel('categorization accuracy') title('Number of hidden neurons vs. accuracy')

It looks like you get the best result around 250 neurons and the best score will be around 0.96 with this basic artificial neural network model.

As you can see, you gain more accuracy if you increase the number of hidden neurons, but then the accuracy decreases at some point (your result may differ a bit due to random initialization of weights). As you increase the number of neurons, your model will be able to capture more features, but if you capture too many features, then you end up overfitting your model to the training data and it won't do well with unseen data. Let's examine the learned weights with 300 hidden neurons. You see more details, but you also see more noise.

net = models{end}; % restore the last model W1 = zeros(sweep(end), 28*28); % pre-allocation W1(:, x1_step1_keep) = net.IW{1}; % reconstruct the full matrix figure % plot images colormap(gray) % set to grayscale for i = 1:25 % preview first 25 samples subplot(5,5,i) % plot them in 6 x 6 grid digit = reshape(W1(i,:), [28,28])'; % row = 28 x 28 image imagesc(digit) % show the image end

You now have some intuition on artificial neural networks - a network automatically learns the relevant features from the inputs and generates a sparse representation that maps to the output labels. What if we use the inputs as the target values? That eliminates the need for training labels and turns this into an unsupervised learning algorithm. This is known as an autoencoder and this becomes a building block of a deep learning network. There is an excellent example of autoencoders on the Training a Deep Neural Network for Digit Classification page in the Neural Network Toolbox documentation, which also uses MNIST dataset. For more details, Stanford provides an excellent UFLDL Tutorial that also uses the same dataset and MATLAB-based starter code.

Beyond understanding the algorithms, there is also a practical question of how to generate the input data in the first place. Someone spent a lot of time to prepare the MNIST dataset to ensure uniform sizing, scaling, contrast, etc. To use the model you built from this dataset in practical applications, you have to be able to repeat the same set of processing on new data. How do you do such preparation yourself?

There is a fun video that shows you how you can solve Sudoku puzzles using a webcam that uses a different character recognition technique. Instead of static images, our colleague Teja Muppirala uses a live video feed in real time to do it and he walks you through the pre-processing steps one by one. You should definitely check it out: Solving a Sudoku Puzzle Using a Webcam.

You got 96% categorization accuracy rate by simply accepting the default settings except for the number of hidden neurons. Not bad for the first try. Since you are using a Kaggle dataset, you can now submit your result to Kaggle.

n = size(sub, 1); % num of samples sub = sub'; % transpose [~, highest] = max(scores); % highest scoring model net = models{highest}; % restore the model Ypred = net(sub); % label probabilities [~, Label] = max(Ypred); % predicted labels Label = Label'; % transpose Label Label(Label == 10) = 0; % change '10' to '0' ImageId = 1:n; ImageId = ImageId'; % image ids writetable(table(ImageId, Label), 'submission.csv');% write to csv

You can now submit the `submission.csv` on Kaggle's entry submission page.

In this example we focused on getting a high level intuition on artificial neural network using a concrete example of handwritten digit recognition. We didn’t go into details such as how the inputs weights and bias units are combined, how activation works, how you train such a network, etc. But you now know enough to use Neural Network Toolbox in MATLAB to participate in a Kaggle competition.

Get
the MATLAB code

Published with MATLAB® R2015a

Pascale works in the MathWorks Paris office and helped coordinate a student competition: the Mission on Mars Robot Challenge. The goal was to design an algorithm that would drive a rover to points of interest as quickly as possible while avoiding obstacles. And while it would have been fun to land a fleet of rovers on the Martian surface, the finals were held in Lyon in a terrain mockup.

Of course before you test in the field, the prudent thing to do is to fine-tune your design in simulation first. MathWorks provided competitors with a starter model in Simulink. As I would expect from any of my coworkers, this model is very well organized and follows many modeling best practices. All the files are managed as a Simulink Project, which makes it easy to share and get others using it quickly. Once you load the project and open the main model, this is what you discover:

The model is comprised on a subsystem that models the rover's control system, dynamics, and sensors. A separate subsystem scores the rover's performance. These subsystems served as a foundation that enabled competition teams to improve on the third subsystem: InputProcessing. In effect, this was the guidance algorithm to the rover's software.

The guidance algorithm for the rover is managed by a Stateflow chart, which contains a state machine that defines the pathing logic. If you're unfamiliar with state machines, we produced a tech talk series on the subject that's worth checking out (I hear the presenter is phenomenal). In a nutshell, state machines are a useful way of expressing a system with different modes of operation that you switch between. When you use Stateflow to design your state machine, you can visualize the mode switching as you simulate:

While this looks pretty, the state machine employs an intentionally poor algorithm. This is a sample of a simulated three-minute attempt to locate all points of interest (marked as circles). The rover doesn't know where the circles are; it has to scan for them with a camera whose range of visibility is shown with the blue trapezoid.

As you can see in the animation, there's no strategy employed to systematically sweep through the field and find points of interest. The rover randomly spins around, moves forward on occasion, and (if it's lucky) finds a circle. Even when that happens, it sometimes loses track of what it detected. The animation shows a case were it identified a point but then drove right past it. The end result of the three minutes is unimpressive. The rover wastes a lot of time covering the same ground and ultimately misses one of the six points of interest.

Think you could design a better algorithm? Well that was the point of the competition. And even though the challenge is over, this is still a fun model to play around with.

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

]]>

**What is simulation metadata?**

The `getSimulationMetadata` method of the `Simulink.SimulationOutput` object returns a `Simulink.SimulationMetadata` object. It has four properties: ** ModelInfo, TimingInfo, UserString, **and

**Timing information**

No more need to insert `tic` and `toc` calls in your model callbacks. The ` TimingInfo` field of the metadata informs you of the time it took for your model to initialize, execute, and terminate:

**Custom information**

If you want to add notes or relevant data, you can use the ** UserString** and

**Now it's your turn**

What are you going to store in the Simulink.SimulatuonOutput ** UserData** field? Let us know by leaving a comment here.

An interactive graphical experiment lets you discover the value of one of the most important numerical quantities in mathematics.

One of my favorite graphical experiments allows you to find the numerical value of $e$. The program is described in Chapter 8 of Experiments with MATLAB and is included in the toolbox available with that book. It was originally called `expgui`; today it is known as `expshow`. Here is a link to the code. I hope you are able to start up MATLAB and actually run `expshow` while you are reading this blog.

Let's say you've forgotten how to find the derivative of $a^t$ with respect to $t$. Be careful. Don't blindly follow your memory that the derivative of $t^n$ is $n t^{n-1}$ and claim that the derivative of $a^t$ is $t a^{t-1}$.

For graphical purposes, we can use an approximate derivative, replacing the tangent by a secant with a step size of $10^{-4}$. Here is the core of `expshow` and the resulting screen shot with $a = 2$

```
t = 0:1/64:2;
h = .0001;
% Compute y = a^t and its approximate derivative.
y = a.^t;
yp = (a.^(t+h) - a.^t)/h;
```

The blue line is the graph of $2^t$. At $t = 1$ it passes through $y = 2$ and at $t = 2$ it hits $y = 4$. The sienna line is the graph of the derivative of $2^t$. The important facts are that it has the same shape as the blue line and lies entirely below it.

Now take your mouse and move the blue line. It's very satisfying to be actually interact with this experiment. Push the blue line until the sienna line moves above it, to somewhere around $a = 3$.

With $a = 3$, at $t = 1$ the graph passes through $y = 3$ and at $t = 2$ it hits $y = 9$. The sienna line is now the graph of the derivative of $3^t$. The important facts are that it still has the same shape as the blue line and now lies entirely above it.

Now you know what to do.

Using the mouse, move the blue line until the two lines lie on top of each other. You have found the only function in the world that is equal to its own derivative and, in the process, discovered that, to three decimal places, the crucial value of the base is 2.718. And you did this without touching the keyboard or typing in any numbers.

Get
the MATLAB code

Published with MATLAB® R2015a

Greg's pick this week is MIPS processor in Simulink by Mikhail.

- Processor simulation in Simulink
- What's so great about this model?
- It's really a tour of Turing
- Wait, I can use a program to write a program that runs a program that I have written in my own language?
- What in the name of Thor's Hammer is HDL?
- This seems confusing...
- What's your take on processor simulation?

Ever wonder how a microprocessor runs a program to perform computations? Better yet, how would you go about building and implementing such a design?

Mikhail's example is simple and based on one published in a book. But it is a well laid out description of the computational process model using Simulink® that supports both the simulation of the program on your desktop machine, as well as the generation of HDL-code to deploy the design onto an FPGA.

To me this is just super cool. Perhaps that says more about me than this File Exchange entry. I'm sure there's some marketing message I should be spouting about how important simulation is to design. But really, what grabbed my attention is

I can follow what was going on and watch the program execute in simulation The model is neat and clean I can see myself writing a little parser or compiler to program it The resulting design can actually be implemented in hardwareUltimately this means *I could create my own processor and programming language.*

I very much like understanding how things work. I have thought a great deal about processors, but perhaps hadn't considered how they actually work in too much detail. After all, my background is mechanical engineering and not computer science.

However... back in university, I took a philosophy course on artificial intelligence, and we spent a good deal of time developing and discussing Turing machines. (These have come back into the spotlight recently with the release of the movie about Alan Turing).

The basic premise behind a Turing machine is that an instruction can read a datum, perform some action with that datum, and write a resulting datum into some location that can store the information.

That is precisely what this model describes. I suppose that's good, because a microprocessor is a Turing machine (well, a finite version of a Turing machine)

It seems like there's some sort of circular dependency here, but the short answer is yes, you can do that. In fact Alan Turing basically proved this notion with the Universal Turing Machine.

You can develop an algorithm that represents fundamental operations. In this case the algorithm can be implemented on hardware by generating HDL from the Simulink model.

Hardware description language (HDL) is a textual language to determine how a piece of electronics will operate by describing the structure of that device.

In terms of silicon devices, it describes data pathways through which bits will flow and be stored. Field programmable gate arrays (FPGAs) are unique silicon fabrications that can be "reconfigured" to a different set of data pathways (as opposed to your computer's processor which has a fixed set of data pathways).

In fact, processor designers and developers use HDL to create and implement processor designs. Often HDL appears as a very low-level language because you are often dealing at the level of individual bits.

It's clear that these concepts are quite abstract, which is why I found the MIPS model so interesting. I could begin to deconstruct how the processor will ultimately parse and respond to an instruction.

I used the fact this algorithm is in Simulink to interrogate different aspects of the processor language design such as,

- What is the current processor instruction at any given time in the program?

I even fantasized about writing my own compiler for this using MATLAB regular expressions or perhaps ANTLR, but writing these blog posts is hard enough!

Have you ever built your own processor? Would you be interested in simulating how a processor works? Do you believe this approach in Simulink is scalable?

Leave your comments here.

Get
the MATLAB code

Published with MATLAB® 8.5

Starting on July 20, 2014 and continuing through this month, the spacecraft acquired almost 800 pictures using three different scientific cameras. The pictures were transmitted to Earth for processing by two independent navigation teams. The PNAV (Project Navigation) team at KinetX had primary navigation responsibility. The PNAV team used MATLAB software. The INAV (Independent Navigation) team at the Jet Propulsion Laboratory used software written in Fortran and C. The PNAV and INAV navigation systems were tested and compared with each other during the mission’s Jupiter flyby and were found to be in agreement.

Image processing and optimization methods were used to determine the location of Pluto and its moons with high precision as soon as possible in the mission. Accurate course corrections early saves much fuel later.

Several factors posed challenges to determining the exact trajectory to Pluto:

- The low resolution of pictures while the spacecraft was still far away limited the precision of the computations.
- The size, shape, and surface brightness variation of Pluto and its moons were not well known. These unknowns had to be estimated and compensated for.
- The effects of spacecraft motion, pointing drift while the camera was collecting light, and the travel time of light all had to be included in the computation.

The PNAV team at KinetX used the MATLAB version of NASA’s SPICE Toolkit, which is used in multiple planetary missions to determine:

- where the spacecraft is located
- how the spacecraft and its instruments are oriented
- the location, size, shape, and orientation of the target being observed
- the events on the spacecraft or ground that might affect scientific observations

Determining the best path to Pluto required three steps:

- Find the center of the stars in the field of view and match them to star catalogs.
- Find the spacecraft attitude, or direction that the spacecraft is pointing.
- Find the center of the target bodies (Pluto and its moons).

Like the overall approach of using two independent navigation programs for the program, each of these analysis steps was solved using two different algorithms, and the results were cross-checked.

In each case, one of the algorithms was based on image processing methods, and the other was based on a nonlinear least squares estimator.

Finding the center of stars was accomplished using matched filtering using previously acquired images of stars in the star catalog.

Spacecraft attitude was determined using an image registration technique that compensated for x- and y-shifts in the image plane as well as rotation.

Target body center-finding used cross-correlation with a model of the target body. As the mission proceeded, the target body models were refined based on new data about the size, shape, and surface brightness variations of Pluto and its moons.

*Image credit: Owen, Dumont, and Jackman, “Optical Navigation Preparations for New Horizons Pluto Flyby”*

*Image credit: APL / New Horizons GeoViz*

*Image credit: NASA/Johns Hopkins University Applied Physics Laboratory/Southwest Research Institute*

“Optical Navigation Preparations for New Horizons Pluto Flyby”

“New Horizons to take new photos of Pluto and Charon, beginning optical navigation campaign”

“An Overview of SPICE: NASA’s Ancillary Data System for Planetary Missions”

]]>**The Problem**

Let's say I have a simple model with a control loop:

If the plant model is direct feedthrough, this will result in an algebraic loop. While Simulink can solve the algebraic loop most of the time, it usually slows down the simulation, and when the solve fails to converge it can lead to errors like this:

**Breaking the loop with a Memory Block**

To break the algebraic loop, you need to insert in the loop a nondirect feedthrough block. The first thing most users think about is a Unit Delay or Memory block.

If the blocks in the algebraic loop have a discrete sample time, inserting a Unit Delay is usually the best solution. Of course this will change the dynamic of the system, this is something you need to evaluate and see if this is acceptable for your application.

If the blocks in the loop have a continuous sample time, what many users try is inserting a Memory block. The Memory block is similar to the Unit Delay block in a sense that it delays its input by one time step, however it works with variable-step signals. Let's see what it does for our model.

At least, now the model simulates to completion and we can look at the results:

However when simulating the model, we quickly notice that it simulates very slowly. If I log data from the model, I can see that it takes more than 500,000 steps to simulate this model for two seconds!

If we plot the steps taken by the variable-step solver, we can see that after the step input, the solver begins to take steps of about 1e-6 seconds and gets stuck there.

Why is it happening? This is because the output of the Memory block is not continuous, and it is driving a block with continuous states, the State-Space block. Every time the output of the Memory block changes, the solver needs to reset, forcing the small step size that we observe. We know that this situation is problematic and we have a Model Advisor check for that: Check for non-continuous signals driving derivative ports

**The Solution: Breaking the Loop using a Transfer Function block**

As suggested by the Model Advisor, the recommended way to break this algebraic loop is to use a continuous block. The one I typically prefer is a first order Transfer Function. Like the Memory block, this will introduce a new dynamic in the system. The trick is to make the time constant of the Transfer Function small enough to not affect the dynamics of the system significantly. In this case, I used 1e-6.

With this change, the model gives similar results, but the simulation completes almost instantly, taking only 633 time steps:

**Now it's your turn**

If you have experiences or suggestions on how to handle algebraic loops, let us know by leaving a comment here.

***** Important Update *****

After publishing this post, a few users contacted me mentioning that some Simulink demos use a Memory block to break algebraic loops. I consequently decided to add this update to highlight the fact that breaking an algebraic loop with a Memory block is problematic only when the loop is continuous. Let's look at a few of those examples and explain why, because the loops are not continuous, it is ok to break them with a Memory block.

sldemo_clutch: In this model, the following pattern is used in the clutch logic:

You can notice that I enabled the port data type and sample time displays to highlight that this loop take a fixed-in-minor-steps sample time, and the data type of the signals involved is boolean. This subsystem implements a discrete combinatorial logic deciding if the clutch should be locked or not depending on two inputs and it's previous state. Since the loop is discrete, the Memory block is the way to go.

sldemo_bounce: In this model, we can see that an algebraic loop is broken by a Memory block:

At first look, this loop has a continuous sample time and is of data type double. So why am I not considering it continuous? Because of when the loop is actually active. Let's look at the logic here. First, we need to note that the Integrator is configured to reset dx/dt when x reaches saturation:

As soon as the Integrator enters the saturation, we want to apply a new velocity that is 80% of the velocity when it entered the saturation, but in opposite direction, to take it out of the saturation. This means that the output of the loop is not used continuously by the Integrator. It is used only for one time step, at a discontinuity, when entering the saturation triggers a zero-crossing event.

I hope those clarification makes it clearer that when I recommend breaking algebraic loops with transfer functions, I am talking about continuous algebraic loops, where (if a Memory block was used) the output of the Memory block would drive a derivative port as can be detected by the Model Advisor check mentioned above.

Jiro's picks this week are Align axes labels in 3D plot by Matthew Arthington and Tools for Axis Label Alignment in 3D Plot by Ligong Han.

When you create a plot, you probably notice a bunch of buttons in the toolbar.

These buttons have been around for a very long time, so you probably have gotten used to these powerful features. They allow you to quickly explore your data in different ways, i.e. by zooming, panning, and rotating.

But did you know that you could combine these interactive tasks with programmatic tasks? For example, you could make MATLAB tell you the current X and Y limits every time you zoom.

plot(rand(10,1)) h = zoom; h.ActionPostCallback = @(o,e) disp(axis);

This topic was covered in a different blog "MATLAB Spoken Here" in this post. As you can see, this is a feature that has been around for a long time.

The two File Exchange submissions by Matthew and Ligong are perfect for combining with this feature for 3D rotation. They allow for automatic alignment of axes labels when you rotate the figures. Matthew's submission came first. It works well and has a lot of good reviews. One limitation was that it only worked with equal aspect ratio and orthographic projection. Ligong was inspired by this and created a version without this limitation.

**Comments**

Give these a try and let us know what you think here or leave a comment for Matthew or Ligong.

Get
the MATLAB code

Published with MATLAB® R2015a

Here I review a few methods for building large MATLAB tables that require scalar operations like this, and I compare the relative speed of these methods.

See documentation on MATLAB tables for more information.

]]>Do you know what that is? Try it now. Go to my recent post, “Displaying a color gamut surface,” and scroll down near the bottom. Just above the comment section you’ll see this:

Go ahead and click on that link. You should see something like this in your browser:

If you have used this little blog feature, I would appreciate you leaving me a comment about it.

Thanks!

]]>I'm happy to welcome back Damian Sheehy as guest blogger. Last time Damian wrote about how Natural Neighbor interpolation addresses FAQs in scattered data interpolation. In this blog he will answer a FAQ on adaptively editing a Delaunay triangulation.

A technical support question that occasionally crops up asks about the best and most efficient way to construct a Delaunay triangulation by adding points to an existing triangulation. We call this operation incremental editing as the goal is to add or remove points in an incremental manner as opposed to recreating from scratch; for example calling `delaunay` multiple times. The documentation for the `delaunayTriangulation` class provides examples that show the syntax that allows you to edit a Delaunay triangulation by adding or removing points. Let's look at a simple 2-D example to highlight the concept.

In this example we will load `trimesh2d.mat` which ships with MATLAB. The file contains points with coordinates (x, y) and triangulation edge constraints that define the boundary of a domain. Let's triangulate the data and take a look at that.

load trimesh2d dt = delaunayTriangulation(x,y,Constraints); tin = find(dt.isInterior()); triplot(dt(tin,:),dt.Points(:,1),dt.Points(:,2)) axis equal

Suppose we want to add the circumcenter points to this triangulation. The circumcenter of a triangle is the center of a circumscribed circle that passes through the vertices of the triangle. The `delaunayTriangulation` class provides a method to compute them - `delaunayTriangulation/circumcenter`. Compute them using this method and add them to the plot. Note, some triangles may share the same circumcenter, so I will call the `uniquetol` function to eliminate the near-duplicates.

cc = dt.circumcenter(tin); cc = uniquetol(cc, 'ByRows',true); hold on plot(cc(:,1),cc(:,2),'.r') hold off

Now, I will add the circumcenter points to the triangulation and plot the result in a new figure.

```
numcc = size(cc,1);
dt.Points(end+(1:numcc),:) = cc;
tin = find(dt.isInterior());
figure
triplot(dt(tin,:),dt.Points(:,1),dt.Points(:,2))
axis equal
```

So that's basically an incremental change we made to the triangulation. When prototyping applications like this at the command line we may find adding a few points to a large triangulation can take almost as long as creating the full triangulation of all the points. Why is that, surely the operation should be more efficient?

There are some practical applications that rely on this expected level of efficiency. For example, an incremental algorithm for 2-D Mesh Generation using Ruppert's Algorithm. Delaunay-based algorithms for reconstructing surfaces from point clouds may also be incremental; triangulating an initial set of points and adaptively adding more points to recover the surface. In fact the additive points in these algorithms are often circumcenters and that's why I chose them in the example. But the question remains, shouldn't an incremental addition of a few points to a large triangulation be more efficient than a complete triangulation of all points. Absolutely, this behavior is honored in the scenario where algorithms written in MATLAB are designed to run most efficiently. If we write the algorithm in a function in a file, the incremental change will be performed efficiently. If we prototype at the command line the performance we get may underperform the efficiency we get from the function-in-a-file format.

Let's run a little example to test this; the following code creates a 3-D Delaunay triangulation of a half-million points and subsequently adds 40K points in 4 incremental updates. Here's the output when the code runs at the command line:

`timeToCreateDelaunay =`

9.7799

`timeToIncrementallyAddPoints =`

10.1267

When I put the code in a function in a file, execution gives the following:

function DelaunayIncrementalTest() tic; dt = delaunayTriangulation(rand(500000,3)); timeToCreateDelaunay = toc tic; dt.Points(end+(1:10000),:) = rand(10000,3); dt.Points(end+(1:10000),:) = rand(10000,3); dt.Points(end+(1:10000),:) = rand(10000,3); dt.Points(end+(1:10000),:) = rand(10000,3); timeToIncrementallyAddPoints = toc end

DelaunayIncrementalTest()

timeToCreateDelaunay = 9.9275 timeToIncrementallyAddPoints = 1.1396

Why the significant performance improvement? MATLAB can execute code more efficiently when it is in the function-in-file format. Loren has a past blog on In-place Operations that highlights the behavior that improves the efficiency here. So to gauge the performance of your MATLAB code it's good to structure your code into functions that reside in files. Then run the performance analyzer and eliminate the bottlenecks.

Does your work involve triangulations or geometric computing? Does your application area require geometric tools or features that are not well supported in core MATLAB? Close the loop and share your experience here; hearing what users' need is the compass that charts our feature enhancements!

Get
the MATLAB code

Published with MATLAB® R2015a

I have just returned from two meetings in Europe, the 26th Biennial Conference on Numerical Analysis at the University of Strathclyde in Glasgow, Scotland, and Sparse Days III in Saint-Girons, France.

The Biennial Conference on Numerical Analysis in Scotland has a long history. The first two meetings were held at the University of St. Andrews in 1965 and 1967. The meetings moved to the University of Dundee in 1969. They were held there every two years, under the leadership of Ron Mitchell, until his retirement in 2007. They moved to the University of Strathclyde in 2009, under the leadership of Alison Ramage. This meeting was the 26th and the 50th anniversary. Although this is a premier conference series in numerical analysis, it was the first time I participated.

The University of Strathclyde is a public research university located in Glasgow, Scotland. The university's name comes from the "Valley of the River Clyde". The school is Scotland's foremost institution of science and engineering.

The conference was held June 23rd through 26th. I was one of a dozen invited speakers -- five from the US, four from continental Europe, and three from the UK.

Almost all of the participants gave half-hour talks. There were 168 talks organized in seven parallel sessions. Many of the talks were in one of twelve minisymposia. Sample minisymposia titles include "Stable and accurate discretisations for convection-dominated problems", "City analytics", "Chebfun: new developments, cool applications and on the horizon", and "Numerical linear algebra for optimsation and data assimilation".

Two of the invited lectures are given special emphasis to honor British numerical analysts. The A. R. Mitchell Lecture honors the University of Dundee's Ron Mitchell who was the dominant force in this conference for forty years.

This year the lecture was given by Mike Giles of the University of Oxford on "Multilevel Monte Carlo Methods". I knew nothing about this subject before his talk and I learned a great deal. His abstract provides a link to a web page featuring a survey paper and MATLAB codes, <http://people.maths.ox.ac.uk/gilesm/acta>

The Fletcher-Powell Lecture honors two mathematicians known for their work in optimization algorithms, including the Davidon-Fletcher-Powell, DFP, formula. Mike Powell, from Cambridge University, had passed away in April. Roger Fletcher, from the University of Dundee, attended the lecture.

This year the lecture was given by my long-time friend Mike Saunders from Stanford. He talked about "Experiments with linear and nonlinear optimization using Quad precision." He has used the GFortran compiler to build a version of his MINOS optimization software with the REAL*16 floating point datatype. He tackled flux balance analysis models of metabolic networks that involve coefficients ranging over 15 or 16 orders of magnitude. Ordinary double precision, that is Fortran's REAL*8, cannot do the job.

The history of the Sparse Days in St Girons conference is itself sparse. There have been just two previous conferences, one in 1994 and one in 2003. The conference is part of a program on "High Performance Linear and Nonlinear Methods for Large Scale Applications" organized by CIMI, the French Centre International de Mathematique et d'Informatique in Toulouse.

Organizers of Sparse Days in St Girons III included my buddies Iain Duff, who splits his time between Rutherford Appleton Laboratory in the UK and the Parallel Algorithms Group at CERFACS in Toulouse; Jack Dongarra, from the University of Tennessee; and Jim Demmel, from the University of California, Berkeley.

St Girons is a picturesque town in the French Pyrenees, not far from the border with Spain. The Tour de France sometimes passes through here. We met in the town's only cinema, which is air conditioned. This turned out to be fortunate planning because a heat wave hit that week, June 28th through July 2nd, with temperatures in the mid 30s C, which is mid 90s F.

The format for this conference was quite different from the one the previous week in Scotland. There were only about half as many attendees, around 100. There were no invited talks and no parallel sessions. Anyone who wanted to give a talk gave one. There were 42 talks, most of them half an hour.

Tim Davis, MathWorks consultant on sparse matrices, who recently moved to Texas A&M, talked about "Sparse SVD, and a GPU-accelerated sparse QR".

Bora Ucar of Ecoles Normales Superieures in Lyon talked about "Two approximation algorithms for bipartite matching on multicore architectures." The analysis of his algorithm happened to involve my old friend, the Lambert W function.

Joost Rommes, from Mentor Graphics in the Netherlands, talked about "Challenges in numerical simulation of electrical networks."

Alex Pothen, from Purdue, talked about "The Virtual Scalpel: real-time matrix computations and finite elements for surgery."

Dan Sorensen, from Rice, talked about "A DEIM induced CUR factorization." This title means that C is some columns of a matrix A, R is some rows of A, and DEIM is an algorithm that finds U so that the product CUR is a good low approximation to A.

Martin Gander, from Universite de Geneve, talked about "Five decades of time parallel integration." I was particularly interested in his talk because I knew something about the subject many years ago, but was not familiar with recent developments. At first, it might seem impossible to parallelize a computation that involves stepping along in time. But Martin gave an overview of four different approaches to how this can be done.

There is another Cleve in this business. Cleve Ashcraft of LSTC, Livermore Software Technology Corporation, is actually a "grandstudent". His PhD advisor at Yale was Stan Eisenstat, who was my PhD advisee when I was a visiting professor at Stanford. At Sparse Days, Cleve talked about "Separability, partitions and coverings.". His talk didn't mention the applications he is deeply involved in at LSTC, which are the developers, among other things, of the crash analysis software LS-DYNA. If you want to see a hint of some serious use of dynamic finite element calculations, take at the LSTC web pages.

This is the group photo from Sparse Days in St Girons III. There are two Cleves. Cleve Ashcraft is the guy with the green shirt and grey beard in the front row. Can you find the other one? (Thanks to Pierre-Henri Cros for the photo, and for handling local arrangements.)

Get
the MATLAB code

Published with MATLAB® R2015a

*Image credit: NASA/JHUAPL/SWRI*

The New Horizons spacecraft reaches its closest approach to Pluto tomorrow!

Check out NASA TV for the countdown starting at 7:30 AM (America/New York), Tuesday, July 14.

Follow @NASANewHorizons on Twitter.

]]>Sean's pick this week is Auto Flush by the Techsource Technical Team.

My pick this week is something that many of us could find useful and have likely thought about at some point during the standard workday: a simple design for the controller behind an automatically flushing toilet.

Techsource's Technical Team has put together a collection of simple Simulink models to showcase getting started with Simulink and an Arduino Board, a low cost embedded target.

This particular one caught my attention for two reasons. First, the main model is called *'UltraSonic_Pee.slx'*. If you want to grab someone's attention, well, a model name like this is a good way to do it. Second, this simple model
is extensible so that I can include my own sensor input and simulate different flushing scenarios. Let's see how this is
done.

Here's the original model:

For my Arduino board, I only have the LEDs and a simple DC motor. However, I have a Lego NXT that has an ultrasonic distance sensor and a DC motor that I can use. I don't know if the DC motor is actually powerful enough to pull down the handle on the urinal, but I'm going to try.

The first thing I am going to do is replace the two Arduino Blocks with a Lego Motor block from the Lego Support Package. Support packages are additional functionality you can freely download to allow MATLAB and Simulink to take advantage of your hardware.

One of the beautiful things about Simulink is the ability to have Variant Subsystems. These allow you to either have different fidelity models or to substitute in different components. My second step is going to be to make the Signal Builder block into a variant subsystem so I can have my Lego distance sensor as an input instead. This is really important for the modeling aspect. I'm fortunate enough right now to have access to the hardware. But if I did not, I could still develop and test my algorithms by simulating the hardware with other signals that I could build or historical data from other tests.

Next I'm going to simulate this in External Mode. This means that the model will be running in Simulink but be grabbing data from the sensors and controlling the motor. By running in External Mode we can use the full power of the MATLAB and Simulink platforms to analyze the results of the model and see it running in semi-real-time. By doing this, I discovered that I needed a couple of gains on the signals to control how close one needs to get to trip the sensor.

Once the model has been tested, we can then embed the controller onto the Lego. One button click in Simulink generates the equivalent C-code, compiles it, and moves it onto the Lego.

Now we're ready to put this to the test!

The looks we got from fellow restroom patrons and the cleaning staff were pretty priceless.

Next I got my friend Adam, who is relocating from Massachusetts to California shortly, to test it out. The control algorithm that Techsource used is to have a preflush and a postflush - not the most environmentally friendly of algorithms but we waste no water in simulation. Since this is no longer simulation, I figured Adam should enjoy this luxury while still being on the east coast.

Unfortunately, the Lego motor was not strong enough to pull the handle down and instead lifted the whole robot up. But hey, not too shabby for an hour's worth of work.

Give it a try and let us know what you think here or leave a comment for the Techsource Technical Team.

Get
the MATLAB code

Published with MATLAB® R2015a

I’ve found these issues can often be addressed with parallel and GPU computing, often with minimal changes to my code. I think this video and others in the video tutorial series are a good introduction to the area. You can also download all the code examples to follow along.

**Algebraic Constraint Violated**

Did you ever run into the following warning?

Based on my experience, if you get this warning as you are building your model, you better address it immediately. As your model grows larger, this kind of problem can become more difficult to debug, and is likely to return when your model gets more complex.

**What does this warning mean?**

In many cases, Simulink models are made of Ordinary Differential Equations (ODE). For those, the model computes the derivative of all states, and the solver, for example ode45, integrates each of them almost independently, ensuring they all respect the tolerances specified in the model configuration.

For some models, especially Simscape models, the equations are not a set of ODE, but a set of Differential Algebraic Equations (DAE).

What does this means? This means that in addition to respecting the tolerances specified in the model configuration, the states must also respect algebraic costraints linking them together.

An example DAE can be found in the MATLAB documentation, and looks like:

As you can see, we do not have an equation that computes the derivative of `y3` as in a ODE. Instead, we have an equation that algebraically links `y3` to `y1` and `y2`.

If the algebraic constraints between the equations of the system are very complex and change very fast, it's possible that Simulink will fail to respect the constraint even when taking the smallest possible step size.

**Finding the root cause**

Here is what I do when that happens. To begin, launch the Simulink command-line Debugger:

Then I set a time breakpoint just before the warning and start the simulation until it reaches that point:

At this point, I want to enable the maximum level of solver tracing, set a breakpoint on solver errors, and continue:

When the breakpoint is hit, typically I see a large amount of failed steps with the comment "Newton iteration failed to converge". The number next to "Ix", is the index of the state failing to respect the algebraic constraint.

You can get its name using the `states` command:

The complex algebraic constraint is probably close to this state. In Hydraulic and Pneumatic domains, this is often caused by what we call a "dry node". For example, in the arrangement below, the flow going through the Check Valve and the Directional Valve is algebraically linked. The flow going through the Check Valve must also go through the Directional Valve, there is nowhere else to go. If the Check Valve cracks, or the Directional Valve command change quickly, this algebraic constraint changes over time.

To workaround this situation, you can reduce the solver tolerances to force it to take smaller steps and better capture the changes in the algebraic constraints. Another option is to break the algebraic constraint. In the Hydraulic domain, this is usually done by inserting a Constant Volume Chamber:

**Conclusion**

I know this method is manual, and not all people are comfortable with the command line debugger, but this has helped me a lot in the past, so I thought I should share. We are working on improving the process for debugging this kind of problem, and I look forward to sharing that with you in the future.

Try this approach and let me know what you think by leaving a comment here.

]]>I hope to see many of you in Québec City!

]]>Sean's pick this week is Testing of Safety Critical Control Systems by Yogananda Jeppu.

With the loss of the SpaceX Falcon 9 last week it seems like an appropriate time to read through Yogananda's compilation of many possible failures and mitigations in a control system. The cause of the failure for the SpaceX Falcon is not yet publicly known so it will be interesting to hear what their engineers discover.

In this document, Yogananda has covered a wide variety of accidents and their causes, possible failures in specific parts of a control system, how to identify and circumvent these potential failures and tips from his experiences.

The most recent update to this slide deck, includes a bit on Simulink Design Verifier and formal methods. I'll be curious to see if he extends it to include examples using Polyspace Code Prover which uses formal methods and static analysis to prove the lack (or presence!) of run-time errors in C/C++ code. For example:

Do you design, test, or research failures in safety critical control systems? If so, are there any other insights that you would like to share?

Give it a read and let us know what you think here or leave a comment for Yogananda.

Get
the MATLAB code

Published with MATLAB® R2015b

I'm happy to welcome Damian Sheehy as this week's guest blogger. Damian works on the development of geometry-related features at MathWorks. He will provide answers to two frequently asked questions; one on scattered data interpolation that he will cover in this blog and the other on Delaunay triangulation that he will cover in the next. Over to you, Damian...

I occasionally get an email from customer support with a title similar to this one: "griddata gives different results. . . ". The support engineers are great, they really know how to choose a good subject line that will get a developer's attention and get a response back to the customer quickly. The subject line could equally well cite `scatteredInterpolant` as it shares the same underlying code as `griddata`. Before I open the email I have a strong suspicion about the cause of the difference. If the first line opens with this: "A customer upgraded to MATLAB R20** and `griddata` gives different results to the previous . . .", then I'm fairly confident what the problem is. I look at the customer's dataset, perform a couple of computations, create a plot, and bingo!, I have a canned response and recommendation at the ready and I turn around the question in a matter of minutes.

So why should `griddata` or `scatteredInterpolant` give different answers after upgrading MATLAB? What has MathWorks done to address this problem? Are there issues with scattered data interpolation that users should be aware of? Yes, there are some subtle behaviors associated with the Nearest Neighbor and Linear interpolation methods for scattered data interpolation. These problems present themselves in specific datasets and the effects may show up as numerical differences after a MATLAB upgrade. I will explain the origin of these problems and the options you have as a user to avoid them altogether. I will also highlight what MathWorks has done to address the problems.

First, let's take a look at the behavior and the data that triggers the problem. To demonstrate, I will choose a simple data set where we have four points at the corners of a square. Each sample point has a different value associated with it and our goal is to compute an interpolated value at some query point within the square. Here's a diagram:

Px = [0; 1; 1; 0]; Py = [0; 0; 1; 1]; V = [10; 1000; 50; 100]; plot(Px, Py, 'or') hold on text(Px+0.02, Py+0.02, {'P1 (10)', 'P2 (1000)', 'P3 (50)', 'P4 (100)'}) pq = [0.5 0.5]; plot(pq(:,1), pq(:,2), '*b') hold off axis equal

First, let's consider the Nearest Neighbor interpolation method. For any query point within the square, the interpolated value is the value associated with the nearest neighbor. The figure shown above illustrates the configuration and sample values in parenthesis. We can see an ambiguity arises when the query point lies at the center of the square. There are four possible interpolation solutions and based on the definition of the method any of the four values is a valid solution. Ideally, we would like to have the same result, no matter what computer MATLAB is running on and no matter what version.

This type of problem can also arise with the Linear interpolation method. To perform linear interpolation, the scattered dataset is first triangulated using a Delaunay triangulation. The interpolated value at a query point is then derived from the values of the vertices of the triangle that enclose the point. But a Delaunay triangulation of this dataset is not unique, the illustration below shows two valid configurations.

subplot(1,2,1); Px = [0; 1; 1; 0; 0; 1]; Py = [0; 0; 1; 1; 0; 1]; pq = [0.5 0.25]; plot(Px, Py, '-b') hold on plot(pq(:,1), pq(:,2),'or') hold off axis equal subplot(1,2,2); Px = [0; 0; 1; 1; 0; 1]; Py = [1; 0; 0; 1; 1; 0]; plot(Px, Py, '-b') hold on plot(pq(:,1), pq(:,2),'or') hold off axis equal

The interpolated result is different in each scenario. The code to demonstrate this is given in the frame below, I have perturbed one data point to flip the diagonal of the triangulation and illustrate the effect.

P = [0 0; 1 0; 1 1; eps (1-eps);]; V = [10; 1000; 50; 100]; pq = [0.5 0.25]; F = scatteredInterpolant(P,V); linearVq = F(pq) P = [0 0; 1 0; 1 1; -eps (1+eps);]; F.Points = P; linearVq = F(pq)

linearVq = 527.5 linearVq = 267.5

The interpolated value at the query point, `linearVq`, is sensitive to how the triangulation edge is created in the tiebreak case. This tiebreak is called a degeneracy and the problem arises in Delaunay triangulations when we have 4 or more cocircular points in 2-D or 5 or more cospherical points in 3-D. Now observe the behavior when we choose the Natural Neighbor interpolation method. This method is also based on an underlying Delaunay triangulation, but it produces the same result in the presence of a diagonal edge swap. Here it is:

```
P = [0 0; 1 0; 1 1; eps (1-eps);];
F.Points = P;
F.Method = 'natural';
naturalVq = F(pq)
P = [0 0; 1 0; 1 1; -eps (1+eps);];
F.Points = P;
naturalVq = F(pq)
```

naturalVq = 397.5 naturalVq = 397.5

Natural Neighbor is also a smoother interpolating function, so it makes a lot of sense to favor Natural Neighbor over Linear. Well that begs the question: shouldn’t the Natural Neighbor method be the default? This would be our preferred choice of default, but this method was added to MATLAB long after the `griddata` Linear method was introduced. Changing a default would have a significant impact, so the change would be more disruptive than the problem it addresses. The Natural Neighbor method is also more computationally expensive, so for large datasets Linear may be preferred for performance reasons.

The example we just reviewed highlights the nature of the problem and gives you a more stable alternative to avoid potential differences from scattered data interpolation after you upgrade MATLAB. For our part here in development, we have long recognized these types of issues create problems for users and we have adopted better underlying algorithms to address them. If you are a long-time user of MATLAB and the `griddata` function, you may recall more annoying past behavior. Prior to R2009a, repeated calls to the function using the now redundant `{'QJ'}` Qhull option gave potentially different results on each call. That problem was resolved in R2009a along with the introduction of Natural Neighbor interpolation for its stability and superior interpolation properties. Since then, improvements in the underlying triangulation algorithms have led to stable and consistent results across all platforms, first for 2-D and then for 3-D. Unfortunately, introducing these improvements meant a potential change in behavior of the upgrade for specific datasets, but that's like having to break eggs to make cake. Looking ahead, the behavior of Delaunay triangulation based functions should be much more stable. Of course, fundamental changes to the coded algorithm are likely to trigger changes like the ones we've seen, though the scope of the problem has been reduced substantially.

I haven't received many support escalations questions related to Natural Neighbor interpolation since it shipped in R2009a. I often wonder if this is because it is working well for users or if users may not be using it. The similar naming to Nearest Neighbor interpolation may cause some to misunderstand the method. How about you, have you used Natural Neighbor interpolation and how has it worked out? Let me know here.

Get
the MATLAB code

Published with MATLAB® R2015a

Augustin (Austin) Dubrulle deserves to be better known in the numerical linear algebra community. His version of the implicit QR algorithm for computing the eigenvalues of a symmetric tridiagonal matrix that was published in a half-page paper in *Numerische Mathematik* in 1970 is faster than Wilkinson's version published earlier. It is still a core algorithm in MATLAB today.

"QR" and "QL" are right- and left-handed, or forward and backward, versions of the same algorithm. We are used to thinking of factoring a matrix into an orthogonal factor, Q, and an upper or right triangular factor, R. This leads to QR algorithms. But for reasons having to do with starting a loop at $1$ rather than $n-1$, the authors of the Handbook Series on Linear Algebra decided to use left triangular and QL algorithms.

Augustin (Austin) Dubrulle is the only guy that I know who improved on an algorithm of Jim Wilkinson. Austin was born in France. He began his career in the early 1960's with IBM at their *Institut de Calcul Scientifique* in Paris. In 1968 he transferred to an IBM center in Houston, and began work on SSP, IBM's Scientific Subroutine Package.

After studying papers and books by Wilkinson and Parlett, Austin derived his own version of the implicit QR algorithm for computing the eigenvalues of a symmetric tridiagonal matrix. He was careful about the use of common subexpressions to save operations. This is especially important when computing eigenvectors because the transformations are applied to an accompanying full matrix.

Austin told me that he wrote up his result in the style of the Linear Algebra series in *Numerische Mathematik* and asked IBM for permission to submit the paper for publication. He was told by company lawyers that, in view of the Justice Department antitrust suit, IBM could not give away software, and that the ALGOL included in the paper was software. Austin tried to make the case that ALGOL was primarily a means of human communication and secondarily a programming language of little practical use, but that argument did not fly.

When Martin and Wilkinson published their implicit algorithm, Austin was encouraged to realize that his version had fewer operations in the inner loop and so would be faster. He arranged to meet Wilkinson at the 1969 University of Michigan numerical analysis summer school. Jim encouraged Austin to submit a note with just the algorithm, no provocative software, to the series.

Here is a recreation of Dubrulle's half-page 1970 paper in *Numerische Mathematik*. This is the entire paper. You can get the original here or purchase it here.

(You can get the Martin and Wilkinson contribution here or purchase it here.)

'********************************************************************************************'

The following is a formulation of the *QL* algorithm with implicit shift which requires fewer operations than the explicit and implicit algorithms described in [1] and [2].

Let $d_i^{(s)}$ $(i=1,...,n)$ and $e_i^{(s)}$ $(i=1,...,n-1)$ be the diagonal and codiagonal elements of the matrix at the $s$ -th iteration and $k_s$ the shift. Then the $(s+1)$ -st iteration can be expressed as follows.

*Acknowledgements*. The author wishes to thank Dr. J.H. Wilkinson for his helpful comments

References

1. Bowdler, H., Martin, R.S., Reinsch, C., Wilkinson, J.H.: The QR and QL algorithms for symmetric matrices. Numerische Mathematik 11, 293-306 (1968).

2. Martin, R. S., Wilkinson, J.H.: The implicit QL algorithm. Numerische Mathematik 12, 377-383 (1968).

A. Dubrulle IBM Corporation Industry Development-Scientific Applications 6900 Fannin Houston, Texas 77025, U.S.A.

'********************************************************************************************'

Dubrulle's algorithm still needs a little work. The computation of the next $e_{i+1}$ as the square root of the sum of squares is dangerous. There is potential for unnecessary underflow or overflow. I discussed this in a post almost two years ago on Pythagorean Addition. I don't suggest actually using the *pythag* algorithm. Instead, do something like this.

if abs(p) >= abs(q) c = q/p; t = sqrt(1+c^2); e(i+1) = t*p; s = 1/t; c = c*s; else s = p/q; t = sqrt(1+s^2); e(i+1) = t*q; c = 1/t; s = s*c; end

When the papers from *Numerische Mathematik* were collected to form the 1971 *Handbook for Automatic Computation, Volume II, Linear Algebra*, almost all of them where reprinted without change. But, despite what its footnote says, Contribution II/4, *The Implicit QL Algorithm*, never appeared in the journal. The paper is the merger of Dubrulle's note and the Martin-Wilkinson contribution that it references. The ALGOL procedures, *imtql1* and *imtql2*, are new.

The authors of Contribution II/4 are A. Dubrulle, R.S.Martin, and J.H.Wilkinson, although the three of them never worked together. Proper credit is given, but I'm afraid that an interesting little bit of history has been lost.

Dubrulle's version of the implicit tridiagonal QR algorithm continues to be important today. We had it in EISPACK. Look at the source for IMTQL1.F. You will see code that is very similar to Austin's note. Except there is a call to a function `PYTHAG` to compute `E(I+1)`.

In LAPACK the *imtql1* and *imtql2* functions are combined into one subroutine named DSTEQR.F Here again you will see code that is very similar to Austin's note. Now `PYTHAG` is named `DLATRG`.

I haven't run any timing experiments myself recently, but I have been told that the method of choice for the real symmetric matrix eigenvalue problem is Cuppen's divide and conquer algorithm, as perfected and implemented in LAPACK's DSTEDC.F. A large problem is recursively broken into smaller ones. When the matrices become small enough, and eigenvectors are desired, DSTEQR is called. This is Dubrulle's version of implicit QR. Then the recursion is unwound to produce the final result.

I could use an appropriate graphic for this post. Here is a picture of the tridiagonal QR algorithm in action, generated by `eigsvdgui` from *Numerical Computing with MATLAB*.

Get
the MATLAB code

Published with MATLAB® R2015a

Jiro's pick this week is Short Path Name on Windows by Jerome Briot.

Thanks, Jerome! Your File Exchange entry helped me big time in my project. I was creating a program that managed a bunch of files, and some of the tasks required me to copy files from one network location to a new location. Sometimes, the new location would have a path name of 261 or more characters. In Windows API, 260 characters is the maximum number of characters allowed in a file path. I can overcome this issue if I map a drive to some path. For example,

\\networklocation\folder1\folder2\folder3\thisfile.txt --> N:\folder3\thisfile.txt

However, I needed to make this a bit more portable so that my team members can use it on their machines. Then I remembered that back in the day, DOS used to abbreviate long path/file names (8.3 filename) to deal with its character limit. After Windows® 95, I've since slowly forgotten about the existence of 8.3 filename. But now, I have a real need for it.

I looked around the Web for the rules for converting a path name to a short path name. That's when I bumped into Jerome's entry. Rather than manually constructing the short file/path names, Jerome's functions uses Windows FileSystemObject (COM) to get the short names. This is obviously much more robust and the correct way.

thisFile = which(mfilename) thisFile_short = fsoGetShortPath(thisFile)

thisFile = C:\Users\jdoke\OneDrive for Business\MyStuff\Work\POTW\2015_06_26_ShortPathName\potw_shortname.m thisFile_short = C:\Users\jdoke\ONEDRI~1\MyStuff\Work\POTW\206E2D~1\POTW_S~1.M

Thanks again, Jerome. It worked perfectly!

**Comments**

Give this a try, let us know what you think here or leave a comment for Jerome.

Get
the MATLAB code

Published with MATLAB® R2015a

Today I'd like to introduce guest blogger Brett Shoelson. Some of you may know Brett through his File Exchange submissions, or through his involvement with the Pick of the Week blog, or from occasional guest posts on Steve’s blog on image processing.

Loren recently told me she had a pending international trip that requires a visa. She had a photograph taken against a background that was unsuitable for use in her visa document, and she asked: "How would I go about changing the background to white?"

"I wonder if this will entail creating a manual mask?", Loren asked.

Great question! I can certainly come up with an automated approach to masking (or "segmenting") the region of interest ("ROI"). For instance:

img = imread('LorenVisaImage.png'); gray = rgb2gray(img); SE = strel('Disk',1,4); morphologicalGradient = imsubtract(imdilate(gray, SE),imerode(gray, SE)); mask = im2bw(morphologicalGradient,0.03); SE = strel('Disk',3,4); mask = imclose(mask, SE); mask = imfill(mask,'holes'); mask = bwareafilt(mask,1); notMask = ~mask; mask = mask | bwpropfilt(notMask,'Area',[-Inf, 5000 - eps(5000)]); showMaskAsOverlay(0.5,mask,'r');

```
%(showMaskAsOverlay is my helper function, available on the MATLAB Central File Exchange.)
```

However, that was complicated by the fact that border portions of Loren's jacket are almost exactly the same color as the background. If she had been standing against a green or blue screen instead, it would have been much easier to automate some chroma key compositing to manipulate the background. Instead, I had to jump through some hoops to get there.

Whenever you are faced with a task like this, it's worth thinking for a moment about the potential return on the investment that automating the task will require. If you need to process many similar images--if this were a frame of a video, for instance--you might indeed want to spend the effort required to automate the segmentation process. On the other hand, if the task is a one-off (as in this case), it might indeed be easier, faster, and potentially more accurate to do things manually. For that, the createMask method of our `imroi` tools facilitates the process.

Creating the code snippet above took a bit of effort. Recreating it with `impoly` took a few minutes:

```
h = impoly(imgca,'closed',false);
```

I much prefer working with `impoly` to working with `imfreehand`, because the former is much easier to adjust after the fact. With `imfreehand`, I find myself tracing borders many times, trying to get it right. But with `impoly`, I can click along, and then adjust individual vertices. I can even delete or add vertices (by pressing the "A" key) as needed to follow difficult contours.

I also prefer working with *open* `imroi` regions for a couple reasons. First, the masks created by `imroi.createMask` are automatically closed by default anyway:

Hmmm. That doesn't *quite* get me the mask I want need. I'm going to modify the `impoly` object by constraining it to the image region:

fcn = makeConstrainToRectFcn('impoly',get(imgca,'XLim'),get(imgca,'YLim')); setPositionConstraintFcn(h,fcn);

Then I'm going to add a couple of vertices, and bring them down to the bottom corners:

Now the mask looks more appropriate for my task:

(I added the red border just to clarify how the mask changed.)

Now, of course, setting the background to white is relatively easy in MATLAB. For a *grayscale* image, we can just use the mask directly:

```
gray = rgb2gray(img);
gray(~mask) = 255;
imshow(gray)
%(The image is of class uint8; 255 is the value of "white" for uint8 images.)
```

First, Loren's image *isn't* grayscale. We'll need to apply that masking function planewise. Second, even for carefully drawn `impoly` regions, the interface between the region we want to keep (i.e., Loren) and the background is pretty blocky, and not very satisfying:

(Okay, maybe it's good enough for a 2"x2" visa picture. But for our purposes here, let's say we're not satisfied with that.)

For the first problem, we can readily break the image into its planewise components, apply the mask, and then reconstruct the color image:

% Break down and mask the planes: r = img(:,:,1); g = img(:,:,2); b = img(:,:,3); r(~mask) = 255; g(~mask) = 255; b(~mask) = 255; % Reconstruct the RGB image: img = cat(3,r,g,b); imshow(img)

To address the second problem (the blocky interface), I'm going to recall a recent discussion I had on the Pick of the Week blog. Capturing the code I discussed in that blog post as a function, and using the same `impoly` I created above, I can create a custom "shell mask" of the Loren-background interface:

```
shellMask = createShellMask(h,'thickness',3);
imshow(shellMask);
```

So why would I do that? Because now I can use the excellent `regionfill` function of the Image Processing Toolbox to interpolate inwards using that mask:

r = regionfill(r,shellMask); g = regionfill(g,shellMask); b = regionfill(b,shellMask); img = cat(3,r,g,b); imshow(img)

Et voila! We have successfully "softened" the interface and created a more natural-looking photograph:

I am interested in hearing whether `createShellMask` is useful for others, and whether I should share it on the File Exchange. Also, the approach I took to the planewise analyses above (masking, and `regionfill`) is fairly generic, and easily converted into a function:

imgout = planewise(fcnhandle,rgbimg,varargin)

If that looks useful to you, let me know and I'll post that as well! Let me know here.

Happy MATLABbing!

Get
the MATLAB code

Published with MATLAB® R2015a

Posted by **Richard Ruff** , June 19, 2015

*Richard is an Application Engineer at MathWorks focused on the Embedded Coder product for code generation, primarily in the Aerospace industry.*

Richard's pick this week is Bus Object Bus Creator by Landon Wagner.

**Pick**

My pick this week is the Simulink library submission for automatically populating the input signal names to a bus creator for a specified bus object: Bus Object Bus Creator.

Do you use bus objects in your Simulink models? Do you spend a lot of time configuring the Bus Creator to match the desired Bus Object? If so, this File Exchange entry is for you.

The description from the entry sums up the benefits nicely:

*When working with bus object types in a design in order to employ strongly "typed" IO the collection of signals into a busCreator associated with a defined bus object type is tedious. After choosing the "Bus:" type in the busCreator "Output data type" the native busCreator does no further favors for you - the user must know or look up how many signals are in the bus type they wish to employ and set "Number of inputs" accordingly. Worse yet, the user must name the individual signals going into the busCreator with the signal names defined in the bus object type. This tool's job is to do that favor for you. (Note that the tool itself is just a masked busCreator running some callback scripts so the inclusion of this tool does not overtly stick out in a design.) Details can be found in BOBCreadme.txt but essentially the tool's mask provides the user with candidate bus object types (From the workspace.) to choose from and upon selection and application of a bus object type the tool provides a busCreator with the necessary number of inputs and with those inputs populated with NAMED stubbed signal lines. These names match the BusElement names in the chosen bus object type. From these named stubbed lines the underlying busCreator inherits the names for "Signals in the bus." Again, this saves the user with having to know these signal names and fill them out themselves.*

A basic example is shown here. Below is a list of the Bus Objects defined in the MATLAB Workspace.

Inspecting the "GuidanceBus" in the Bus Editor, we can see the elements of the bus object: "Ve", "Xe", etc.

Opening the "BusObjectBusCreator" library, we see it contains a single block.

Adding this block to a model and opening the dialog, we can select from the dropdown list of available bus objects. Once we apply the selection, we can see the block is updated to reflect the selected bus object. The number of input ports is altered to reflect the number of elements in the bus and the signals feeding the inputs are labeled to match the bus elements.

**Comments**

Give it a try and let us know what you think here or leave a comment for Landon.

Get
the MATLAB code

Published with MATLAB® R2015a

Have you been interested in data science competitions, but not sure where to begin? Today's guest blogger, Toshi Takeuchi, would like to give a quick tutorial on how to get started with Kaggle using MATLAB.

- The Titanic Competition on Kaggle
- Data Import and Preview
- Establishing the Baseline
- Back to Examining the Data
- Exploratory Data Analysis and Visualization
- Feature Engineering
- Your Secret Weapon - Classification Learner
- Random Forest and Boosted Trees
- Model Evaluation
- Create a Submission File
- Conclusion - Let's Give It a Try

MATLAB is no stranger to competition - the MATLAB Programming Contest continued for over a decade. When it comes to data science competitions, Kaggle is currently one of the most popular destinations and it offers a number of "Getting Started 101" projects you can try before you take on a real one. One of those is Titanic: Machine Learning from Disaster.

The goal of the competition is to predict the survival outcomes for the ill-fated Titanic passengers. You use the training data to build your predictive model and you submit the predicted survival outcomes for the test data. Your score is determined by the prediction accuracy.

Don't worry if you don't rank well on this one. There are entries with a 1.00000 score in the leaderboard, but they either seriously overfit their models to the test data, or perhaps even cheated, given that the full dataset is available from the other sources. That is not only pointless, but also raises serious questions - what kind of standards of conduct must data scientists meet to produce trustworthy results?

So just think of this as a way to do a practice run on Kaggle before you take on a real challenge.

If you haven't done so, sign up with Kaggle - it's free. Then navigate to the Titanic data page to download the following files:

`train.csv`- the training data`test.csv`- the test data

We begin by importing the data into tables in MATLAB. Let's check the imported data. I am assuming that you have downloaded the CSV files into the current folder.

Train = readtable('train.csv','Format','%f%f%f%q%C%f%f%f%q%f%q%C'); Test = readtable('test.csv','Format','%f%f%q%C%f%f%f%q%f%q%C'); disp(Train(1:5,[2:3 5:8 10:11]))

Survived Pclass Sex Age SibSp Parch Fare Cabin ________ ______ ______ ___ _____ _____ ______ ______ 0 3 male 22 1 0 7.25 '' 1 1 female 38 1 0 71.283 'C85' 1 3 female 26 0 0 7.925 '' 1 1 female 35 1 0 53.1 'C123' 0 3 male 35 0 0 8.05 ''

`Train` contains the column `Survived`, and it is the response variable that denotes the survival outcome of the passengers:

1 - Survived 0 - Didn't survive

When you downloaded the data from Kaggle, you probably noticed that additional files were also available - `gendermodel`, `genderclassmodel`, etc. These are simple predictive models that determined the outcome based on the gender or gender and class. When you tabulate the survival outcome by gender, you see that 74.2% of women survived.

disp(grpstats(Train(:,{'Survived','Sex'}), 'Sex'))

Sex GroupCount mean_Survived ______ __________ _____________ female female 314 0.74204 male male 577 0.18891

If we predict all women to survive and all men not to, then our overall accuracy would be 78.68% because we would be correct for women who actually survived as well as men who didn't. This is the baseline Gender Model. Our predictive model needs to do better than that on the training data. Kaggle's leaderboard shows that the score of this model on the test data is 0.76555.

gendermdl = grpstats(Train(:,{'Survived','Sex'}), {'Survived','Sex'}) all_female = (gendermdl.GroupCount('0_male') + gendermdl.GroupCount('1_female'))... / sum(gendermdl.GroupCount)

gendermdl = Survived Sex GroupCount ________ ______ __________ 0_female 0 female 81 0_male 0 male 468 1_female 1 female 233 1_male 1 male 109 all_female = 0.78676

When we looked at `Train`, you probably noticed that some values were missing in the variable `Cabin`. Let's see if we have other variables with missing data. We also want to check if there are any strange values. For example, it would be strange to see 0 in `Fare`. When we make changes to `Train`, we also have to apply the same changes to `Test`.

Train.Fare(Train.Fare == 0) = NaN; % treat 0 fare as NaN Test.Fare(Test.Fare == 0) = NaN; % treat 0 fare as NaN vars = Train.Properties.VariableNames; % extract column names figure imagesc(ismissing(Train)) ax = gca; ax.XTick = 1:12; ax.XTickLabel = vars; ax.XTickLabelRotation = 90; title('Missing Values')

We have 177 passengers with unknown age. There are several ways to deal with missing values. Sometimes you simply remove them, but let's use the average, 29.6991, for simplicity in this case.

avgAge = nanmean(Train.Age) % get average age Train.Age(isnan(Train.Age)) = avgAge; % replace NaN with the average Test.Age(isnan(Test.Age)) = avgAge; % replace NaN with the average

avgAge = 29.699

We have 15 passengers associated with unknown fares. We know their classes, and it is reasonable to assume that fares varied by passenger class.

fare = grpstats(Train(:,{'Pclass','Fare'}),'Pclass'); % get class average disp(fare) for i = 1:height(fare) % for each |Pclass| % apply the class average to missing values Train.Fare(Train.Pclass == i & isnan(Train.Fare)) = fare.mean_Fare(i); Test.Fare(Test.Pclass == i & isnan(Test.Fare)) = fare.mean_Fare(i); end

Pclass GroupCount mean_Fare ______ __________ _________ 1 1 216 86.149 2 2 184 21.359 3 3 491 13.788

With regards to `Cabin`, you notice that some passengers had multiple cabins and they are all in the first class. We will treat missing values as 0. Some third class cabin numbers are irregular and we need to handle those exceptions.

% tokenize the text string by white space train_cabins = cellfun(@strsplit, Train.Cabin, 'UniformOutput', false); test_cabins = cellfun(@strsplit, Test.Cabin, 'UniformOutput', false); % count the number of tokens Train.nCabins = cellfun(@length, train_cabins); Test.nCabins = cellfun(@length, test_cabins); % deal with exceptions - only the first class people had multiple cabins Train.nCabins(Train.Pclass ~= 1 & Train.nCabins > 1,:) = 1; Test.nCabins(Test.Pclass ~= 1 & Test.nCabins > 1,:) = 1; % if |Cabin| is empty, then |nCabins| should be 0 Train.nCabins(cellfun(@isempty, Train.Cabin)) = 0; Test.nCabins(cellfun(@isempty, Test.Cabin)) = 0;

For two passengers, we don't know their port of embarkation. We will use the most frequent value, `S` (Southampton), from this variable to fill in the missing values. We also want to turn this into a numeric variable for later use.

% get most frequent value freqVal = mode(Train.Embarked); % apply it to missling value Train.Embarked(isundefined(Train.Embarked)) = freqVal; Test.Embarked(isundefined(Test.Embarked)) = freqVal; % convert the data type from categorical to double Train.Embarked = double(Train.Embarked); Test.Embarked = double(Test.Embarked);

Let's also turn `Sex` into a numeric variable for later use.

Train.Sex = double(Train.Sex); Test.Sex = double(Test.Sex);

Let's remove variables that we don't plan to use, because they contain too many missing values or unique values.

Train(:,{'Name','Ticket','Cabin'}) = []; Test(:,{'Name','Ticket','Cabin'}) = [];

At this point, we can begin further exploration of the data by visualizing the distribution of variables. This is a time consuming but very important step. To keep it simple, I will just use one example - `Age`. The histogram shows that you have a higher survival rate for agess under 5, and a very low survival rate for ages above 65.

figure histogram(Train.Age(Train.Survived == 0)) % age histogram of non-survivers hold on histogram(Train.Age(Train.Survived == 1)) % age histogram of survivers hold off legend('Didn''t Survive', 'Survived') title('The Titanic Passenger Age Distribution')

How can you take advantage of this visualization? We can create a new variable called `AgeGroup` using discretize() to group values into separate bins like `child`, `teen`, etc.

% group values into separate bins Train.AgeGroup = double(discretize(Train.Age, [0:10:20 65 80], ... 'categorical',{'child','teen','adult','senior'})); Test.AgeGroup = double(discretize(Test.Age, [0:10:20 65 80], ... 'categorical',{'child','teen','adult','senior'}));

Creating such a new variable by processing existing variables is called **feature engineering** and it is a critical step to perform well with the competition and it is where your creativity really comes in. We had already created a new variable `nCabins` to deal with missing data, but often you do this as part of exploratory data analysis. Let's also look at `Fare`.

figure histogram(Train.Fare(Train.Survived == 0)); % fare histogram of non-survivers hold on histogram(Train.Fare(Train.Survived == 1),0:10:520) % fare histogram of survivers hold off legend('Didn''t Survive', 'Survived') title('The Titanic Passenger Fare Distribution') % group values into separate bins Train.FareRange = double(discretize(Train.Fare, [0:10:30, 100, 520], ... 'categorical',{'<10','10-20','20-30','30-100','>100'})); Test.FareRange = double(discretize(Test.Fare, [0:10:30, 100, 520], ... 'categorical',{'<10','10-20','20-30','30-100','>100'}));

The Classification Learner app is a new GUI-based MATLAB app that was introduced in R2015a in Statistics and Machine Learning Toolbox. This will be your secret weapon to try out different algorithms very quickly. Let's launch it!

classificationLearner

- Click on
`Import Data` - Select
`Train`in**Step 1**in`Set Up Classification`dialog box - In
**Step 2**, change the "Import as" value for`PassengerId`to "Do not import", and`Survived`to "Response". All other variables should be already marked as`Predictor`. - In
**Step 3**, just leave it as is to`Cross Validation`.

At this point, we are ready to apply some machine learning algorithms on the dataset. One of the popular algorithms on Kaggle is an ensemble method called **Random Forest**, and it is available as `Bagged Trees` in the app. Let's try that by selecting it from the `classifier` menu and clicking on the `Train` button.

When finished, you can open the `Confusion Matrix` tab. You see that this model achieved 83.7% overall accuracy, which is better than the Gender Model baseline.

`Boosted Trees` is another family of ensemble methods popular among Kaggle participants. You can easily try various options and compare the results in the app. It seems Random Forest is the clear winner here.

You can save the trained model into the workspace by clicking on `Export Model` in the app. If you save the model as `trainedClassifier`, then you can use it on `Test` as follows.

yfit = predict(trainedClassifier, Test{:,trainedClassifier.PredictorNames})

You can also generate a Random Forest model programmatically using `TreeBagger`. Let's adjust the formatting of the data to satisfy its requirements and split the training data into subsets for holdout cross validation.

Y_train = Train.Survived; % slice response variable X_train = Train(:,3:end); % select predictor variables vars = X_train.Properties.VariableNames; % get variable names X_train = table2array(X_train); % convert to a numeric matrix X_test = table2array(Test(:,2:end)); % convert to a numeric matrix categoricalPredictors = {'Pclass', 'Sex', 'Embarked', 'AgeGroup', 'FareRange'}; rng(1); % for reproducibility c = cvpartition(Y_train,'holdout', 0.30); % 30%-holdout cross validation

Now we can train a Random Forest model and get the out-of-bag sampling accuracy metric, which is similar to the error metric from k-fold cross validation. You can generate random indices from the cvpartition object `c` to partition the dataset for training.

% generate a Random Forest model from the partitioned data RF = TreeBagger(200, X_train(training(c),:), Y_train(training(c)),... 'PredictorNames', vars, 'Method','classification',... 'CategoricalPredictors', categoricalPredictors, 'oobvarimp', 'on'); % compute the out-of-bag accuracy oobAccuracy = 1 - oobError(RF, 'mode', 'ensemble')

oobAccuracy = 0.82212

One of the benefits of Random Forest is its feature importance metric, which represents the change in prediction error with or without the presence of a given variable in the out-of-bag sampling process.

[~,order] = sort(RF.OOBPermutedVarDeltaError); % sort the metrics figure barh(RF.OOBPermutedVarDeltaError(order)) % horizontal bar chart title('Feature Importance Metric') ax = gca; ax.YTickLabel = vars(order); % variable names as labels

As expected `Sex` has the most predictive power, but `nCabins`, an engineered feature we came up with, also made a significant contribution. This is why feature engineering is important to do well in the competition! We also used fairly naive ways to fill missing values; you can also be much more creative there.

To get a sense of how well this model actually performs, we want to check it against the holdout data. The accuracy drops significantly against unseen data, and that's what we expect to see when we submit our prediction to Kaggle.

[Yfit, Yscore] = predict(RF, X_train(test(c),:)); % use holdout data cfm = confusionmat(Y_train(test(c)), str2double(Yfit)); % confusion matrix cvAccuracy = sum(cfm(logical(eye(2))))/length(Yfit) % compute accuracy

cvAccuracy = 0.79401

When you tweak your features and modify your parameters, it is useful to use a `perfcurve` plot (performance curve or receiver operating characteristic plot) to compare the performance. Here is an example.

posClass = strcmp(RF.ClassNames,'1'); % get the index of the positive class curves = zeros(2,1); labels = cell(2,1);% pre-allocated variables [rocX, rocY, ~, auc] = perfcurve(Y_train(test(c)),Yscore(:,posClass),'1'); figure curves(1) = plot(rocX, rocY); % use the perfcurve output to plot labels{1} = sprintf('Random Forest - AUC: %.1f%%', auc*100); curves(end) = refline(1,0); set(curves(end),'Color','r'); labels{end} = 'Reference Line - A random classifier'; xlabel('False Positive Rate') ylabel('True Positive Rate') title('ROC Plot') legend(curves, labels, 'Location', 'SouthEast')

To enter your submission to the Kaggle competition, all you have to do is to upload a CSV file. You just need the `PassengerId` and `Survived` columns for submission, and you populate the `Survived` with 1s and 0s. We are going to use the Random Forest model we built to populate this variable.

PassengerId = Test.PassengerId; % extract Passenger Ids Survived = predict(RF, X_test); % generate response variable Survived = str2double(Survived); % convert to double submission = table(PassengerId,Survived); % combine them into a table disp(submission(1:5,:)) % preview the table writetable(submission,'submission.csv') % write to a CSV file

PassengerId Survived ___________ ________ 892 0 893 0 894 0 895 0 896 0

When you upload the submission CSV file, you should see your score immediately, and that would be around the 0.7940 range, putting you within the top 800. I'm pretty sure you are seeing a lot of room for improvement. For example, I just used averages for filling missing values in `Fare` but perhaps you can do better than that given the importance of the feature. Maybe you can come up with better engineered features from the variables I glossed over.

If you want to learn more about how you canget started with Kaggle using MATLAB, please visit our Kaggle page and check out more tutorials and resources. Good luck, and let us know your results here!

Get
the MATLAB code

Published with MATLAB® R2015a

The 1960 programming language ALGOL 60 influenced the design of many subsequent languages, including C and a miniature language called PL/0. I based the design of the first MATLAB on PL/0.

*Algol* is a bright, eclipsing three-star system in the constellation *Perseus*. It provides the name of a family of "algorithmic" programming languages developed very early in the computer era by an international group of, primarily, mathematicians.

I do not have room here to describe ALGOL. Let me just say that ALGOL has profoundly affected my professional life.

In the first published report in 1958, the programming language was originally called IAL, for International Algebraic Language. There were eight authors of the report, including John Backus, an American who headed the IBM team that developed the original Fortran; Friedrich Bauer, a German who was the subject of my previous blog post; and Heinz Rutishauser, a Swiss numerical analyst who invented the LR algorithm.

We had two of versions of ALGOL 58 at Stanford when I was a grad student there in the early '60s. Burroughs had their BALGOL on the B220. We used that for the beginning programming course for about a year. Then two students, Larry Breed and Roger Moore, implemented BALGOL on the new campus IBM mainframe, the 7094. We called it SUBALGOL and used it heavily, including in the programming courses, for a couple of years.

One of my all-time favorite acronyms is JOVIAL, for Jules Own Version of International Algebraic Language. It was developed beginning in 1959 by a group headed by Jules Schwartz at SDC, System Development Corporation. JOVIAL and it successors were in use for almost fifty years by US military projects.

ALGOL, short for Algorithmic Language, was presented in a 1960 report with thirteen authors and edited by Danish computer scientist Peter Naur. The report summarized almost everything that was known about the design of programming languages at the time, which was barely ten years after the development of the first stored program computers. Wikipedia has this to say about ALGOL's influence.

It was the standard method for algorithm description used by the ACM in textbooks and academic sources for more than thirty years. In the sense that most modern languages are "algol-like", it was arguably the most successful of the four high-level programming languages with which it was roughly contemporary: Fortran, Lisp, and COBOL.

The ACM began publishing collections of scientific algorithms in 1960. The collected algorithms appeared in *Communications of the ACM* (*CACM*) until 1975 when *Transactions on Mathematic Software* (*TOMS*) began publication and took over the collection. Initially, all the contributions were required to be in ALGOL 60. Fortran was permitted after it was standardized in 1968. The last ALGOL submission was in 1977. In more recent years, even some MATLAB programs have been accepted.

In the early years of the ACM algorithms section, Stanford played an important role in the editorial process. George Forsythe, who was the founding chairman of Stanford's Computer Science Department and my Ph. D. advisor, was also President of ACM at one time and editor of the algorithms section at another.

My grad school buddy, Bill McKeeman, contributed a couple of early algorithms, in ALGOL. One is algorithm 135, a linear equation solver, "Crout with equilibration and iteration". The Crout algorithm computes the LU factorization from the inner product description of the matrix elements. Another is algorithm 145, "Adaptive numerical integration by Simpson's rule". This one uses recursion even though our BALGOL compiler at Stanford at the time could not handle recursion and Bill could not actually run his program. When we did acquire a Burroughs B5000 with a full ALGOL compiler, Bill verified that his recursive program would work within hours of the machine's installation. (Many years later, one of McKeeman's crowning achievements is the current version of the MATLAB `why` function, which also uses recursion.)

My first ever professional publication is "Remark on Certification of Matrix Inversion Procedures" in the Algorithms section of the July 1963 *CACM*. It is in response to a contribution by Peter Naur, the editor of the ALGOL 60 report. Naur, who was not a numerical analyst, had made the common mistake of comparing the computed inverse of the floating point approximation to the Hilbert matrix to the known inverse of the exact Hilbert matrix. The difficulty is that the effect of the initial floating point approximation is larger than the effect of the rounding errors in the inversion.

For an example of an ALGOL procedure, let's take a look at SOLVE, a program from *Computer Solution of Linear Algebraic Systems*, the 1967 textbook by Forsythe and myself. We have programs in this book in ALGOL, Fortran, and PL/1. This program is preceded by DECOMPOSE, a longer program that computes the *LU* decomposition. Reading these programs today, I'm surprised to see that I used a global array to store the pivot indices. These programs are in the publication format which uses bold face characters for the key words.

In the late 1960s the journal *Numerische Mathematik* ran a series of research papers about numerical linear algebra that includes ALGOL procedures. Altogether there are 29 papers by a combined total of 19 authors. Sixteen of the papers are by J. H. Wilkinson and coauthors. In 1970 all of the papers, with a few modifications, were collected in a volume edited by Wilkinson and C. Reinsch titled "Handbook for Automatic Computation, Volume 2: Linear Algebra". (I'm not sure if Volume 1 of this intended larger project was ever written or published.)

ALGOL was always intended to be a language for research and publication, not for actual execution. IBM dominated the computer market in the '60s and '70s, certainly in the US, and IBM did not support ALGOL. Moreover, in its standard form, ALGOL 60 had no provision for input and output. So it was always expected that the algorithms published by the ACM and *Numerische Mathematik* would have to be modified or translated to some other language before they could actually be used.

The NATS Project, the National Activity to Test Software, was a project at Argonne National Laboratory whose focus was the testing, certification, distribution, and support of mathematical software. Most of the investigators were Argonne staff. I visited Argonne in the summer. The Algol procedures in the *Numerische Mathematik* linear algebra series were organized in two parts. The second part dealt with matrix eigenvalue problems. We translated the procedures in part two to Fortran, carefully preserving their structure and numerical properties. We then tested and timed them at about twenty university and national laboratory sites, and made the collection publically available. "Certification" consisted of a carefully worded offer to respond to any reports of poor or incorrect performance.

The result was EISPACK, for Matrix Eigensystem Routines. We published the Users' Guide in three parts in the 1970s. *Numerische Mathematik* and the Wilkinson/Reinsch *Handbook* had been published by the German publisher Springer-Verlag and we regarded EISPACK as a translation project, so published our guides with them as well.

EISPACK was followed at Argonne by another project, LINPACK, to produce a package for solving linear equations and related problems. Together, the two packages would provide a fairly complete set of routines for computations involving dense matrices. LINPACK was not a NATS project and did not have immediate ALGOL antecedents.

Niklaus Wirth is a Swiss computer scientist who received his PhD from Berkeley in 1963, was on the Stanford faculty at the time the Computer Science department was created in 1965, and returned to the ETH in Switzerland in 1968.

While he was at Stanford, Wirth and Tony Hoare proposed some modifications and simplifications of ALGOL 60 that made it easier to compile and practical for actual execution. The IFIP working group on algorithmic languages felt that the proposal was not a sufficient advance over ALGOL 60, but Wirth implemented it nevertheless. The result became known as ALGOL W, "W" for Wirth, and was used in many programming courses for a number of years.

Wirth is probably best known for his language PASCAL, another descendant of ALGOL. Introduced in 1970, PASCAL is one of the first languages to provide structured programming and data structures. PASCAL proved to be very popular on personal computers in the 1980s.

In 1976 Wirth published a textbook, *Algorithms + Data Structures = Programs*. The last chapter is "Language Structures and Compilers". Wirth gives a formal definition of a miniature programming language that he calls PL/0. He describes how to parse PL/0 and generate code for a hypothetical machine.

PL/0 has only one data type, integer. Wirth defines PL/0 with the seven syntax diagrams for program, block, statement, condition, expression, term, and factor. Here are three of those diagrams.

These definitions are recursive. An expression consists of terms separated by additive operators. A term consists of factors separated by multiplicative operators. And a factor consists of identifiers, numbers, or expressions in parentheses.

The syntax diagram for statement introduces **begin**, **end**, **if**, **then**, **while**, **do**, and **call**. The diagram for block introduces **procedure**. And the diagram for condition involves relational operators. So, despite the simplicity of the definitions, there is a rich structure.

At the time, I was at professor at the University of New Mexico, teaching courses in linear algebra and numerical analysis. I wanted my students to have access to LINPACK and EISPACK without writing Fortran programs. Almost as a hobby, I used Wirth's chapter to develop a simple Matrix Laboratory. Consequently, the syntax of the original MATLAB is that of PL/0. I will discuss this in more detail in my next blog post.

Get
the MATLAB code

Published with MATLAB® R2015a

For more information about programmatically accessing the MATLAB Editor type:

`>>help matlab.desktop.editor`

Unfortunately, not all models can run in Rapid Accelerator mode. In order to generate the stand-alone executable required by the Rapid Accelerator mode, the model must respect certain criteria.

When your model is not compatible with Rapid Accelerator, it is always worth checking if Accelerator mode can speedup your simulation time. However Accelerator Mode does not have an option equivalent to `RapidAcceleratorUpToDateCheck = 'off'`, and consequently every time the simulation starts, Simulink verifies if the model changed.

Today I want to share a trick to skip the initialization in accelerator mode.

**The Problem**

To see the effect of initialization time properly, we need a large model. So I put together a model made of almost 40,000 blocks:

This model simulates almost ten times faster in Accelerator mode compared to Normal mode. The accelerator target takes almost 30 minutes to generate code, but once the code has been generated, it takes about 1 minute to complete one simulation.

However, if we simulate for zero second, we see that just the time to initialize and terminate the model takes about 17 seconds. In other words, a quarter of the time needed to run this model is spent during the initialization.

Let see how we can improve that.

**The Solution: Accelerated Model Reference**

As you are probably aware, model reference allows you to run a model in accelerator mode within another model. I created a model with only one Model block and referenced my large model.

When doing the same test as before, we get:

Oups... the initialization time increased to 32 seconds. In addition to compiling the model, Simulink needs to do some additional work related to model referencing. But don't worry, the story is not over yet!

By default, new models have the Rebuild option set to `If any changes detected`. In this mode, the structural checksum of the referenced model is always computed to ensure that the model has not been modified since the accelerator mex-file has been generated.

Let's open the Model Referencing section of the Model Configuration and set the Rebuild option to `If any changes in known dependencies detected`.

And run our test to measure the initialization time:

In this case, Simulink only verifies that the model.slx file and its dependencies have not changed since the MEX file was generated. If this is the case, the initialization of the referenced model is skipped entirely.

**Now it's your turn**

Are you taking advantage of accelerated model reference to speed up the initialization of your models? Let us know by leaving a comment here.

]]>A few years ago, Jiro wrote a popular post for making pretty plots on this blog. We also host a blog specifically about graphics by Mike. And with the R2014b release of MATLAB came an updated graphics system that Dave described last year in a 3 part series: Part 1, Part 2, and Part 3.

Even with that, I continue to hear questions about how to accomplish certain tasks, such as using a symbol to indicate degrees. This post contains a collection of a few tips that may help you update your plots to match more closely what you are trying to convey.

Suppose we have some temperature data to plot and want to add tick marks that include the degree symbol. Here's some ways to achieve this.

Mock up data for the first of each month this year.

t1 = datetime(2014,1:12,1); temp = [0 2 12 11 15 25 23 27 25 24 12 8];

h = plot(t1,temp,':*'); ax = h.Parent; title('A Year of Temperatures on the 1st of the Month') ylabel('Degrees Celcius ^{\circ}')

ytl = ax.YTickLabel

ytl = '0' '5' '10' '15' '20' '25' '30'

```
ytld = strcat(ytl,'^{\circ}')
```

ytld = '0^{\circ}' '5^{\circ}' '10^{\circ}' '15^{\circ}' '20^{\circ}' '25^{\circ}' '30^{\circ}'

ax.YTickLabel = ytld;

First replace the labels with the original values. Then update the Y axis label. If you want to use the € symbol and it is not on your keyboard, you can use `char(8364)`; $ is char(36).

ax.YTickLabel = ytl; ylabel('Cost in €') title('Income by Month')

Since we want to set all the labels at once, we want to return the results back into a cell array, hence the `false` setting for `'UniformOutput'`.

ax.YTickLabel = cellfun(@(x)sprintf('€ %s',x),ax.YTickLabel,'UniformOutput',false)

ax = Axes (Income by Month) with properties: XLim: [7.3559e+05 7.3594e+05] YLim: [0 30] XScale: 'linear' YScale: 'linear' GridLineStyle: '-' Position: [0.13 0.11 0.775 0.815] Units: 'normalized' Use GET to show all properties

Just create the latex expression for the math, and call `text`, with the `'interpreter'` parameter set to `'latex'`.

plot(magic(3)) t1 = text(1.25,5,'$$\frac{a}{b}$$'); t1.FontSize = 18; t1.Interpreter = 'latex'; t2 = text(2.2,8,'$$\sqrt{3}$$','FontSize',32,'Interpreter','latex');

Do you annotate your plots frequently? What do you need to show? Let us know here.

Get
the MATLAB code

Published with MATLAB® R2015a

Fritz Bauer, eminent German computer scientist and last surviving member of the organizing committee of the 1964 Gatlingburg Conference on Numerical Algebra, passed away on March 26 at the age of 90.

Execute these commands in any MATLAB since version 4.0, 1992.

load gatlin image(X) colormap(map) axis off

This is the organizing committee of the 1964 Conference on Numerical Algebra in Gatlinburg, Tennessee. They are J. H. Wilkinson, Wallace Givens, George Forsythe, Alston Householder, Peter Henrici, and Friedrich Bauer. Bauer, on the far right, was the youngest member of the committee.

I obtained the 8-by-10, black-and-white, glossy photo when I attended the conference as a grad student. It was my first big professional meeting and the first of what has become an important series of meetings for me and MATLAB. I kept the photo in my files until 1992 when we could handle images in MATLAB. I scanned in the photo and it became one of our first example images.

All six of the men in the photo have made contributions that have had a lasting impact on numerical analysis, computer science, and, ultimately, MATLAB. And all six have influenced my life personally.

Credit: <http://upload.wikimedia.org/wikipedia/commons/c/c7/FriedrichLudwigBauer.jpg>

Friedrich Bauer, 1924-2015.

My friend Walter Gander has written an obituary of Bauer for the Communications of the ACM that I reference below. Bauer was born on June 10, 1924, in Regensburg, Germany. After World War II, he studied mathematics and physics at Ludwig Maximillians Universitat in Munich, where he received his Ph.D. in 1951.

Bauer was involved in the early development of computers in Germany. These included STANISLAUS, a relay based computer, in 1951, and PERM, a stored program electronic computer, in 1952-56.

Bauer held a professorship in mathematics at Mainz from 1958 until 1963. He then returned to Munich and the Munich Institute of Technology, where he spent the remainder of his professional career. He advised 39 PhD students before he retired in 1989.

Bauer, together with his colleague Klaus Samelson, invented the stack for use in parsing and evaluating algebraic expressions. Today this is also known as a LIFO, for Last In First Out, data structure.

Bauer was an early advocate for the recognition of *computer science* as an independent discipline in German universities. He also advocated the notion of *software engineering* and, in 1972, suggested a definition.

Establishment and use of sound engineering principles to obtain, economically, software that is reliable and works on real machines efficiently.

This definition of software engineering is now universally quoted.

Bauer was one of the principal authors of the reports on the programming languages International Algebraic Language, IAL, also known as Algol 58, and Algorithmic Language 1960, Algol 60. Wilkinson's research on algorithms for matrix eigenvalues was published in *Numerische Mathematik* in Algol. Equally important, Algol led directly to Niklaus Wirth's pedagogical programming language PL/0, which led to the design of MATLAB. So Bauer had a strong influence on the design of the MATLAB language. I want to tell that story in my next blog post.

The numerical linear algebra community represented by our Gatlinburg photo was very closely knit. Bauer visited Oak Ridge several times. He visited Stanford for a quarter in 1967, stayed in Gene Golub's home, and gave a course where he lectured on the theory of norms. Bauer and Householder were close friends. In fact, Householder married Bauer's wife's older sister, Heidi.

One of my favorite results from the numerical analysis that underlies matrix computation is the *Bauer-Fike Theorem*. It tells how close a computed eigenvalue is to an exact one. You need to be able to estimate the condition of the eigenvectors.

**Bauer-Fike Theorem**. Let $A$ be a diagonalizable matrix and let $V$ be its matrix of eigenvectors. Let $\mu$ and $x$ be an *approximate* eigenvalue and eigenvector with corresponding *residual*

$$ r = A x - \mu x $$

Then there exists $\lambda$, an exact eigenvalue of $A$, such that

$$ | \lambda - \mu | \le \kappa(V) \frac{||r||}{||x||} $$

where

$$ \kappa(V) = ||V|| ||V^{-1}|| $$

is the *condition number* of the eigenvector matrix.

There is a proof in Wikipedia.

I didn't know Fritz Bauer well. I only saw him for a few days every three years at the Gatlinburg/Householder meetings. But the memories are vivid.

In 1996, Householder XIII was in Pontresina, Switzerland. The traditional Wednesday afternoon "break" consisted of an excursion to the Morteratsch Glacier. The adventurous among the world's leading numerical analysts took off on a hike down the glacier, back to the base. I did not want to be left out, although the hiking boots I had hauled to Europe were not in good shape. Halfway down the glacier, a few of us were falling behind. Here comes Fritz and a couple of others. They had taken an longer, more difficult route to inspect a waterfall. He was over seventy years old at the time, but he looked great. He is an avid hiker and was in terrific shape. He was wearing the traditional Alpine lederhosen and first-rate hiking boots.

That was almost twenty years ago, but that's how I'll remember him. Fritz Bauer -- Computer Science Renaissance Man and Bavarian Alpinist.

Walter Gander, "The Life of a Computer Pioneer", Blog@CACM, April 13, 2015, <http://cacm.acm.org/blogs/blog-cacm/185577-the-life-of-a-computer-pioneer/fulltext>

Wikipedia, "Bauer-Fike Theorem", <http://en.wikipedia.org/wiki/Bauer%E2%80%93Fike_theorem>

Get
the MATLAB code

Published with MATLAB® R2015a

**The Problem**

I often get questions similar to the following from users:

I need to model a Pilot-Operated Check Valve. The block should look like the Pilot-Operated Check Valve included with SimHydraulics:

However the behavior I am looking for is slightly different. The only data I have from the supplier is the following figure. When the pressure on the pilot port, multiplied by the pilot ratio, is smaller that pressure on port B, the valve should behave like curve 4, allowing flow only form A to B. When the pressure on the pilot port becomes larger, the valve should go in piloted mode and let flow in both directions, following the characteristic of the dashed curve

**Overview of the solution**

Since all we have are flow-pressure curves, we are going to use an approach similar to the one used in the Fixed Orifice Empirical.

This approach consists of sensing pressures and imposing a flow rate. For a simple orifice, the implementation looks like:

For our Pilot Valve, we will need to sense different pressures and we will need to use two Lookup Tables, one for each mode.

**Getting the Data**

The first thing we need is to do is get data out of the datasheet. For that, I like to use a File Exchange submission titled Data Thief by Adnan.

This submission is very easy to use. If your datasheet is in PDF format, take a screenshot of it and save it as a image file, like a PNG file. Then you can call the Data Thief function by passing it the name of the image file, and the limit values of the curve. In my example, the max value of the pressure axis is 28 bar, the origin is at 0 bar, and 0 l/min, and the maximum flow is 150 l/min. A figure will popup, where you click on the maximum y, origin, maximum x, and then points you want to extract. When you are done, hit enter and the function will return the x and y values in the specified output variables.

Once I get the data, I need to prepare the data so I can cover the full range. In the piloted case, I need to mirror the curve to allow the flow in both direction. In the non-piloted case, I need to set the flow to zero for all the range below the cracking pressure. The code looks like:

This gives us the data we need for our model.

**Option 1: Using blocks**

Now we can use this data in two ways. If you prefer to connect blocks graphically, this option is for you. Using blocks like Pressure Sensor, a Simscape Lookup Table, and a Flow Source, and a few others from the Physical Signals section of the library, we can come up with the following:

When specifying the `x` and `y` values for the lookup table, be careful with the units. The Lookup Table needs to receive and output values in the standard MKS system, here this means Pascals and meter cube per second, while the data we got from the datasheet is in bar and l/min.

One more things to note, I inserted a very small orifice in parallel with the flow source. This is to ensure that our custom valve never generate exactly zero flow. This would behave badly numerically.

**Option 2: Custom Simscape Component**

A second option that will help managing the units more easily is a Simscape composite component. Using this approach, in the `components` section, we declare which blocks we want to use, and in the `connections` section, we define how they are connected together. In a setup section, we can use the `value` function to specify that the values passed to the lookup table block should be in Pascal and m^3/sec. By doing so, the user can specify values in any units he wants, and we take care of the conversion automatically.

The code looks like the following:

**The Results**

To test the valve, I created a model that exercises the valve in its whole range of validity.

We can see that the flow can go in both directions when the valve is piloted, and that it cracks at 10 bar when non-piloted.

**Now it's your turn**

If you are interested in this topic, I also recommend this MATLAB Central File Exchange submission. It contains lots of ressources to model hydraulic components based on datasheets.

Let us know how you model custom components in SimHydraulics by leaving a comment here.

]]>When I wanted to show you a code snippet in the past, I displayed the code in the text of my blog post.

To do so, I often used the function `type`, or sometimes `dbtype` (if I wanted to refer to line numbers). In this way, I did not require you to download code from elsewhere to see what I was discussing.

Now using R2015a I can take advantage of a new feature of the markup language for publish. Using the `include` markup allows me to include external content. How do I do this? With the delimiters

<include>fileOfInterest.m</include>

Not only that, but the content is included with the proper syntax highlighting for MATLAB code!

Here's some code from an older post of mine of methods for computing Fibonacci numbers.

function f = fibrec(n) %FIBREC Recursive function for n Fibonacci numbers. % Minimize the error checking so we don't bog down in it. I have included it %in comments instead. % if ~isscalar(n) | ~isreal(n) | n<0 | fix(n)~=n % error('ArtBlog:fibrec:MBRealPosInt','N must be a real positive integer') % end if n == 1, f = 1; % First element is 1. return; elseif n == 2 f = [1 1]; % First two elements are 1. else % Call fibrec with previous result and compute next one from it. fm1 = fibrec(n-1); f = [fm1 fm1(end)+fm1(end-1)]; end

Let us know how here.

Get
the MATLAB code

Published with MATLAB® R2015a

**Introduction**

Recently, I read an article in Wired (TM) magazine that posed a question that I am sure many of you have confronted -- what is the most efficient way to get passengers onto an airplane. While this article highlights a theoretical (and deterministic) approach to doing this, I set about trying to model this situation and try some variations of this strategy and others that capture at least a few more elements of randomness in the boarding process.

**Generating a queue of passengers**

I used MATLAB to generate a queue of passengers utilizing many different boarding strategies. These strategies are shown in the figure below. The most common one used by airlines is the one labeled 'backToFrontZone'. A colored circle in each plot is more bright than another when that passenger has a higher probability of being in the passenger queue earlier. I use the 'datasample' function in MATLAB to generate a sample passenger queue with those weighted probabilities. Each decrease in shade of green represents a 10 fold decrease in probability of a darker shade dot (i.e. passenger) jumping ahead of the brighter shade dot in the queue.

close all; clear all % Boarding types dictate how you want to order passengers for boarding boardingTypes = {'frontToBackZone', 'backToFrontZone', 'random', ... 'alternateRowZone', 'windowToAisleZone', 'windowToAisleAlternateZone'}; %boardingTypes = {'windowToAisleAlternateZone'}; numRows = 20; % Number of rows in the plane numSeatsPerRow = 6; % Number of seats per row crossZoneFactor = 10; % Factor that controls how likely people from different % zones mix. Higher numbers means mixing less likely expRndMeanStowTime = 1/4; % Mean time (mins) for stowing bags expRndMeanGetUpSitTime = 1/6; % Mean time (mins) for getting up and sitting down % when passenger needs to let someone into % their row seatOrders = []; idealOrders = []; for k = 1:length(boardingTypes) fprintf('Generating %s\n', boardingTypes{k}); % Use MATLAB to produce the passenger order [seatOrder, idealOrder] = airseatSetupSim(boardingTypes{k}, 'repeatable', ... numRows, numSeatsPerRow, ... crossZoneFactor, expRndMeanStowTime, ... expRndMeanGetUpSitTime); seatOrders = [seatOrders, seatOrder]; %#ok idealOrders = [idealOrders, idealOrder]; %#ok end airseatDrawBoardingSchemes(boardingTypes, seatOrders);

The six strategies include:

**frontToBackZone:**Boarding from front of aircraft to back broken into 4 zones.**backToFrontZone:**Boarding from back of aircraft to front broken into 4 zones.**random:**Boarding happens across the aircraft randomly**alternateRowZone:**Boarding (back to front) is done in four zones and each zone only includes alternate rows.**windowToAisleZone:**Boarding broken into 3 zones as ‘window’, ‘middle’, and ‘aisle’**windowToAisleAlternateZone:**Hybrid of 4 and 5 with 6 zones where alternate rows in window, middle, and aisle are treated as separate zones.

**Run simulations**

To model the aircraft cabin, I used a model built in SimEvents. SimEvents is discrete-event simulation library for the Simulink platform. To begin, we generate passenger entities and assign them attributes like row and seat number, along with a random amount of time it will take to stow their bag and sit down.

Then using a library I can model one row of seats and then repeat that row as many times as I want to form the cabin. In each row, I use the Attribute Function block, where I can write MATLAB code accessing and modifying attributes of the passenger entities:

I connected 20 of those in series and got interesting results I will describe below.

To make the simulation more visually appealing, I created a custom visualization, using new features available in R2015a. In MATLAB, I created a simevents.CustomObserverInterface object, which I connected to my model using simevents.connectObserver. In the object, I can specify which blocks to observe using getBlocksToObserve, and then I can code any type of visualization in the entityAdvance to be updated when an entity advances.

Here is what the custom observer object looks like. Note that I removed the actual plotting code to help keeping the focus on the SimEvent feature.

classdef airseatViz < simevents.CustomObserverInterface % AIRSEATVIZ Helper class that draws results of airplane seating model. % properties (Access=private) mFig; mPatches = cell(1, 20); mNumRows; mNumSeats; end % --------------------------------------------------------------------- methods (Access=public) function this = airseatViz(numRows, numSeats) this.mNumRows = numRows; this.mNumSeats = numSeats; end function fig = getFigure(this) fig = this.mFig; end function blks = getBlocksToObserve(this) %#ok<MANU> blks1 = find_system('airseat', 'FollowLinks', 'on', ... 'BlockType', 'SingleServer'); blks2 = find_system('airseat', 'FollowLinks', 'on', ... 'BlockType', 'EntitySink'); blks = union(blks1,blks2); end function p = getPace(this) %#ok<MANU> p = 5; end function initialize(this, ~) % Some code to initialiaze the visualization end function entityAdvance(this, entity, ~, to) enStruct = sedb.eninfo(entity); currRow = enStruct.Attributes.CurrentSeatLoc; destRow = enStruct.Attributes.Row; % Some code to update the visualization when entity advance end function entityDestroy(this, entity, ~) % Code to update the isualization when the entity is destroyed end end end % classdef

Here is what the animation looks like for the case where boarding begins at the back of the plane:

**Results**

I plotted the results for each scheme, along with the ideal case. I define ideal as a perfectly behaved group that boards exactly in the order prescribed by that scheme with no randomness.

The results do validate that the current process used by airlines can only get worse if they decide to board front to back! I was surprised that some schemes work better than the truly random scheme -- but perhaps I should have expected that given the article I started with. It looks like even with the randomness thrown into the deterministic model presented above -- represented by the scheme labeled 'windowToAisleAlternateZone', it performed the best. Maybe when I am in the airport next time things would have changed for the better?

Here are the results for the six schemes described above, normalized by the back to front scheme, which is the standard used by most companies.

**Now it's your turn**

Download Mani's project here and try coming up with a better boarding scheme.

This is the third in a multi-part series on the MATLAB random number generators. MATLAB has used variants of George Marsaglia's ziggurat algorithm to generate normally distributed random numbers for almost twenty years.

It's important to have a memorable name for an algorithm. Ziggurats are ancient Mesopotamian terraced temple mounds that, mathematically, are two-dimensional step functions. A one-dimensional ziggurat underlies Marsaglia's algorithm. I'm not sure if we would still be using the algorithm today if Marsaglia had called it the "step function algorithm".

The probability density function, or *pdf*, of the normal distribution is the bell-shaped curve

$$ f(x) = c e^{-x^2/2} $$

where $c = 1/(2\pi)^{1/2}$ is a normalizing constant that we can ignore. If we were to generate random points $(x,y)$, uniformly distributed in the plane, and reject any of them that do not fall under this curve, the remaining $x$'s form our desired normal distribution. The ziggurat algorithm covers the area under the pdf by a slightly larger area with $n$ sections. The following figure has $n = 6$; the actual code we use today has $n = 256$. The choice of $n$ affects the speed, but not the accuracy of the code.

z = zigplot(6)

z = 0 0.8288 1.1713 1.4696 1.7819 2.1761

The top $n-1$ sections of the ziggurat are rectangles. The bottom section is a rectangle together with an infinite tail under the graph of $f(x)$. The right-hand edges of the rectangles are at the points $z_k$ shown with circles in the figure. With $f(z_1) = 1$ and $f(z_{n+1}) = 0$, the height of the $k$ th section is $f(z_k) - f(z_{k+1})$. The key idea is to choose the $z_k$'s so that all $n$ sections, including the unbounded one on the bottom, have the same area.

There are other algorithms that approximate the area under the pdf with rectangles. The distinguishing features of Marsaglia's algorithm are the facts that the rectangles are horizontal and have equal areas.

For a specified number, $n$, of sections, it is possible to solve a transcendental equation to find $z_n$, the point where the infinite tail meets the last rectangular section. In our picture with $n = 6$, it turns out that $z_n = 2.18$. In an actual code with $n = 256$, $z_n = 3.6542$. Once $z_n$ is known, it is easy to compute the common area of the sections and the other right-hand endpoints $z_k$. It is also possible to compute $\sigma_k = z_{k-1}/z_k$, which is the fraction of each section that lies underneath the section above it. Let's call these fractional sections the *core* of the ziggurat. The right-hand edge of the core is the dotted line in our figure. The remaining portions of the rectangles, to the right of the dotted lines in the area covered by the graph of $f(x)$, are the *tips*. The computation of the $z_k$ 's and $\sigma_k$ 's is done only once and the values included in the header of the source code.

With this initialization, normally distributed random numbers can be computed very quickly. The key portion of the code computes a single random integer, $j$, between $1$ and $n$ and a single uniformly distributed random number, $u$, between $-1$ and $1$. A check is then made to see if $u$ falls in the core of the $j$ th section. If it does, then we know that $u z_j$ is the $x$-coordinate of a point under the pdf and this value can be returned as one sample from the normal distribution. The code looks something like this.

j = randi(256); u = 2*rand-1; if abs(u) < sigma(j) r = u*z(j); else r = tip_computation end

Most of the $\sigma_j$'s are greater than 0.98, and the test is true over 98.5% of the time. One normal random number can usually be computed from one random integer, one random uniform, an if-test, and a multiplication. The point determined by $j$ and $u$ will fall in a tip region less than 1.5% of the time. This happens if $j = 1$ because the top section has no core, if $j$ is between $2$ and $n-1$ and the random point is in one of the tips, or if $j = n$ and the point is in the infinite tail. In these cases, the more costly tip computation is required.

It is important to realize that, even though the ziggurat step function only approximates the probability density function, the resulting distribution is exactly normal. Decreasing $n$ decreases the amount of storage required for the tables and increases the fraction of time that extra computation is required, but does not affect the accuracy. Even with $n = 6$, we would have to do the more costly corrections over 25% of the time, instead of less than 1.5%, but we would still get an exact normal distribution.

Details of the ziggurat implementation have varied over the years. The original 1984 paper by Marsaglia and Tsang included a fairly elaborate transformation algorithm for handling the tips. We used this algorithm for several releases of MATLAB and I reproduced its behavior in the program `randntx` in *Numerical Computing with MATLAB*. That method and that code are now obsolete.

The 2000 paper by Marsaglia and Tsang available at the `jstatsoft` link given below has a simpler rejection algorithm for use in the tips. We have been using this in more recent releases of MATLAB, including current ones.

Marsaglia and Tsang were anxious to make their normal generator as fast as their uniform generator. But they were a little too anxious. Their original code made one call to a 32-bit uniform generator. They used the high order 7 bits for the vertical index $j$ into the ziggurat and then reused all 32 bits to get the horizontal abscissa $u$. Later investigators, including Jurgen Doornik, found this correlation led to failures of certain statistical tests.

We now make two calls to the 32-bit Mersenne Twister generator (during sequential computation). We take 8 bits to get $j$ and then 52 of the remaining 56 bits to get a double precision $u$.

How does this affect the timing? Allocate a long vector.

clear m = 25e6 x = zeros(m,1);

m = 25000000

Generate a random uniform vector and a random normal vector. Then compare the two execution times.

tic, x = rand(m,1); tu = toc tic, x = randn(m,1); tn = toc ratio = tu/tn

tu = 0.3416 tn = 0.4520 ratio = 0.7558

So, random uniform execution time is about three-quarters of the random normal execution time.

Thanks to Peter Perkins for his help on the entire series about the MATLAB random number generators.

This post on the ziggurat is adapted from section 9.3 of *Numerical Computing with MATLAB*.

George Marsaglia and W. W. Tsang, "A fast, easily implemented method for sampling from decreasing or symmetric unimodal density functions." *SIAM Journal on Scientific and Statistical Computing* 5, 349-359, 1984. <http://epubs.siam.org/doi/pdf/10.1137/0905026>

George Marsaglia and W. W. Tsang, "The ziggurat method for generating random variables." *Journal of Statistical Software* 5, 1-7, 2000. <http://www.jstatsoft.org/v05/i08>

Jurgen A. Doornik, "An Improved Ziggurat Method to Generate Normal Random Samples." PDF, Nutfield College, Oxford, 2005. <http://www.doornik.com/research/ziggurat.pdf>

Cleve Moler, *Numerical Computing with MATLAB*, SIAM 2004, <http://dx.doi.org/10.1137/1.9780898717952> Also available at: <http://www.mathworks.com/moler/random.pdf>

Get
the MATLAB code

Published with MATLAB® R2014a

- The player works on devices that do not have Flash such as iPad, and most other mobile devices
- The player can go full screen, which is useful for me as I will often be capturing a larger area of the screen in my videos than Doug did
- The video auto-plays if you click on the thumbnails on the blog index page

A friend just pointed out to me a really cool article: Turns Out Satellites Work Great for Mapping Earthquakes. It's about mapping earthquakes using satellite data. This sounded intriguing because I know earth scientists use MATLAB after the fact to analyze seismic data, but I was less certain what exactly they might do with satellite data and earthquakes. The article, centered on work by Bill Barnett from U. of Iowa, is very interesting. Rather than wait for data to be processed well after events occur, Bill demonstrates what can be done with the data in a much shorter time period, allowing first responders to expedite their decision making regarding where, when, and how to respond to damaging events.

I contacted Bill and verified that he uses MATLAB to perform calculations such as how much the earth moved, what kind of earthquake mechanism is in play, and how large the active fault dimensions are. He shares his code here. This sort of analysis may encourage data providers to make relevant data available as quickly as possible so scientists can collaborate with first responders as soon as possible.

Do you have an interesting story to tell with your data? I bet you do! Let us know here.

Get
the MATLAB code

Published with MATLAB® R2015a

For more information, see the documentation on automatic variable and function highlighting.

- Variable and Subfunction Highlighting in R2010b (4:22)
- Automatic Function and Variable Highlighting (Documentation)
- Renaming Variables (Documentation)

This is the second of a multi-part series about the MATLAB random number generators. If you ask for `help rng`, you will get lots of information, including the fact that there are three modern generators.

'twister' Mersenne Twister 'combRecursive' Combined Multiple Recursive 'multFibonacci' Multiplicative Lagged Fibonacci

My previous post was about `twister`. Today's post is about the other two, `'combRecursive'` and `'multFibonacci'`, which are designed for parallel computing.

I frequently use the card game Blackjack to demonstrate parallel computing. At the same time I can demonstrate the random number generators. I regard Blackjack as a financial instrument, not unlike the stock of a publicly traded company. Simulating the size of an investment as a function of time is a typical application of the Monte Carlo technique.

Begin by opening a pool of workers.

parpool;

Starting parallel pool (parpool) using the 'local' profile ... connected to 2 workers.

Four players each play twenty thousand hands of Blackjack.

p = 4; n = 20000;

Collect the results in this array.

B = zeros(n,p);

A parallel `for` loop executes a Blackjack program that knows nothing about parallel computing.

parfor k = 1:p B(:,k) = cumsum(blackjacksim(n)); end

Plot the results.

plot(B) title(['Blackjack , $ ', int2str(sum(B(end,:)))]) xlabel('Number of hands') ylabel('Stake, ($)') axis([0 n -2500 2500])

The card dealing calls the random number generator. It is essential that the different parallel workers have different, independent streams of random numbers. The default MATLAB generator `twister` does not offer this feature. Simply starting `twister` with different seeds on different workers does not provide statistically independent streams. So we turn to the other generators. To see which one is in use here, run a small `spmd`, "single program, multiple data" block.

spmd r = rng s = r.State' format short x = rand(1,7) end

Lab 1: r = Type: 'combRecursive' Seed: 0 State: [12x1 uint32] s = Columns 1 through 6 1720035765 2052922678 1637499698 3048064580 1173461082 2391850890 Columns 7 through 12 1862757735 2368998908 1385613640 1660833332 146924518 3104031825 x = 0.8789 0.6969 0.0409 0.4609 0.7528 0.2871 0.5241 Lab 2: r = Type: 'combRecursive' Seed: 0 State: [12x1 uint32] s = Columns 1 through 6 323405913 3817048408 3712601073 1070773748 1552739185 3267875480 Columns 7 through 12 1594297407 2533167732 3377045245 3413340742 2651847732 1248925296 x = 0.1072 0.3194 0.1048 0.6623 0.0878 0.3692 0.8035

We see that we have two workers, that they are both using the `combRecursive` generator, that they have the same `seed`, but different `states`, so they are generating different random numbers.

Also known as `mrg32k3a`. A 32-bit combined multiple recursive generator (CMRG), due to Pierre L'Ecuyer, at the Universite de Montreal, and his colleagues, described in the papers referenced below. This generator is similar to the CMRG implemented in the RngStreams package. It has a period of $2^{127}$, and supports up to $2^{63}$ parallel streams, via sequence splitting, and $2^{51}$ substreams each of length $2^{76}$. Here is a link to the C source code. combmrg2.c

The state of the backbone generator is a 2-by-3 array S that evolves at each step according to the linear recurrence expressed succinctly in MATLAB by

m1 = 2^32 - 209; m2 = 2^32 - 22853; x1 = mod(1403580*S(1,2)-810728*S(1,3),m1); x2 = mod(527612*S(2,1)-1370589*S(2,3),m2); S = [x1 S(1,1) S(1,2)) x2 S(2,1) S(2,2)];

A single precision random real `u` is then produced by

z = mod(x1-x2,m1); if z > 0, u = z/(m1+1); else, u = m1/(m1+1); end

The important feature of this generator is that it is possible to create different initial states for each worker in a parallel pool so that the resulting streams of random numbers are statistically independent.

Also known as `mlfg6331_64`. A 64-bit multiplicative lagged Fibonacci generator (MLFG), developed by Michael Mascagni and Ashok Srinivasan at Florida State University. This generator, which has lags $l=63$ and $k=31$, is similar to the MLFG implemented in the SPRNG package. It has a period of approximately $2^{124}$. It supports up to $2^{61}$ parallel streams, via parameterization, and $2^{51}$ substreams each of length $2^{72}$.

The state of this generator is a length 63 vector of 64-bit integers S. The recurrence relation is

$$ x_n = x_{n-k} \times x_{n-l} (mod 2^{64}) $$

Each random double precision value is created using one 64-bit integer from the generator; the possible values are all multiples of $2^{-53}$ strictly within the interval (0,1).

Again, the important feature of this generator is that it is possible to create different initial states for each worker in a parallel pool so that the resulting streams of random numbers are statistically independent.

Which one should you use? Most of the time, stick with the default and you'll be OK. You will get `'twister'` in serial computations and `'combRecursive'` on the workers in a parallel pool. You can use

rng('shuffle')

at the beginning of a session if you want different sequences of random numbers in different sessions. Otherwise, don't worry about setting the generator or the seed.

If you want to experiment, you can use `rng` to try different generators and different starting seeds on your computation. If you find a problem where it makes a significant difference, please let us know.

Thanks again to Peter Perkins.

Pierre L'Ecuyer, R. Simard, E. J. Chen, and W. D. Kelton. "An Objected-Oriented Random-Number Package with Many Long Streams and Substreams." Operations Research, 50(6):1073-1075. 2002. <http://www.iro.umontreal.ca/~lecuyer/myftp/papers/streams00.pdf>

Pierre L'Ecuyer. "Good Parameters and Implementations for Combined Multiple Recursive Random Number Generators." Operations Research 47(1):159-164. 1999. <http://dx.doi.org/10.1287/opre.47.1.159>

Michael Mascagni and Ashok Srinivasan. "Parameterizing Parallel Multiplicative Lagged-Fibonacci Generators." Parallel Computing, 30: 899-916. 2004. <http://www.cs.fsu.edu/~asriniva/papers/mlfg.ps>

Michael Mascagni and Ashok Srinivasan. "SPRNG: A Scalable Library for Pseudorandom Number Generation." ACM Transactions on Mathematical Software, Vol 26 436-461. 2000. <http://www.cs.fsu.edu/~asriniva/papers/sprngacm.ps>

Get
the MATLAB code

Published with MATLAB® R2014b

**Undo Parameter Changes**

For as long as I can remember, in Simulink it has always been possible to undo graphical editing, like moving blocks or adding/deleting them. In R2015a, undo also affects parameter changes.

As you can see in the following recording, when undoing a parameter change, the affected block is quickly highlighted to indicate you that it's the one that changed

Since I edit and debug Simulink models all day, this feature has already saved me a lot of clicks and time.

**Now it's your turn**

Do you like this enhancement? Let us know by leaving a comment here.

]]>Feel free to make requests for topics in the comments area at any time.

]]>

Do you watch movies on Netflix? Binge-watch TV series? Do you use their movie recommendations? Today's guest blogger, Toshi Takeuchi, shares an interesting blog post he saw about how Netflix uses machine learning for movie recommendations.

Back in 2006 Netflix announced a famed machine learning and data mining competition "Netflix Prize" with a $1 million award, finally claimed in 2009. It was a turning post that led to Kaggle and other data science compeititions we see today.

With all the publicity and media attention it got, was it really worth $1 million for Netflix? What did they do with the winning solutions in the end? I came across a very interesting *insider* blog post "Netflix Recommendations: Beyond the 5 stars" that reveals practical insights about what really matters not just for recommender systems but also generally for any real world commercial machine learning applications.

The goal of the NetFlix Prize was to crowdsource a movie recommendation algorithm that delivers 10%+ improvement in prediction accuracy over the existing system. If you use Netflix, you see movies listed under "movies you may like" or "more movies like so-and-so", etc. These days such recommendations are a huge part of internet retail businesses.

It is probably useful to study a very simple example recommendation system based on a well known algorithm called Collaborative Filtering. Here is a toy dataset of movie ratings from 6 fictitious users (columns) for 6 movies released in 2014 (rows).

movies = {'Big Hero 6','Birdman','Boyhood','Gone Girl','The LEGO Movie','Whiplash'}; users = {'Kevin','Jay','Ross','Spencer','Megan', 'Scott'}; ratings = [1.0, 4.0, 3.0, 2.0, NaN, 1.0; 5.0, 1.0, 1.0, 4.0, 5.0, NaN; NaN, 2.0, 2.0, 5.0, 4.0, 5.0; 5.0, NaN, 3.0, 5.0, 4.0, 4.0; 3.0, 2.0, NaN, 3.0, 3.0, 3.0; 4.0, 3.0, 3.0, NaN, 4.0, 4.0];

The idea behind Collaborative Filtering is that you can use the ratings from users who share similar tastes to predict ratings for unrated items. To get an intuition, let's compare the ratings by pairs of users over movies they both rated. The plot of ratings represents their preference space. The best-fit line should go up to the right if the relationship is positive, and it should go down if not.

figure subplot(2,1,1) scatter(ratings(:,1),ratings(:,2),'filled') lsline xlim([0 6]); ylim([0 6]) title('Movie Preference Space by Two Users') xlabel('Kevin''s ratings'); ylabel('Jay''s ratings') for i = 1:size(ratings,1) text(ratings(i,1)+0.05,ratings(i,2),movies{i}) end subplot(2,1,2) scatter(ratings(:,1),ratings(:,4),'filled') lsline xlim([0 6]); ylim([0 6]) xlabel('Kevin''s ratings'); ylabel('Spencer''s ratings') for i = 1:size(ratings,1) text(ratings(i,1)+0.05,ratings(i,4),movies{i}) end

By looking at the slope of the best-fit lines, you can tell that Kevin and Jay don't share similar tastes because their ratings are negatively correlated. Kevin and Spencer, on the other hand, seem to like similar movies.

One popular measure of similarity in Collaborative Filtering is Pearson's correlation coefficient (cosine similarity is another one). It ranges from 1 to -1 where 1 is positive correlation, 0 is no correlation, and -1 is negative correlation. We compute the pairwise correlation of users using rows with no missing values. What are the similarity scores between Kevin and Jay and Kevin and Spencer?

sims = corr(ratings, 'rows', 'pairwise'); fprintf('Similarity between Kevin and Jay: %.2f\n',sims(1,2)) fprintf('Similarity between Kevin and Spencer: %.2f\n',sims(1,4))

Similarity between Kevin and Jay: -0.83 Similarity between Kevin and Spencer: 0.94

Because Kevin and Jay have very different tastes, their similarity is negative. Kevin and Spencer, on the other hand, share highly similar tastes. Users who share similar tastes are called neighbors and we can predict ratings of unrated items by combining their existing ratings for other items. But we need to find those neighbors first. Let's find the neighbors for Kevin.

sims = sims - eye(length(users)); % set self-correlations to 0 kevin_corrs = sims(1,:); [ngh_corr, ngh_idx] = sort(kevin_corrs,'descend'); ngh_corr

ngh_corr = 0.9661 0.9439 0.8528 0 -0.4402 -0.8315

Kevin has three neighbors who have a high correlation with him. We can use their ratings and correlation scores to predict Kevin's ratings. The weighted average method is a basic approach to make predictios. Because the rating scale can be different among individuals, we need to use mean-centered ratings rather than raw ratings. Kevin hasn't rated 'Boyhood' yet. Would he like it?

kevin_mu = nanmean(ratings(:,1)); % Kevin's average rating ngh_corr(4:end) = []; % drop non-neighbors ngh_idx(4:end) = []; % drop non-neighbors ngh_mu = nanmean(ratings(:,ngh_idx),1); % neighbor average ratings Predicted = nan(length(movies),1); % initialize an accumulator for i = 1:length(movies) % loop over movies ngh_r = ratings(i,ngh_idx); % neighbor ratings for the movie isRated = ~isnan(ngh_r); % only use neighbors who rated meanCentered =... % mean centered weighted average (ngh_r(isRated) - ngh_mu(isRated)) * ngh_corr(isRated)'... / sum(ngh_corr(isRated)); Predicted(i) = kevin_mu + meanCentered; % add Kevin's average end Actual = ratings(:,1); % Kevin's actual ratings table(Actual, Predicted,'RowNames',movies) % compare them to predicted fprintf('Predicted rating for "%s": %.d\n',movies{3},round(Predicted(3)))

ans = Actual Predicted ______ _________ Big Hero 6 1 1.4965 Birdman 5 4.1797 Boyhood NaN 4.5695 Gone Girl 5 4.2198 The LEGO Movie 3 2.8781 Whiplash 4 3.9187 Predicted rating for "Boyhood": 5

Looks like Kevin would rate 'Boyhood' as a 5-star movie. Now, how accurate was our prediction? The metric used in the Netflix Prize was Root Mean Square Error (RMSE). Let's apply it to this case.

RMSE = sqrt(nanmean((Predicted - Actual).^2))

RMSE = 0.5567

Now you saw how basic Collaborative Filtering worked, but an actual commercial system is naturally much more complex. For example, you don't usually see raw predicted ratings on the web pages. Recommendations are typically delivered as a top-N list. If so, is RMSE really a meaningful metric? For a competition, Netflix had to pick a single number metric to determine the winner. But choosing this metric had its consequences.

The goal of the competition was to crowdsource improved algorithms. $1M for 3 years of R&D? One of the teams spent more than 2000 hours of work to deliver 8.43% improvement. Combined, a lot of people spent enormous amounts of time for a slim hope of the prize. Did Netflix use the solutions in the end?

They did adopt one solution with 8.43% improvement, but they had to overcome its limitations first.

- The number of ratings in the competition dataset was 100 million, but the actual production system had over 5 billion
- The competition dataset was static, but the number of ratings in the production system keeps growing (4 million ratings per day when the blog post was written)

They didn't adopt the grand prize solution that achieved 10% improvement.

- Additional accuracy gains that we measured did not justify the engineering effort needed to bring them into a production environment
- The Netflix business model changed from DVD rental to streaming and that changed the way data is collected and recommendations are delivered.

Why additional accuracy gains may not be worth the effort? For example, you could improve the RMSE by closing the prediction gaps in lower ratings - movies people would hate. Does that help end users? Does that increase revenue for Netflix? Probably not. What you can see here is that, in a production system, **scalability** and **adaptability** to changing business needs are bigger challenges than RMSE.

Had Netflix chosen for the competition some metrics more aligned with the needs of the production system, they might have had an easier time adopting the resulting solutions.

You often learn more from what didn't go well than what went well. Netflix was not able to take full advantage of the winning solutions, but they certainly appear to have learned good lessons based on how they now operate now. I would say they got well more than $1M's worth from this bold experiment. Let's see how Netflix is taking advantage of lessons learned.

When Netlix talks about their current system, it is notable what they highlight.

- "75% of what people watch is from some sort of recommendation"
- "continuously optimizing the member experience and have measured significant gains in member satisfaction"

What they now care about is usage, user experience, user satisfaction, and user retention. Naturally they are also well aligned with the bottomline of Netflix as well. The second bullet point refers to A/B testing Netflix is conducting on the live production system. That means they are constantly changing the system...

"Coming up with a software architecture that handles large volumes of existing data, is responsive to user interactions, and makes it easy to experiment with new recommendation approaches is not a trivial task." writes the Netflix blogger in "System Architectures for Personalization and Recommendation". You can also see more detail in this presentation.

One of the techniques used in the winning solutions in the Netflix Prize was an ensemble method called Feature-Weighted Linear Stacking. Netflix apparently adopted a form of linear stacking technique to combine the predictions from multiple predictive models to produce final recommendations. The blog post gives an example: if `u` = user, `v` = video item, `p` = popularity, `r` = predicted rating, `b` = intercept, and `w1` and `w2` are respective weights, then a simple combination of popularity score and predicted ratings by Collaborative Filtering can be expressed as:

$$f_{rank}(u,m) = w1 \cdot p(v) + w2 \cdot r(u,v) + b$$

And you can just add more terms to this linear equation as you develop more predictive models, each running on its subsystem... This is a very flexible architecture and makes it easy to run an A/B test - it is just a matter of changing weights!

Netflix uses three layers of service to achieve this - offline, nearline and online to overcome this challenge.

- Offline to process data - pre-computes the time-consuming steps in a batch process.
- Online to process requests - responds to user action instantaneously by taking advantage of the Offlne and Nearline outputs
- Nearline to process events - bridges the two sub-systems by pre-computing frequently performed operations and caching the result in advance of active user action

As a simple example for an offline process, let's use the MovieLens dataset with 100K ratings from 943 users on 1682 movies.

data_dir = 'ml-100k'; if exist(data_dir,'dir') ~= 7 unzip('http://files.grouplens.org/datasets/movielens/ml-100k.zip') end data = readtable(fullfile(data_dir,'u.data'),'FileType','text',... 'ReadVariableNames',false,'Format','%f%f%f%f'); data.Properties.VariableNames = {'user_id','movie_id','rating','timestamp'}; sp_mur = sparse(data.movie_id,data.user_id,data.rating); [m,n] = size(sp_mur); figure spy(sp_mur) title('Movie-User Matrix Sparsity Pattern') xlabel('users') ylabel('movies')

Older movies in lower row indices have fairly dense ratings but the newer movies with higher row indices are really sparse. You can also see that older users in lower column indices don’t rate newer movies at all – it seems they have stopped using the service.

You can take advantage of sparsity by computing directly on the sparse matrix. For example, you can pre-compute the Pearson correlation scores and mean centered ratings, and, once we have an active user, they can be used to compute the neighborhood and generate recommendations for that user.

It turns out you need to update the pre-computation of neighborhood and mean-centered ratings fairly frequently because most users don't rate many movies and one new rating can shift values a lot. For this reason, the item-based approach is used more commonly than user-based approach we studied earlier. Only big difference is that you compute similarity based on movies rather than users. It is just a matter of transposing the mean-centered matrix, but you need to update this less frequently because item-based scores are more stable.

meanRatings = sum(sp_mur,2)./sum(sp_mur~=0,2); [i,j,v] = find(sp_mur); meanCentered = sparse(i,j,v - meanRatings(i),m,n); sims = corr(meanCentered', 'rows', 'pairwise');

We can still do this in-memory with a 100K dataset, but for larger datasets, MATLAB can run MapReduce jobs with MATLAB Distributed Computing Server on Hadoop clusters for offline jobs. MATLAB Production Server enables rapid deployment of MATLAB code in Production environment and it may be a good choice for nearline or online uses where concurrent A/B testing is performed, because you avoid dual implementation and enable rapid system update based on the test result.

For an example of how you use MapReduce in MATLAB, check out "Scaling Market Basket Analysis with MapReduce". Incidentally, Collaborative Filtering and Association Rule Mining are related concepts. In Market Basket Analysis, we consider transactions as basic units. In Collaborative Filtering, we consider users as basic units. You may be able to repurpose the MapReduce code from that post with appropriate modifications.

To learn more about MATLAB capabilities, please consult the following resources.

As noted earlier, Collaborative Filtering and Market Basket Analysis are closely related applications of data mining and machine learning. Both aim to learn from the customer behaviors, but Market Basket Analysis aims for high frequency transactions while Collaborative Filtering enables personalized recommendations. Naturally, you cannot restock your store shevles for individual shoppers, but you can in ecommerce. You can also think of Latent Semantic Analysis as a related technique applied in text analytics domain.

Just as Market Basket Analysis found its way into web usage data mining, you can probably use Colaborative Filtering for the web data analysis. There may be applications in other fields as well. Any thoughts about your creative use of "Collaborative Filtering" with your data? Let us know here.

Get
the MATLAB code

Published with MATLAB® R2015a

This is the first of a multi-part series about the MATLAB random number generators.

If you issue the following commands at any point in any recent version of MATLAB, you will always get this plot.

```
rng default
hist(randn(10000,1),100)
```

The `rng` command controls the random number generator that is used by the `rand`, `randn`, and `randi` functions. When called with the `default` parameter, `rng` resets the generator to the condition that it has when a fresh MATLAB is started. To see what that condition is, just call `rng` by itself.

rng

ans = Type: 'twister' Seed: 0 State: [625x1 uint32]

We see that the default generator is `'twister'`, the default seed is 0, and that the state is a length 625 vector of unsigned 32-bit integers.

If you ask for `help rng`, you will get lots of information, including the fact that there are three modern generators.

Generator Description ------------------------------------------------------ 'twister' Mersenne Twister 'combRecursive' Combined Multiple Recursive 'multFibonacci' Multiplicative Lagged Fibonacci

The remainder of today's post is about `'twister'`. I will cover the others in future posts.

Mersenne Twister is, by far, today's most popular pseudorandom number generator. It is used by every widely distributed mathematical software package. It has been available as an option in MATLAB since it was invented and has been the default for almost a decade.

Mersenne Twister was developed by professors Makoto Matsumoto and Takuji Nishimura of Hiroshima University almost twenty years ago. Here is their home page. The C source code is available here.

Mersenne primes are primes of the form *2^p - 1* where *p* itself is prime. They are named after a French friar who studied them in the early 17th century. We learn from Wikipedia that the largest known prime number is the Mersenne prime with *p* equal to 57,885,161. The Mersenne Twister has *p* equal to 19937. This is tiny as far as Mersenne primes go, but huge as far as random number generators are concerned.

Matsumoto explains how the name "Mersenne Twister" originated in the Mersenne Twister Home Page.

MT was firstly named "Primitive Twisted Generalized Feedback Shift Register Sequence" by a historical reason. Makoto: Prof. Knuth said in his letter "the name is mouthful." Takuji: ........ a few days later Makoto: Hi, Takkun, How about "Mersenne Twister?" Since it uses Mersenne primes, and it shows that it has its ancestor Twisted GFSR. Takuji: Well. Makoto: It sounds like a jet coaster, so it sounds quite fast, easy to remember and easy to pronounce. Moreover, although it is a secret, it hides in its name the initials of the inventors. Takuji: ....... Makoto: Come on, let's go with MT! Takuji: ....well, affirmative. Later, we got a letter from Prof. Knuth saying "it sounds a nice name." :-)

The integer portion of the Mersenne twister algorithm does not involve any arithmetic in the sense of addition, subtraction, multiplication or division. All the operations are shifts, and's, or's, and xor's.

All the elements of the state, except the last, are unsigned 32 bit random integers that form a cache which is carefully generated upon startup. This generation is triggered by a `seed`, a single integer that initiates the whole process.

The last element of the state is a pointer into the cache. Each request for a random integer causes an element to be withdrawn from the cache and the pointer incremented. The element is "tempered" with additional logical operations to improve the randomness. When the pointer reaches the end of the cache, the cache is refilled with another 623 elements.

The algorithm is analyzed by investigating the group theoretic properties of the permutation and tempering operations. The parameters have been chosen so that the period is the Mersenne prime 2^19937-1. This period is much longer than any other random number generator proposed before or since and is one of the reasons for MT's popularity.

By design, the results generated satisfy an equidistribution property in a 623-dimensional cube.

Here is the function in the Mersenne Twister source code that converts a pair of random uint32s into a random double. You can see that it takes the top 27 bits of one int and the top 26 bits of the other, sticks them together, and multiplies by what MATLAB would call `eps/2`. This is the only place in the code where floating point arithmetic is involved.

double genrand_res53(void) { unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6; return(a*67108864.0+b)*(1.0/9007199254740992.0); }

The result would be zero in the highly unlikely event that both `a` and `b` are zero. If that happens, the MATLAB interface rejects the result and calls this function again. So the smallest double results when `a` equals zero and `b` equals 1. The largest double results when both `a` and `b` are all 1's. Consequently, the output from `rand` is in the closed interval

$$ 2^{-53} \leq x \leq 1-2^{-53} $$

For more about random number seeds, streams, and state, see Peter Perkins, guest blogger in Loren's blog. And, of course, see the documentation.

Thanks to Peter Perkins for the work he has done on our random number suite over the years, and for enlightening me.

M. Matsumoto and T. Nishimura. "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudorandom Number Generator." ACM Transactions on Modeling and Computer Simulation, 8(1):3-30. 1998. Available online at: <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf>.

Get
the MATLAB code

Published with MATLAB® R2014b

There are 42 muscles in the human face. When different combinations of these muscles are flexed to varying degrees, a single face can take on millions of different positions. While even the most skilled experts can’t discern what someone is thinking based on their facial expressions, computers can build models based on the aggregation of billions of data points.

MATLAB already has a host of tools to build machine learning models and process big data. It also has the ability to identify and track human faces and facial features. Bringing all these technologies together, we’ve built a revolutionary tool that will change the way you use MATLAB. It’s called Face Coder, and early trials have shown promising results.

To use Face Coder, all you have to do is set up a webcam and look at your computer screen. Instead of having to type all your code out by hand, a slow and error-prone process, MATLAB will interpret your facial expressions and automatically generate the code you intend to write.

Most of our Beta feedback has been very positive:

- “Face Coder has changed the way I write code. I can stare at my screen for hours thinking about MATLAB and never have to wait for my fingers to catch up to my thoughts.”
- “With Face Coder, I feel a much closer connection to my computer. It’s almost as if MATLAB can read my mind.”
- “I hope Face Coder will come out on MATLAB Mobile. I love to code, but I also love to maintain an active lifestyle. My dream is to someday be able to write code while I’m out cycling or kayaking.”

But other feedback has shown that we still have some work left to do:

- “I often get distracted while I’m working. With Face Coder, I sometimes refocus on my work only to find that MATLAB has typed out my thoughts about bills I need to pay or what I want to eat for lunch.”
- “Much of my work is classified, and it would be dangerous if unauthorized persons gained access to my MATLAB code. Therefore, I’m concerned that this feature will be a security vulnerability. I’ll need to wear sunglasses or a mask while programming to prevent anyone from recording and then interpreting my face.”

Though Face Coder is not yet available, all the technologies it’s built on are available in current MATLAB tools. These include:

]]>Learning by watching videos has become standard technique for students of everything from piano to surgery. And increasingly people are learning to code by watching other people code.

On this site we have lots of clever people solving MATLAB problems on Cody. I would love to see somebody post a video to YouTube that shows them actually in the process of solving a Cody problem. Or maybe just describing after-the-fact how they did solve it. Or even analyzing, sports-commentator style, the various techniques used to solve a particular problem. If you make a “How I Solved It” video for Cody leave a link here and tell us about it!

And just to put my money where my mouth is, I have provided a humble example here: How I Solved It: MATLAB Pangrams (Cody Problem 27).

]]>*by David Wey*

Have you ever wondered who created the file you’re looking at on the File Exchange, or who that nice person was that answered your question on MATLAB Answers? Maybe you thought to yourself *I wonder what else that person has done?* As of this week answering these questions got a lot easier.

**NEW PROFILE PAGES**

Our latest update to MATLAB Central includes new profile pages. These pages aggregate content from across all the MATLAB Central areas, thus easing the pain of having separate profile pages for each area (Answers, Cody, File Exchange, and so on). On these profile pages you’ll be able to see how active the person has been and how recently they’ve participated in the community. You will see all their answers, files, Cody problems, trends, and links.

To view someone’s new profile page, click on any author link wherever you see it, whether in Cody or the File Exchange. The popup displays a summary of information related to that area of MATLAB Central (such as the File Exchange or Answers). Click again on the user name to view the profile page or click on the other labels to see their submissions as a filtered search.

We hope you’ll enjoy the new profile pages and we welcome your thoughts and reactions. Check out some of your favorite community member pages, as well as your own, and let us know what you think in the comments below.

]]>The first game I thought of was Connect Four, since it’s a straightforward extension of Tic Tac Toe. And sure enough, it’s not to hard to convert our Tic Tac Toe engine into a Connect Four engine.

**CONNECT FOUR**

In Connect Four, your board is a 6-by-7 grid, and your goal is to be the first to put four pieces in a line. The game board is oriented vertically, so gravity causes your pieces to drop to the bottom of the grid. This means that you have at most seven legal moves. This is good for us, since it constrains the number of possible moves, which in turn makes it quicker for us to play hundreds of games.

The game harness botMoves doesn’t change at all. We just make a new class, ConnectFour, that knows the rules and can draw the Connect Four game board. And when I say “knows the rules”, I mean exactly two things.

- What legal moves can I make right now?
- If the game is over, who won?

That’s it.

Finding legal moves is easy. We just look for any column that has zeroes left in it. But how do you determine if someone has won the game? We need to scan the board for four-in-a-row of either color in any of four directions: vertical, horizontal, and two diagonals. This seemed like it was going to call for some careful thinking, until I remembered that there was a Cody problem about this very topic! See Problem 90. Connect Four Win Checker.

I shopped around in the answers to Problem 90 and found a nice piece of code authored by long time Cody champ @bmtran: Solution 2314.

I wired everything together, and it plays a pretty mean game. Here’s an animation of my Connect Four bot playing itself.

**CONNECT FOUR VARIANTS**

The Monte Carlo approach means that we haven’t had to put a lot of special game-specific smarts into the code. This gives us great freedom to explore variants on our basic Connect Four game. Suppose, instead of racing to get four in a row, players were trying to be the first to make an “L” shape. In order to make this work, I just had to change a few lines in the isGameOver method. And rather than clone and modify the entire ConnectFour class, I built a FourGameBase class from which I could quickly subclass variants. This base class knows how to draw the board and pick legal moves. The winning conditions are handled by the subclasses.

Here’s my ConnectEl game bot playing itself.

We can imagine a whole family of Connect Tetris variants: be the first to make a square, or a T-shape, or an S-shape. Not all of these variants are interesting as games. I built square and T-shape versions of the game, and in both cases defensive play was so easy that, as with Tic Tac Toe, all the games ended in ties. But this ability to rapidly create and play games leads to another observation. We can use our Monte Carlo code to mutate games in search of interesting variants. An interesting game is one in which either side has a reasonable chance of winning. An uninteresting game consistently ends in a tie or victory for a particular side. So we can effectively use bots to do our play-testing for us. The optimal strategies simply emerge from the Monte Carlo soup.

I want to mention one last Connect Four variant before I move on: Four Corners. In Four Corners, you try to force your opponent to be the first to put his pieces into the four corners of a rectangle.

In this case, I used code from a Cody answer by Alfonso Nieto-Castanon to test for the victory condition: Spot the rectangle.

Here is my bot slugging it out with itself in Four Corners.

This is a fun game to play against the bot. It’s surprisingly tricky to spot all the potential rectangles on the board.

**REVERSI**

For my last example, I wrote a program that does a passable job of playing Reversi. As before, I used Cody to diminish the coding effort by posting two problems.

- Problem 2538. Find the Next Legal Move in Reversi (thanks Tim!)
- Problem 2565. Determine the Result of a Move in Reversi (thanks Jan!)

Here’s the animation of the game in action.

**WHAT NEXT?**

If you, like me, are intrigued with the possibilities of Monte Carlo gaming, I encourage you to take a look at my code on GitHub. Maybe you can add some new games to the list that I’ve already created. There are plenty of games that would work with this harness. Or you can improve the code that’s already there, making it speedier or more efficient. Because of the brute force nature of the Monte Carlo technique, it would be an ideal place to experiment with parallel programming. Or you might make it more interactive, with a click-to-move interface.

I’ll stop here, but this is clearly a rich space to explore. It’s a powerful reminder that computers are very different from brains. Sometimes you can make up for not being clever by being stupid very fast.

]]>**COMPUTERS, CHESS, AND GO**

I read an article in IEEE Spectrum about computer programs that play Go (AIs Have Mastered Chess. Will Go Be Next?). If you review the history of game-playing computers, you’ll see that chess programs improved steadily until eventually they could beat the best human players. Go programs, on the other hand, have been stuck at a level of play that was nowhere close to the best human. Why is that?

The basic element of a game-playing program is the look-ahead. Essentially, the program says “If I move here, is that better or worse than if I move there?” In chess, this is straightforward to evaluate. But in Go, this basic look-ahead strategy doesn’t work so well. It’s much harder to evaluate whether one board position is stronger than another.

But recently, Go programs have started to get much better. What happened?

**TWO IDIOTS FINISH THE GAME**

Go programs have improved by applying a Monte Carlo technique. It’s nothing like how a human plays, but it works remarkably well. And it only works because we can ask the computer to do a lot of dumb stuff very quickly. I call it “Two Idiots Finish the Game”.

Consider the following situation. You’ve reached a critical point in the game. We’ll call it position X. You’re considering

move A and move B. Which one should you make? Now instead of looking just one move ahead, play the game all the way to completion. But there’s an obvious problem with this. If you’re not smart enough to figure out your next move, how can you play an entire game? Simple: just ask two idiots to make random (but legal) moves until one of them wins. Then return the game to position X and have them play again. And again. And again and again. Sometimes they start with move A, and sometimes B. After your speedy but not-so-clever friends have played a few thousand games, examine the record. Is an idiot (with an idiot for an opponent) more likely to win with move A or move B? Those simulated games will give you the answer. Here’s the amazing thing: the idiot’s best move is your best move too. Don’t ask one clever mouse to solve the maze. Release ten thousand stupid mice and follow the lucky ones. This is what cheap computation buys you.

What’s beautiful about this approach is that it’s completely free of strategy. You don’t need to build up special knowledge structures about any particular game. You just need to know what moves are legal and how the game ends.

**TIC TAC TOE**

As soon as I read about this technique, I wanted to try it in MATLAB. So let’s make a program that can play Tic Tac Toe (also known as Naughts and Crosses). I’ve written Tic Tac Toe programs in MATLAB before. I’ve tried to make them clever and I’ve tried to make them learn. It’s not that hard. What’s fun about this Monte Carlo approach is that, with minimal effort I can teach it a new game. In fact, it makes playing lots of games easy. With a little object-oriented programming, you can write a generic game-playing harness. Then you just need to plug in some code that knows a few rules, and presto! You’ve got an instant game-playing program.

Here’s what I did. I made a class called TicTacToe that knows the rules of the game and how to draw the board. Then I wrote a function called botMoves that can look at the game object and make the next move. The separation is very clean. All of the Monte Carlo logic mentioned above lives in botMoves.

I only need a short script to have the bot play itself.

game = TicTacToe; nSimulatedGames = 1000; while ~game.isGameOver botMoves(game,nSimulatedGames); end

The variable nSimulatedGames refers to the number of simulated games we’ll ask our idiot friends to play for each potential move. Here’s an animation of what it looks like in action.

As it happens, the computer always ties itself. That’s actually good news, since Tic Tac Toe is unwinnable if your opponent is the least bit clever. So our bot is smart enough to prevent itself from winning. A little play-testing shows that it’s smart enough to avoid losing to a human too. But if we prefer, we can make the program less competitive by lowering the number of simulated games it plays. If I only let it run ten simulated games for each possible move, I can beat it easily.

I haven’t displayed much of my code here in the blog, but you can get your hands on it at this GitHub repository: Monte-Carlo-Games. Here is the TicTacToe class, and here is the botMoves function.

**NEXT WEEK**

This is the first of a two-part post. Next time we’ll show how quickly we can adapt our simple Tic Tac Toe harness for other games. We’ll also bring a community element into our programming. We’ll use Cody to source some of the tricky parts of our coding effort!

*by Chaitanya Chitale*

Are you new to MATLAB and looking for a way to get started? If so, check out the new MATLAB Onramp course.

MATLAB Onramp is an online course that provides a brief introduction to the MATLAB language. The course gives you hands-on MATLAB experience via the use of an integrated, web-based version of MATLAB, as shown below.

As you work through the course you can try out different solutions to the problems asked and get immediate feedback.

The MATLAB Onramp course covers the basics of importing data, manipulating arrays, creating visualizations, and much more.

The course takes approximately 2 hours to complete. However, you can take the course at your own pace. After starting it, you can leave and come back to it at any time.

The MATLAB Onramp course is complimentary with your purchase of MATLAB. To get started, head over to the MATLAB Academy site: https://matlabacademy.mathworks.com.

You can also access the course from MATLAB R2014b by selecting Help → MATLAB Academy from the Home tab.

Interested in other training? Explore the other training courses that are offered by MathWorks.

]]>*by Wendy Fullam*

Have you ever found yourself thinking, “I wish there was a place where I could go to see a bunch of MATLAB examples…”?

MathWorks has just what you’re looking for. We just launched a big new feature on the Support page: MATLAB Examples. On it, you can discover thousands of code examples for MATLAB and Simulink: data manipulation, signal processing, machine learning, statistics, finance, computational biology, finance, you name it. Maybe you’ll even discover a few things that you didn’t know were possible.

MATLAB Examples is a single destination to find high-quality code examples, including those authored by both MathWorks staff and contributors to the File Exchange. Here’s what the main page looks like.

What can you do here?

You can browse by topic area:

You can search by keyword, and filter the results:

After finding an example that you want to use, you can quickly see what product contains that example:

or access the code:

You can also find related examples:

The Help Team at MathWorks hopes you enjoy the new interface and welcomes your thoughts and reactions. Check out the new area and then let us know what you think in the comments below!

]]>by Toshi Takeuchi

One of the new R2014b features that deserves your attention is Git integration. Git is a source control system (also known as version control or source code management system) that enables collaborative software development. Why does that matter to you? Programming is an essential skill in many technical fields even outside computer science, and some universities now offer software carpentry workshops to enhance coding skills for researchers. Source control is one of those essential skills in software carpentry.

Until now, you may have tinkered alone with code you needed for your project. However, there are other people who may be working on similar problems and they may be writing similar programs. Source control enables you to work with other people so that you don’t have to do it all alone. Collaboration lets you be more productive in other aspects of your project.

Even if you don’t care about such collaboration, wouldn’t it be cool to share your personal project and see other people using it? They may even fix bugs and improve your code for you!

GitHub is one of the most popular websites that host Git repositories. The best place to share your MATLAB projects is File Exchange because of its popularity with the MATLAB user community. And guess what – File Exchange is integrated with GitHub! Now you see the connection?

What is a Git repository? A repo (repository) is a directory that holds your source code and any associated files. Local repos are on your local drive, and the remote repos are on GitHub or other hosts, and you sync the local repos to remote repos as you write your code. You can start with either the local or remote repos, but in this example I am going to start with a remote repo.

The process looks like this for a single developer:

- Create or fork a repo on GitHub
- Clone the repo to your local drive – this is your local repo
- Add your files to the local repo
- Sync your local repo to remote repo
- Repeat this process to keep your source code in sync as you write more code
- Share your GitHub repo on File Exchange

When you talk about Git, you cannot go without mentioning “forking”. In the simplest terms forking means copying someone else’s public repo on the remote server, rather than starting a repo from scratch. In practice forking is used as a way to contribute to the existing projects or to start a new project using an existing project as a starting point. Once you make changes to your forked project, you can send a merge request to the original developer, and your changes may be accepted and incorporated into the main project.

Forking enables a flexible distributed style of collaboration and number of forks you have on your project acts as a measure of popularity – similar to the count of likes or followers on Facebook or Twitter. The social aspect of forking is an interesting topic on its own, but we need to skip it for this post.

Signing up with GitHub is easy – just click the *sign up* button on the homepage and follow instructions. A free account would do for now. You also need to download and install Git. Even though GitHub has GUI apps for Windows and Mac, you need to set up the command-line tool for use with MATLAB. You also want to follow these instructions for registering binary files with Git.

Creating a repo on GitHub is very easy – just follow these instructions. From this point on I will assume you named your repo *Hello-World* and initialized it with a README file. Please note that you can only create a public repo with a free account.

Until recently, you needed to use the command line tool for this step, but starting with R2014b we can just use MATLAB’s Current Folder window. No more Git commands like `git init`, `git status`, `git add`, or `git commit`!

Open your copy of MATLAB and create an empty folder. Right-clicking the empty space in the Current Folder window to bring up a contextual menu, and select *Source Control > Manage Files*.

This will open a new dialog box: *Manage Files Using Source Control*.

- For
*Select control integration*, choose*Git* - For
*Repository path*, click*Change*

You now see a new dialog box: *Select a Repository*. Copy and paste the URL of the remote repo you just created. You can find the URL in the right sidebar of your new repo on GitHub.

You choose either SSH or HTTPS depending on how you setup your authentication on GitHub.

Click *Validate*. You may be asked for your login password for authentication. You can close the dialog box when your path is validated.

Back in *Manage Files* dialog box, the *sandbox* should be already set to your current folder. All you need to do now is hit *Retrieve*.

You have now successfully cloned the remote repo to your local drive. Check your Current Folder window. You should see just one file – README.md, but with a green circle next to it. This is just a text file but you can apply wiki-like syntax called markdown to make it appear like a regular web page when viewed on GitHub. README serves as the front page of your repo on GitHub.

Let’s add a new MATLAB script file *helloworld.m*. It will appear with a blank circle – it means it is not added to Git source control yet. To add it to Git, right-click on the file and select *Add to Git*. The empty circle changes to “+” symbol. When you modify a file already under source control, the symbol becomes a blue square.

You can continue editing files as you like, but at some point, you want to take a snapshot of the edits you made. That’s when you do a commit. You can select any empty space in the Current Folder window to bring up the contextual menu and select *Commit to Git Repository*. This will bring up a dialog box where you can add your comment about the changes you made since the last commit. Comments will be helpful to keep to track of your changes and revert back to earlier commits if necessary.

When you commit, the snapshot is saved in the local repo, but it is also a good idea to mirror the changes to the remote repo as well. To do so, bring up the contextual menu by right-clicking an empty space in the Current Folder window and select *Push*. That will push your changes to the remote repo. You may need to enter your password.

The real power of source control comes from the ability to create multiple branches from your project. By default, you have a branch called “master” in your repo. You can create a new branch from the master branch, makes changes, and then merge those changes back to the master. This mechanism is used for working on new experimental features without affecting the working code on the master. You can branch and merge in MATLAB but the details are beyond the scope of this post.

If you have been curious about Git but put off by its complicated command line interface, Git integration in R2014b makes Git much more pleasant and approachable. I hope this quick introduction motivates you to take advantage of this new feature. When you do, please don’t forget to post your project to the File Exchange. To learn more about Git, it is actually helpful to start with reviewing the underlying concept about how Git works.

]]>- Connect MATLAB Mobile to your computer with the MATLAB Connector. For more information on how to do this, refer to the getting started instructions here. This feature is only supported on MATLAB R2014a and later, so make sure you are on a compatible version.
- Install the MATLAB Support Package for Android Sensors. Choose
*Add-ons*from the MATLAB Toolstrip, and then choose*Get Hardware Support Packages*. This will open the support package installer. Choose*Android Sensors*from the list and follow the instructions. - To establish communication between the sensors on your device and MATLAB, create a mobiledev object, as follows:

m = mobiledev;

Once you have completed the 3 steps from the above section, go to MATLAB Mobile and turn on the accelerometer. You should see something akin to this:

You can also enable the sensor directly from MATLAB, by executing the following command:

m.AccelerationSensorEnabled = 1;

Did you notice the enabled

m.Logging = 1;You can verify this in MATLAB, note the Current Sensor Values in the result:

m

m = mobiledev with properties: Connected: 1 Logging: 1 InitialTimestamp: '02-Oct-2014 21:53:26.707' AccelerationSensorEnabled: 1 (20 Logged values) AngularVelocitySensorEnabled: 0 MagneticSensorEnabled: 0 OrientationSensorEnabled: 0 PositionSensorEnabled: 0 Current Sensor Values: Acceleration: [0.2631 5.9226 8.1850] (m/s^2)

Walk around your campus/home/floor with your device. Once you are satisfied, stop sending this data to MATLAB. You can either tap on the

m.Logging = 0;To retrieve the data, use the accellog variable:

[a, t] = accellog(m);

Once you have retrieved the logged acceleration data, you can plot it in MATLAB:

plot(t, a); legend('X', 'Y', 'Z'); xlabel('Relative time (s)'); ylabel('Acceleration (m/s^2)');

Calculate the magnitude to convert your X, Y and Z vectors to scalar values. Then, plot it.

x = a(:,1); y = a(:,2); z = a(:,3); % Calculate and plot magnitude. mag = sqrt(sum(x.^2 + y.^2 + z.^2, 2)); plot(t, mag); xlabel('Time (s)'); ylabel('Acceleration (m/s^2)');

To remove constant effects such as gravity, you can subtract the mean from this data.

% Accounting for gravity. magNoG = mag - mean(mag); % Plot magnitude. plot(t, magNoG); xlabel('Time (s)'); ylabel('Acceleration (m/s^2)');

The plotted data is now centered on zero, and shows peaks which correspond to a step taken while walking.

To determine the number of steps taken, you can use to FINDPEAKS function from Signal Processing Toolbox. In this example, we are treating only peaks with a minimum height above one standard deviation as a step. This threshold should be tuned experimentally to match a person’s level of movement while walking, hardness of floor surfaces etc.

% Use FINDPEAKS to determine the local maxima. minPeakHeight = std(magNoG); [pks, locs] = findpeaks(magNoG, 'MINPEAKHEIGHT', minPeakHeight);The number of steps taken is the number of peaks:

numSteps = numel(pks)

numSteps = 15Finally, you can also identify these locations on your plot of acceleration magnitude data:

hold on; % Place a red marker on the locations that correspond to peaks. plot(t(locs), pks, 'r', 'Marker', 'v', 'LineStyle', 'none'); title('Counting Steps'); xlabel('Time (s)'); ylabel('Acceleration Magnitude, Gravity Removed (m/s^2)'); hold off;

Once you are done, make sure you turn off the acceleration sensor and clear the mobiledev object.

```
m.AccelerationSensorEnabled = 0;
clear m;
```

- Acquiring Data from Sensors - MATLAB Mobile
- Android Sensor Support from MATLAB
- Documentation for mobiledev Object

Let us know by leaving a comment below.

Published with MATLAB® R2014b

]]>