Today's post is by Sean de Wolski and me. We spent a few days this week at training for the start of the year, held in a Sheraton hotel.
At the completion of our stay, we each got a note in our rooms. Here is a sample.
Actually we read this in with the App from the Image Acquisition Toolbox and cropped the image (imcrop from Image Processing Toolbox) to just display the mysterious binary text.
Using functionality from the Computer Vision System Toolbox, we then converted our binary text using the ocr (Optical Character Recognition) function.
Icrop = imread('Sher.png');
bintext = ocr(Icrop)
str = bintext.Words;
str = [str{:}];
bintext = ocrText with properties: Text: '010101000110100001100001011011100110101100100...' CharacterBoundingBoxes: [879x4 double] CharacterConfidences: [879x1 single] Words: {14x1 cell} WordBoundingBoxes: [14x4 double] WordConfidences: [14x1 single]
Next guess was how the text was encoded  4 bits, 8 bits? We figured 8 was most likely.
newstr = reshape(str',8,[])';
Let's look at the first few lines.
snip = newstr(1:3,:)
snip = 01010100 01101000 01100001
Could these be text?
bin2dec(snip)
ans = 84 104 97
These could be ASCII characters!
Have you ever misread an O for a 0? Or an l for a 1? We have! So we decided to program defensively for those and any other weird swaps that might confuse the OCR.
newstr(newstr == 'O') = '0'; newstr(newstr ~='0') = '1';
Finally let's write out the result.
message = char(bin2dec(newstr))'; disp(message(1:74)) disp(message(75:end))
Thank you for your continued loyalty to Sheraton Boston. It is a pleasure to host Mathworks year after year.
Get
the MATLAB code
Published with MATLAB® R2014b
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.
]]>I promised earlier to post about divergent colormaps. A divergent colormap is usually constructed by concatenating two colormaps together that have different color schemes.
Here is an example of a divergent colormap from www.ColorBrewer.org by Cynthia A. Brewer, Geography, Pennsylvania State University.
rgb = [ ... 94 79 162 50 136 189 102 194 165 171 221 164 230 245 152 255 255 191 254 224 139 253 174 97 244 109 67 213 62 79 158 1 66 ] / 255; b = repmat(linspace(0,1,200),20,1); imshow(b,[],'InitialMagnification','fit') colormap(rgb)
This colormap has the brightest color in this middle, which is typical for divergent colormaps.
A divergent colormap is used to compare data values to a reference value in a way that visually highlights whether values are above or below the reference.
Let me show you an example using sea surface temperatures. I learned today that 28 degrees Celsius has a specific meaning in hurricane modeling and prediction. It is considered to be the minimum temperature for hurricane formation.
Using the NASA EOSDIS Reverb tool, I found and downloaded sea surface temperature from the Aqua satellite for August 1, 2011. The data arrived as a netCDF file.
amsre = ncread('20110801AMSREREMSSL2P_GRIDDED_25amsre_20110801rtv01.nc',... 'sea_surface_temperature'); size(amsre)
ans = 1440 720 2
The two planes contain data from two different orbits.
orbit1 = flipud(amsre(:,:,1)'); orbit2 = flipud(amsre(:,:,2)'); imshow(orbit1,[],'InitialMagnification','fit') title('Orbit 1')
imshow(orbit2,[],'InitialMagnification','fit') title('Orbit 2')
The black regions are NaNs representing missing data. Combine the two orbits by taking the maximum.
sst = max(orbit1,orbit2); imshow(sst,[],'InitialMagnification','fit')
Now grab a subset of the data that includes a portion of the Atlantic Ocean.
ssta = sst(180:400,330:700); imshow(ssta,[],'InitialMagnification','fit')
Apply the divergent colormap.
colormap(rgb) colorbar
OK, I see one problem right off the bat. It looks like our data is in Kelvin. Let's fix that.
ssta = ssta  273.15; imshow(ssta,[],'InitialMagnification','fit') colormap(rgb) colorbar
Now I want to make use of our divergent colormap to highlight the temperature of interest, 28 degrees. To do that, I'll use the CLim property of the axes object. The CLim property is a twoelement vector. The first element tells us the data value that maps to the lowest colormap color, and the second element tells us the data value that maps to the highest colormap color. MATLAB computes these values automatically based on the range of the data.
get(gca,'CLim')
ans = 1.9500 33.6000
To get the middle colormap color to represent 28 degrees, we need to pick our two CLim values so that they are equally far away from 28 degrees. After some experimentation, I settled on a range from 22 degrees to 34 degrees.
set(gca,'CLim',[22 34])
That's looking better, but it isn't helpful to have the missing data displayed using the bottom color of the colormap. I will fix that problem using this procedure.
Here are just the first two steps.
mask = isnan(ssta); graymask = 0.5 * mask; graymask = repmat(graymask,1,1,3); hold on h = imshow(graymask); hold off
The missing values are displayed as gray, but we can't see the sea surface temperatures. Fix that by using the NaNmask as the overlay image's AlphaData.
h.AlphaData = mask; % This syntax requires MATLAB R2014b.
Temperatures close to the threshold are displayed using yellow, the brightest color in the colormap and the one right in the middle. Temperatures above the threshold are displayed using oranges and reds, while temperatures below are displayed using greens and blues.
Before I go, I want to give a shout to the new Developer Zone blog that just started today. If you are into advanced software developer topics, I encourage you to head over there and take a look.
Get
the MATLAB code
Published with MATLAB® R2014b
If $n$ is odd, replace $n$ by $3n+1$, if not, replace $n$ by $n/2$. Repeat. A famous conjecture made by Lothar Collatz is that no matter what value of $n$ is chosen to start, the process eventually terminates at $n=1$. Do not expect a proof, or a counterexample, in this blog.
You might want to download function threenplus1.m from my textbook Numerical Computing with MATLAB. This not only produces the graphs in this post, it also has a pair of uicontrols that allow you to increase and decrease the value of n.
Here is the core of threenplus1. The generated sequence is collected in a vector y. We don't know ahead of time how long y is going to be. In fact, that's the point of the computation. Does the while loop terminate and, if so, how long is y?. The ability of MATLAB to grow a vector an element at a time comes in very handy here.
dbtype 44:52 threenplus1
44 y = n; 45 while n > 1 46 if rem(n,2)==0 47 n = n/2; 48 else 49 n = 3*n+1; 50 end 51 y = [y n]; 52 end
A good example is provided by starting with $n = 7$. Here is the sequence collected in y. Its length is 17 and its max is 52.
7 22 11 34 17 52 26 13 40 20 10 5 16 8 4 2 1
Here is the plot. The yaxis is logarithmic in anticipation of some fairly large values in later runs. Notice that as soon as the sequence hits a power of 2, it plunges to termination.
It is interesting to see how the graphs produced by threenplus1 change from one value of $n$ to the next. Here are the first ten elements of y for $n$ = 108, 109, and 110.
108 54 27 82 41 124 62 31 94 47 . . . 109 328 164 82 41 124 62 31 94 47 . . . 110 55 166 83 250 125 376 188 94 47 . . .
After the first eight steps all three sequences reach the same value, so they are identical after that. The length is 114 and the max is 9232. Here are the three graphs, superimposed on one another. The similarity between the three sequences is very apparent when you use the increment n buttons provided by threenplus1.
Prime numbers of the form $n = 2^p1$ where $p$ itself is prime are known as Mersenne primes. But that's a whole 'nother story. For now, suffice it to say that such numbers are candidates to produce 3n+1 sequences with a high maximum value. Here is $n = 2^{13}1$. The sequence reaches 6,810,136 before it eventually terminates in 159 steps.
I have to admit that this triggers a bug in some versions of MATLAB. If you see specious characters in the formatting of the legend on the yaxis, insert this line near the end of threenplus1.
set(gca,'yticklabel',int2str(get(gca,'ytick')'))
Check out the two references below about orbits associated with the Collatz conjecture.
When I Google "three n plus 1", the first link is to an excellent Wikipedia article on the "Collatz conjecture". According to Wikipedia, the famous German mathematican Lothar Collatz first made the conjecture, in 1937, that the process terminates for any starting value. I am already familiar with Collatz's name because of his work in numerical analysis and approximation theory, including finite difference methods for differential operators. My mentor John Todd had some contact with Collatz shortly after the end of World War II that could be the subject of a blog post some day.
Many others have also been interested in the 3n+1 problem. Stanislaw Ulam (see my previous blog post), Shiszuo Kakutani, Sir Bryn Thwaites, Helmut Hasse, Paul Erdos, John Horton Conway, Douglas Hofstadter, and Richard Guy. But Wikipedia quotes Paul Erdos as saying about the Collatz conjecture: "Mathematics may not be ready for its solution."
Wikipedia, Collatz Conjecture, <http://en.wikipedia.org/wiki/Collatz_conjecture>
Jason Davies, Collatz Conjecture Orbits, <http://www.jasondavies.com/collatzgraph>
xkcd, Collatz Conjecture Orbits, <http://xkcd.com/710>
Get
the MATLAB code
Published with MATLAB® R2014b
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 6by7 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.
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 fourinarow 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 gamespecific 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 Tshape, or an Sshape. Not all of these variants are interesting as games. I built square and Tshape 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 playtesting 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 NietoCastanon 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.
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 clicktomove 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.
]]>Simulation Target
If you click play in the following model:
You will notice that a folder named slprj with a lot of code and 4 MEXfiles will be created:
Simulink Preferences
If you do not want the slprj and MEXfiles to be generated in your current directory, it is possible to go in the Simulink Preferences and specify a Simulation Cache Folder and a Code Generation Folder to be used instead of the local directory.
Overriding the Simulink Preferences
Now, let's say I need to often switch between different projects, and each of those projects need to use different folders for code generation. In that case, it is possible to use Simulink.fileGenControl. For example, if I include Simulink Project in my workflow, the startup shortcut for my project could look like:
project = simulinkproject; projectRoot = project.RootFolder; myCacheFolder = fullfile(projectRoot, 'work'); myCodeFolder = fullfile(projectRoot, 'code'); Simulink.fileGenControl('set',... 'CacheFolder', myCacheFolder,... 'CodeGenFolder', myCodeFolder,... 'createDir', true)
That way, all the code and MEXfiles generated by the model will go in those folders and the current directory will remain clean.
Once you are done working on this project, you can restore the values stored in the Simulink Preferences using this line:
Simulink.fileGenControl('reset')
Don't Clean Up: Faster Initialization
Simulink checks the cached files at initialization to make sure they are up to date and match the model. This prevents regenerating files, and results in faster initialization. If you regularly work with the same models, keeping these derived files around can help you save time when switching between projects.
Now it's your turn
Let us know how you use this feature by leaving a comment here
]]>Sean's pick this week is A thin MATLAB wrapper for the Git source control system by Manu Raghavan.
I just began working on my first collaborative project where we're using Git as the Source Control system.
Starting in MATLAB R2014b, MATLAB supports Git and SVN source control directly from the Current Folder browser.
The Current Folder has support for a common subset of Git functionality. However, some of the more advanced operations, like merging, are not supported. This is where Manu's submission comes in handy. It provides a simple function that wraps around the Git OS command to give you all of the functionality in a similar mechanism to those who have used command line Git before.
Let's take a look at it working when I fixed a bug in my generateReport function. The blue square indicates that the file has changed since it was last commited.
And then to commit the fix:
All green, good to go!
As an aside  Additionally, I learned about the more() function from looking through Manu's code as well!
Give it a try and let us know what you think here or leave a comment for Manu.
Also, what do you think of the new changes to MATLAB Central profile? For example, here is my colleague, Brett's page:
Please let us know your feelings on these changes as well below too!
Get
the MATLAB code
Published with MATLAB® R2014b
Today I’d like to introduce guest blogger Ken Atwell who works for the MATLAB Development team here at MathWorks. Today, Ken will be discussing with you the MapReduce programming technique now available in the R2014b release of MATLAB. MapReduce provides a way to process large amounts of filebased data on a single computer in MATLAB. For very large data sets, the same MATLAB code written using MapReduce can also be run on the "big data" platform, Hadoop®.
The data set we will be using consists of records containing flight performance metrics for USA domestic airlines for the years 1987 through 2008. Each year consists of a separate file. If you have experimented with "big data" before, you may already be familiar with this data set. The full data set can be downloaded from here. A small subset of the data set, airlinesmall.csv, is also included with MATLAB® to allow you to run this and other examples without downloading the entire data set.
MapReduce is a programming technique that is used to "divide and conquer" big data problems. In MATLAB, the mapreduce function requires three input arguments:
This is an oversimplification to some extent, since the output of a call to the map function can be shuffled and combined in interesting ways before being passed to the reduce function. This will be examined later.
Let's look at a straightforward example to illustrate how mapreduce is used. In this example we want to find the longest flight time out of all the flights recorded in the entire airline data set. To do this we will:
datastore is used to access collections of tabular text files stored on a local disk, or the Hadoop® Distributed File System (HDFS™). It is also the mechanism for providing data in a chunkwise manner to the map calls when using mapreduce. A previous blog post explained how datastore works and how it is used for reading collections of data that are too large to fit in memory.
Let's create our datastore and first preview the data set. This allows us to peek into the data set, and identify the format of the data and the columns that contain the data we are interested in. The preview will normally provide a small chunk of data that contains all the columns present in the data set, although I will only show a handful of columns to make our article more readable.
ds = datastore('airlinesmall.csv', 'DatastoreType', 'tabulartext', ... 'TreatAsMissing', 'NA'); data = preview(ds); data(:,[1,2,3,9,12])
ans = Year Month DayofMonth UniqueCarrier ActualElapsedTime ____ _____ __________ _____________ _________________ 1987 10 21 'PS' 53 1987 10 26 'PS' 63 1987 10 23 'PS' 83 1987 10 23 'PS' 59 1987 10 22 'PS' 77 1987 10 28 'PS' 61 1987 10 8 'PS' 84 1987 10 10 'PS' 155
For the computation we are looking to perform, we only need to look at the "ActualElapsedTime" column: it contains the actual flight time. Let's configure our datastore to only provide this one column to our map function.
ds.SelectedVariableNames = {'ActualElapsedTime'};
Now we will write our map function (MaxTimeMapper.m). Three input arguments will be provided to our map function by mapreduce:
Our map function finds the maximum value in the table 'data', and saves a single key ('MaxElapsedTime') and corresponding value to the intermediate data storage object. We will now go ahead and save the following map function (MaxTimeMapper.m) to our current folder.
function MaxTimeMapper(data, ~, intermediateValuesOut) maxTime = max(data{:,:}); add(intermediateValuesOut, 'MaxElapsedTime', maxTime); end
Next, we create the reduce function (MaxTimeReducer.m). Three input arguments will also be provided to our reduce function by mapreduce:
Our reduce function will receive a list of values associated with the 'MaxElapsedTime' key which were generated by the calls to our map function. The reduce function will iterate through these values to find the maximum. We create the following reduce function (MaxTimeReducer.m) and save it to our current folder.
function MaxTimeReducer(~, intermediateValuesIn, finalValuesOut) maxElapsedTime = inf; while hasnext(intermediateValuesIn) maxElapsedTime = max(maxElapsedTime, getnext(intermediateValuesIn)); end add(finalValuesOut, 'MaxElapsedTime', maxElapsedTime); end
Once the map and reduce functions are written and saved in our current folder, we can call mapreduce and reference the datastore, map function, and reduce function to run our calculations on our data. The readall function is used here to display the results of the MapReduce algorithm.
result = mapreduce(ds, @MaxTimeMapper, @MaxTimeReducer); readall(result)
Parallel mapreduce execution on the local cluster: ******************************** * MAPREDUCE PROGRESS * ******************************** Map 0% Reduce 0% Map 100% Reduce 100% ans = Key Value ________________ ______ 'MaxElapsedTime' [1650]
If Parallel Computing Toolbox is available, MATLAB will automatically start a pool and parallelize execution of the map functions. Since the number of calls to the map function by mapreduce corresponds to the number of chunks in the datastore, parallel execution of map calls may speed up overall execution.
The use of keys is an important feature of mapreduce. Each call to the map function can add intermediate results to one or more named "buckets", called keys.
If the map function adds values to multiple keys, this leads to multiple calls to the reduce function, with each reduce call working on only one key's intermediate values. The mapreduce function automatically manages this data movement between the map and reduce phases of the algorithm.
The behavior of the map function in this case is more complex and will illustrate the benefits to using multiple keys to store your intermediate data. For every airline found in the input data, use the add function to add a vector of values. This vector is a count of the number of flights for that carrier for each day in the 21+ years of data. The carrier code is the key for this vector of values. This ensures that all of the data for each carrier will be grouped together when mapreduce passes it to the reduce function.
Here's our new map function (CountFlightsMapper.m).
function CountFlightsMapper(data, ~, intermediateValuesOut)
dayNumber = days((datetime(data.Year, data.Month, data.DayofMonth)  datetime(1987,10,1)))+1;
daysSinceEpoch = days(datetime(2008,12,31)  datetime(19875,10,1))+1;
[airlineName, ~, airlineIndex] = unique(data.UniqueCarrier);
for i = 1:numel(airlineName) dayTotals = accumarray(dayNumber(airlineIndex==i), 1, [daysSinceEpoch, 1]); add(intermediateValuesOut, airlineName{i}, dayTotals); end end
The reduce function is less complex. It simply iterates over the intermediate values and adds the vectors together. At completion, it outputs the values in this accumulated vector. Note that the reduce function does not need to sort or examine the intemediateKeysIn values; each call to the reduce function by mapreduce only passes the values for one airline.
Here's our new reduce function (CountFlightsReducer.m).
function CountFlightsReducer(intermediateKeysIn, intermediateValuesIn, finalValuesOut)
daysSinceEpoch = days(datetime(2008,12,31)  datetime(1987,10,1))+1;
dayArray = zeros(daysSinceEpoch, 1);
while hasnext(intermediateValuesIn) dayArray = dayArray + getnext(intermediateValuesIn); end add(finalValuesOut, intermediateKeysIn, dayArray); end
To run this new analysis, just reset the datastore and select the active variables of interest. In this case we want the date (year, month, day) and the airline carrier name. Once the map and reduce functions are written and saved in our current folder, mapreduce is called referencing the updated datastore, map function, and reduce function.
reset(ds); ds.SelectedVariableNames = {'Year', 'Month', 'DayofMonth', 'UniqueCarrier'}; result = mapreduce(ds, @CountFlightsMapper, @CountFlightsReducer); result = result.readall();
Parallel mapreduce execution on the local cluster: ******************************** * MAPREDUCE PROGRESS * ******************************** Map 0% Reduce 0% Map 100% Reduce 50% Map 100% Reduce 100%
Up until this point we have been only analyzing the sample data set (airlinesmall.csv). To look at the number of flights each day over the full dataset, let's load the results of running our new MapReduce algorithm on the entire data set.
load airlineResults
Before looking at the number of flights per day for the top 7 carriers, let us first apply a filter to the data to smooth out the effects of weekend travel. This would otherwise clutter the visualization.
lines = result.Value; lines = horzcat(lines{:}); [~,sortOrder] = sort(sum(lines), 'descend'); lines = lines(:,sortOrder(1:7)); result = result(sortOrder(1:7),:); lines(lines==0) = nan; for carrier=1:size(lines,2) lines(:,carrier) = filter(repmat(1/7, [7 1]), 1, lines(:,carrier)); end
figure('Position',[1 1 800 800]); plot(datetime(1987,10,1):caldays(1):datetime(2008,12,31),lines) title ('Domestic airline flights per day per carrier') xlabel('Date') ylabel('Flights per day (7day moving average)') try carrierList = readtable('carriers.csv.txt', 'ReadRowNames', true); result.Key = cellfun(@(x) carrierList(x, :).Description, result.Key); catch end legend(result.Key, 'Location', 'SouthOutside')
The interesting characteristic highlighted by the plot is the growth of Southwest Airlines (WN) during this time period.
The code we have just created was done so on our desktop computer and it allows us to analyze data that normally would not fit into the memory of our machine. But what if our data was stored on the "big data" platform, Hadoop® and the data set is too large to reasonably transfer to our desktop computer?
Using Parallel Computing Toolbox along with MATLAB® Distributed Computing Server™, this same code can be run on a remote Hadoop cluster. The map and reduce functions would remain unchanged, but two configuration changes are required:
After these configuration changes, the algorithm can then be run on Hadoop.
MATLAB Compiler can be used to deploy applications created using MapReduce, and it can also generate Hadoop specific components for use in production settings.
Learn more about MapReduce in MATLAB with our algorithm examples and let us know how you might see MapReduce being used to analyze your big data here.
Get
the MATLAB code
Published with MATLAB® R2014b
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 gameplaying 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 gameplaying program is the lookahead. 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 lookahead 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. So 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 notsoclever 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 objectoriented programming, you can write a generic gameplaying harness. Then you just need to plug in some code that knows a few rules, and presto! You’ve got an instant gameplaying 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 playtesting 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: MonteCarloGames. Here is the TicTacToe class, and here is the botMoves function.
NEXT WEEK
This is the first of a twopart 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!
Jiro's pick this week is num2words by Stephen Cobeldick.
Just recently, I witnessed a typo on a check payment. The amount was for "$19050" (unfortunately, it wasn't for me) but the written text read "One thousand nine hundred fifty dollars". Many people might wonder why someone would make such a simple mistake. Did the person miscount the number of zeroes?
I think I have an idea. This was written in Japan. In Japan, and in many other Asian countries, large numbers are grouped in myriads (every 10000). is called man (万), is called oku (億), is called cho (兆). However in Western counting, large numbers are grouped in 1000s  thousand, million, billion, trillion, etc. So, it's likely that the person writing the check wrote "19050", read it as "1 man 9050", and incorrectly translated to "1 thousand 9 hundred 50". Maybe...
Anyway, all of this could have been avoided if this person had access to Stephen's num2words! Just let his function spell out the number in words.
num2words(19050)
ans = nineteen thousand and fifty
This is the simplest use case, but num2words comes with a wealth of options to customize how the results are returned.
num2words(19050,'type','cheque')
ans = nineteen thousand and fifty dollars only
Stephen includes a detailed help in the function, so you can learn about all the options just by typing
help num2words
What gave me a little chuckle was that when I looked inside his function, I could see that he made sure that his function can count up to Uncentillion which is . So if you need to spell out a number larger than , you might be out of luck. :)
Comments
Give this a try, and let us know what you think here or leave a comment for Stephen.
Get
the MATLAB code
Published with MATLAB® R2014b
The prime spiral was discovered by Stanislaw Ulam in 1963, and featured on the cover of Scientific American in March, 1964.
Stanislaw Ulam (19091984) was an eminent PolishAmerican mathematician who spent most of his career at Los Alamos, but who also had connections with several American universities, especially the University of Colorado. His many illustrious colleagues included von Neumann, Teller, Metropolis, Erdos, Richtmyer, Rota, and Kac.
Ulam had a wide range of mathematical interests. He was deeply involved in the early development of computers at Los Alamos. He is one of the primary developers of the Monte Carlo method. Later in his career he made contributions to mathematical biology and medicine. He wrote popular books, including "Adventures of a Mathematician" and "Analogies Between Analogies".
As the story goes, Ulam was bored during a seminar in 1963 and began to doodle in his notebook. He numbered the integer lattice points in the plane in a spiral fashion and then highlighted the prime numbers. He was surprised to see that with this numbering the primes were not distributed randomly throughout the plane, but appeared to fall on diagonal segments.
Here is the 10by10 prime spiral, showing the spiral numbering and the 25 primes among the first 100 integers. The most prominent diagonal segment includes the primes 19, 7, 23, 47, and 79.
Ulam's numbering at any point $(i,j)$ in the plane is a piecewise quadratic function of $i$ and $j$. For example, the values along the upper and lower halves of the antidiagonal, or their offsets by one, (depending upon whether the size is odd or even), are perfect squares. So there are certainly not any primes along these diagonals.
This reminds me of the quadratics that generate primes. The most famous is Euler's polynomial,
$$ p(n) = n^2n+41 $$
This generates primes for $n$ from 1 through 40, but, obviously, not for $n$ equal to 41.
n = 1:41; p = n.^2n+41; isprime(p)
ans = Columns 1 through 13 1 1 1 1 1 1 1 1 1 1 1 1 1 Columns 14 through 26 1 1 1 1 1 1 1 1 1 1 1 1 1 Columns 27 through 39 1 1 1 1 1 1 1 1 1 1 1 1 1 Columns 40 through 41 1 0
The Euler primes themselves don't show up in the Ulam prime spiral. But each of the diagonal segments corresponds to a quadratic polynomial, although probably one that generates fewer than 40 primes.
The MATLAB demos directory has a spiral function that generates Ulam's numbering. I'm afraid that the code is ancient, and pretty bad. It's not too hard to write better code. This generates successively larger spirals by rotating one of a given size 180 degrees and then appending another column at the right and another row at the bottom to produce the next larger size.
type improved_spiral
function S = improved_spiral(n) % SPIRAL SPIRAL(n) is an nbyn matrix with elements ranging % from 1 to n^2 in a rectangular spiral pattern. S = []; for m = 1:n S = rot90(S,2); S(m,m) = 0; p = m^2m+1; v = (m1:1:0); S(:,m) = pv'; S(m,:) = p+v; end
Once we have the spiral numbering, the prime spiral is our old friend, the spy plot.
type primespiral_core.m
function primespiral(n) % PRIMESPIRAL PRIMESPIRAL(n) is Ulam's prime spiral. P = zeros(n,n); P(primes(n^2)) = 1; S = improved_spiral(n); P = P(S); spy(P)
Here are prime spirals of three different sizes.
The function primespiral in my textbook Numerical Computing with MATLAB includes a provision for starting the numbering at a value larger than one. Animating this feature emphasizes the diagonal nature of the prime locations. The title on this plot is the starting value of the numbering.
Get
the MATLAB code
Published with MATLAB® R2014b
And so...another year over, a new one just begun.
2014 saw some very big changes. MATLAB has changedsubtly. If you haven't tried R2014b yet, try itwe think you'll like it. The new MATLAB Graphics systemlongawaited! is live, and while there may be some growing pains, we're excited about the new things it will enable us to do. In this forum, we welcomed Sean as a lead blogger, and Brett stepped aside to pursue other projects. (He's still around, contributing now and then ;) !)
And once again, the year has seen big growth of MATLAB Central, and lots of big additions to the File Exchange. So without further ado, here is the retrospective of 2014's Picks of the Week.
You can read the original blog posts for each of these picks by clicking on the names of the files and following the links to the "Pick of the Week" tag. As always, your comments are welcome.
Get
the MATLAB code
Published with MATLAB® R2014b
Title/Author/Summary/Notes  Image 

__________________________________________________________________________  
Smooth 3D bezier curves with implicit control points
Uses Hobby's algorithm to plot smooth curves in 3D through specified control points


__________________________________________________________________________  
Sphere Fit (least squared)
Fits a sphere to a set of noisy data. Does not require a wide arc or many points.


__________________________________________________________________________  
Texas Hold'Em Poker Analyzer
Analyzes the win frequency and expected return of a hand in a poker game called Texas Hold'Em.


__________________________________________________________________________  
URL Filter
Scrape one or more numbers off of a web page


__________________________________________________________________________  
International Flight Route Planning Simulator using Aerospace Blockset
Example that shows how you can visualize airline routes on the globe.


__________________________________________________________________________  
Bellagio Fountain Simulation
Simulation of the Bellagio Fountain (100 nozzles)


__________________________________________________________________________  
Polygon simplification
Simplifies a closed polygon by reducing the number of vertices to a specified number


__________________________________________________________________________  
Worker Object Wrapper
Simplifies managing resources such as large data within PARFOR loops and SPMD blocks


__________________________________________________________________________  
Draw a Polar Dendrogram
Draws a polar dendrogram.


__________________________________________________________________________  
Convert seconds to human readable string(s)
Convert an amount of seconds to hours/days/weeks/etc.


__________________________________________________________________________  
Straightforward COPY and PASTE functions
Allows very simple manual exchange with other applications through the clipboard.


__________________________________________________________________________  
Packing Santa's Sleigh
This is the solution by team alfnie to the Mathworks/Kaggle Packing Santa's Sleigh competition


__________________________________________________________________________  
2048 MATLAB Edition
This is a MATLAB implementation of the 2048 game


__________________________________________________________________________  
Communications System Toolbox Support Package for RTLSDR Radio
MathWorks Communications System Toolbox Team Design and prototype softwaredefined radio (SDR) systems using RTLSDR with MATLAB and Simulink


__________________________________________________________________________  
kml2struct
Import a .kml file as a series of shapefile structs. Like kml_shapefile, but more stable.


__________________________________________________________________________  
BOT
Block diagram Optimization Tool


__________________________________________________________________________  
ISSTracker V1.0.0
International Space Station Live Tracking GUI


__________________________________________________________________________  
Wing Designer
Wing Designer computes aircraft performance measures from wing and engine parameters.


__________________________________________________________________________  
Automatic enveloping, variance change and activity detection with Hilbert Transform
Automatic Signal Segmentation and activity detection with Hilbert Transform and smoothing.


__________________________________________________________________________  
cspy.m
cspy: a colorcoded version of sparse matrix viewer spy


__________________________________________________________________________  
MultipleColon
multiple colonintervals


__________________________________________________________________________  
Configurable Simulink Model for DCDC Converters with PWM PI Control
A Simulink model configurable to buck, boost and buckboost DCDC converters with PWM PI control


__________________________________________________________________________  
Useful Figure Management Utilities
3 Useful figure management utilities.


__________________________________________________________________________  
CMEX Programming Tutorial Examples
This submission contains several C scripts that help as an introduction to CMEX programming.


__________________________________________________________________________  
module  encapsulate MATLAB package into a name space module
Modules are handy way to access packages in MATLAB.


__________________________________________________________________________  
legappend
Append new entries to an existing legend


__________________________________________________________________________  
parTicToc
This is a utility for timing PARFOR loops.


__________________________________________________________________________  
ThingSpeak support from MATLAB
MathWorks Internet of Things Team Prototype Internet of Things applications using ThingSpeak and MATLAB


__________________________________________________________________________  
PianoTuner
Graphical aid for tuning the middle octave of a piano


__________________________________________________________________________  
Reference Creator
Find intext citations and create the reference list!


__________________________________________________________________________  
Image Blur Metric
Measure the blur level of still image.


__________________________________________________________________________  
Blind image quality assessment through anisotropy
This measure discriminates blur and Gaussian noise.


__________________________________________________________________________  
Noise Level Estimation from a Single Image
It can precisely estimate noise level from a single image.


__________________________________________________________________________  
Simscape Language: Nonlinear rotational spring (Torque == kx^2)
Simscape language example showing how to use a MATLAB function in the equations section.


__________________________________________________________________________  
SimRF Models of Analog Devices RF Agile Transceivers
Simulate and verify agile RF transceivers and predict impact of RF imperfections using test signals


__________________________________________________________________________  
Colormouse
Mouse control of the color axis


__________________________________________________________________________  
makeinstall
Create a single install file for your toolbox distribution.


__________________________________________________________________________  
Makeinstall Technology
The technology behind the Makeinstall tool.


__________________________________________________________________________  
QuadSim
Parameter driven Simulink model for quadcopter simulation and control system design


__________________________________________________________________________  
Display progress, pause or stop a timeconsuming loop
Allows you to display a loop’s progress, pause (and also change related code) or stop it.


__________________________________________________________________________  
Simulink Dual Gravity Drained Tank
Process control simulation tutorial for PID tuning


__________________________________________________________________________  
votebar
function to plot 'electionlike' 3d bars


__________________________________________________________________________  
Fast Noise Estimation in Images
Estimate the standard deviation of the noise in a grayscale image.


__________________________________________________________________________  
Round with significant digits
Rounds towards the nearest number with N significant digits.


__________________________________________________________________________  
NaturalOrder Row Sort
Naturalorder row sort of a cell array of strings, with customizable numeric format.


__________________________________________________________________________  
NaturalOrder Filename Sort
Naturalorder sort of filenames or filepaths, with customizable numeric format.


__________________________________________________________________________  
Customizable NaturalOrder Sort
Naturalorder sort of a cell array of strings, with customizable numeric format.


__________________________________________________________________________  
Growth  Yet another Obfuscated MATLAB code
Obfuscated MATLAB code where the source looks like its output, with a functional programming demo.


__________________________________________________________________________  
Distribute figures
This function can distribute all open figures on different parts of the screen(s).


__________________________________________________________________________  
Koch curve
Interface for generating the Koch curve

Guest blogger, Kelly Luetkemeyer, who is a software developer at MathWorks, returns with an article on accessing RESTful web services using MATLAB. Kelly's previous articles included Tracking a Hurricane using Web Map Service and Visualizing the Gulf of Mexico Oil Slick using Web Map Service. Kelly was the lead developer of the new RESTful web service interface feature in MATLAB.
In R2104b, MATLAB introduced three new features to access RESTful web services:
To illustrate the use of these features, let's do a thought experiment and imagine that you are interested in applying for a software development job at MathWorks. This blog post illustrates how you can use RESTful web services to help in your job search. In particular it shows you how to:
The first step in your job search is to obtain a listing of all available jobs for your desired location. Job listings are available in an RSS feed, and you can use the function webread to read the feed. The RSS feed is an XML format. Specify 'xmldom' content type in order to return a Java Document Object Model for easy access to the data. Read the RSS feed and parse the data to create a listing of jobs of interest.
options = weboptions('ContentType','xmldom'); rssURL = 'http://www.mathworks.com/company/jobs/opportunities/rss.xml'; dom = webread(rssURL, options);
Retrieve all the items from the DOM that are tagged 'item'.
items = dom.getElementsByTagName('item');
itemsLength = items.getLength;
Create a table to contain the information from the RSS feed. Map the XML tag names to the table variable names by using two cell arrays. Initialize the table using cell2table.
tagNames = {'id', 'location', 'title', 'job_function', 'job_type'}; tableVariableNames = {'Reference', 'Location', 'Title', 'Department', 'Type'}; emptyItems = cell(itemsLength,length(tableVariableNames)); jobsTable = cell2table(emptyItems,'VariableNames',tableVariableNames);
Loop through each item and add the information to the table. Use zerobased indexing to obtain the node since the data is contained in a Java object. Use onebased indexing for the table array.
for k=0:itemsLength1 item = items.item(k); for n = 1:length(tagNames) node = item.getElementsByTagName(tagNames{n}); jobsTable{k+1,n} = {char(node.item(0).getTextContent)}; end end
Select your desired location.
myLocation = 'USMANatick';
rows = strcmpi(jobsTable.Location, myLocation);
natickJobs = jobsTable(rows,:);
natickJobs.Location = [];
Select your desired department.
myDepartment = 'software development';
rows = strcmpi(natickJobs.Department, myDepartment);
natickSoftwareJobs = natickJobs(rows,:);
natickSoftwareJobs.Department = [];
Display a table of possible jobs that match your request.
natickSoftwareJobs
natickSoftwareJobs = Reference Title Type _________ _________________________________________________________________________ ___________________ '14070' 'MATLAB Editor Software Engineer (Senior) (14070MCAR)' 'Experienced' '14064' 'Software Engineer  Graphical Platform (14064BWAL)' 'Experienced' '14049' 'Release Engineer (14049BHIL)' 'Experienced' '14028' 'Software Engineer (14028MCOL)' 'Experienced' '14027' 'Software Engineer (14027MCOL)' 'Experienced' '14023' 'Software Developer (14023MCOL)' 'Experienced' '13866' 'Controls Systems Engineer (13866MCAR)' 'Experienced' '13859' 'Software Engineer (13859MCOL)' 'Experienced' '13857' 'Web User Interface Developer (13857MCOL)' 'Experienced' '13856' 'Senior Software Engineer (13856MCOL)' 'Experienced' '13855' 'Software Engineer (13855MCOL)' 'Experienced' '13854' 'Software Engineer  GUI Components (13854MCOL)' 'Experienced' '13853' 'Build and Release Engineer  Perl / Shell (13853MCOL)' 'Experienced' '13851' 'Code Replacement Library Software Engineer (13851GMAR)' 'Experienced' '13848' 'Senior Performance Engineer (13848MCAR)' 'Experienced' '13793' 'Software Engineer  MATLAB Development Tools (13793MCAR)' 'Experienced' '13792' 'Software Engineer  Desktop / Toolbox Integration (13792MCAR)' 'Experienced' '13791' 'Autonomous Controls Systems Engineer (13791MCAR)' 'Experienced' '13789' 'Software Engineer  MATLAB Hardware (13789MCAR)' 'Experienced' '13779' 'Web Application Developer (13779JJUS)' 'Experienced' '13778' 'Software Engineer  Cloud Storage (13778JJUS)' 'Experienced' '13777' 'Senior Software Engineer  Web Services (13777JJUS)' 'Experienced' '13776' 'Senior Software Engineer  MATLAB Online (13776JJUS)' 'Experienced' '13775' 'Senior Technical Lead  MATLAB Web Services (13775JJUS)' 'Experienced' '13774' 'Java Software Engineer  Support Software Installer (13774JJUS)' 'Experienced' '13773' 'Java Software Engineer (13773JJUS)' 'Experienced' '13772' 'Mobile Software Engineer (13772JJUS)' 'Experienced' '13752' 'Software Engineer  Numerical Optimization (13752SMAR)' 'Experienced' '13751' 'Principal Software Engineer  Data Munging (13751SMAR)' 'Experienced' '13750' 'Senior Software Engineer  MATLAB Math (13750SMAR)' 'Experienced' '13747' 'Image Algorithms Software Engineer (13747SMAR)' 'Experienced' '13746' 'Software Engineer  Data Analysis (13746SMAR)' 'Experienced' '13745' 'Data Analysis Tools Developer (13745SMAR)' 'Experienced' '13743' 'Senior Software Engineer  Statistics and Machine Learning (13743SMAR)' 'Experienced' '13742' 'Senior MATLAB Web UI Engineer (13742SMAR)' 'Experienced' '13734' 'Senior Software Engineer  Application Deployment (13734MCAR)' 'Experienced' '13733' 'Software Engineer  Compiler Componentization (13733MCAR)' 'Experienced' '13732' 'UI Developer  Toolboxes & Addons (13732MCAR)' 'Experienced' '13730' 'Software Engineer  C+ (13730MCAR)' 'Experienced' '13729' 'Software Engineer  C++ (13729MCAR)' 'Experienced' '13689' 'Simulink Developer Engineer (13689KCAR)' 'Internships' '13668' 'CoOp: Design Automation Software Development (13668KCAR)' 'Internships' '13660' 'C++ Software Engineer Algorithms and Code Generation (13660GMAR)' 'Experienced' '13651' 'Compiler Engineer for FPGA Design (13651KCAR)' 'Internships' '13621' 'C++ Software Engineer  User Interfaces and Visualization (13621BWAL)' 'Experienced' '13613' 'MATLAB Web Applications Engineer  Intern (13613JJUS)' 'Internships' '13603' 'C++ Software Engineer Intern (13603KCAR)' 'Internships' '13594' 'Senior MATLAB Test Frameworks Engineer (13594BHIL)' 'Experienced' '13593' 'Software Engineer for Design Automation (13593KKOT)' 'Internships' '13577' 'Data Tools Development Intern (13577KCAR)' 'Internships' '13569' 'Senior Compiler Engineer (13569GMAR)' 'Experienced' '13566' 'Simulink Compiler Engineering Intern (13566KCAR)' 'Internships' '13563' 'LRM Updates to Verilog and SystemVerilog Emitters (13563KCAR)' 'Internships' '13558' 'Software Engineer GUI and Hardware Connectivity (13558GMAR)' 'Experienced' '13545' 'C++ Software Engineer  Graphical User Interfaces (13545BWAL)' 'Experienced' '13538' 'OAuth / LTI Integration Software Intern (13538KCAR)' 'Internships' '13532' 'Software Scope Infrastructure Internship (13532KCAR)' 'Temps/Consultants' '13531' 'Software Scope Infrastructure Internship (13531KCAR)' 'Temps/Consultants' '13530' 'FPGA Prototyping Intern (13530KCAR)' 'Internships' '13529' 'HDL Verifier Intern (13529KCAR)' 'Internships' '13518' 'Communications and Signal Processing HDL Development (13518KCAR)' 'Internships' '13516' 'FPGA Prototyping for Image, Video and Computer Vision (13516KCAR)' 'Internships' '13515' 'Audio Intern (13515KCAR)' 'Internships' '13508' 'C++ Software Developer  Discrete Event Simulation (13508GMAR)' 'Experienced' '13473' 'C++ Developer (13473MCAR)' 'Experienced' '13466' 'JavaScript Software Engineer (13466SMAR)' 'Experienced' '13422' 'Release Engineering Intern (13422KCAR)' 'Internships' '13412' 'Web Applications Development Intern (13412KCAR)' 'Internships' '13403' 'Documentation Toolsmith (13403JJUS)' 'Experienced' '13393' 'Robotics Software Developer (13393GMAR)' 'Experienced' '12383' 'C++ Software Engineer  Compiler Engineer for MATLAB (12383MCAR)' 'Experienced' '12382' 'Senior Compiler Engineer LLVM (12382GMAR)' 'Experienced' '12379' 'Senior Software Engineer  Internet Services Platform (12379JJUS)' 'Experienced' '12377' 'Senior Software Engineer  Graphics / OpenGL (12377SMAR)' 'Experienced' '12369' 'Build System Software Engineer (12369BHIL)' 'Experienced' '12365' 'Test Infrastructure Software Developer (12365BHIL)' 'Experienced' '12362' 'Software Engineer  Data Analysis Tools (12362SMAR)' 'Experienced' '12358' 'Web UI Developer  MATLAB Editor (12358MCAR)' 'Experienced' '12356' 'Software Developer  Graphic Export and Printing (12356SMAR)' 'Experienced' '12352' 'JavaScript Build Infrastructure Engineer (12352BHIL)' 'Experienced' '12351' 'Software Engineer  Maven Build and Deployment (12351BHIL)' 'Experienced' '12342' 'Web Application Developer (12342JJUS)' 'Experienced' '12324' 'Compiler Engineer  Parallel Computing (12324GMAR)' 'Experienced' '12322' 'Web UI Software Engineer (12322SMAR)' 'Experienced' '12320' 'Senior Web Application Developer (12320JJUS)' 'Experienced' '12317' 'Senior Software Engineer  Enterprise Solutions (12317MCAR)' 'Experienced' '12316' 'C++ Software Engineer  Simulink Architecture Group (12316GMAR)' 'Experienced' '12315' 'Computer Vision Systems Software Engineer (12315GMAR)' 'Experienced' '12312' 'Web Application Developer  Groovy on Grails/Java (12312BHIL)' 'Experienced' '12308' 'Analog and Mixed Signal Simulation Developer (12308GMAR)' 'Experienced' '12305' 'Compiler Engineer (12305GMAR)' 'Experienced' '12304' 'Software Engineer  Data Streaming and Visualization (12304BWAL)' 'Experienced' '12301' 'C++ Algorithm Optimization Software Engineer (12301GMAR)' 'Experienced' '12300' 'Senior Compiler Developer  Autogenerate C from Simulink (12300GMAR)' 'Experienced' '12284' 'C++ Algorithm Developer  Stateflow (12284GMAR)' 'Experienced' '12275' 'C++ Code Generation Software Engineer (12275GMAR)' 'Experienced' '12240' 'Java Software Engineer Distributed Application Development (12240BHIL)' 'Experienced' '12229' 'Senior Engineer, Requirements Modeling and Management (12229GMAR)' 'Experienced' '12227' 'Senior Engineer, Model Static Analysis (12227GMAR)' 'Experienced' '12217' 'Release Engineering Intern (12217BHIL)' 'Internships' '12174' 'Internal Tools Developer (12174BHIL)' 'Experienced' '12159' 'Release Engineer (12159BHIL)' 'Experienced' '12155' 'Web Software Test Infrastructure and Architecture Engineer (12155BHIL)' 'Experienced' '12072' 'Web GUI Engineer (12072SMAR)' 'Experienced' '12042' 'Web Application Developer  Commerce Systems (12042JJUS)' 'Experienced' '12029' 'Senior Software Engineer  Graphics / OpenGL (12029SMAR)' 'Experienced' '11941' 'Numerical Simulation Software Engineer (11941GMAR)' 'Experienced' '11883' 'Senior Software Engineer, Program Analysis/Formal Methods (11883GMAR)' 'Experienced' '11715' 'Product Release Engineer (11715BHIL)' 'New Graduate' '11669' 'Web UI Software Engineer (11669JJUS)' 'Experienced' '11506' 'Support Package for Mobile devices (11506SMAR)' 'Internships' '11471' 'DevOps Engineer  Cloud/SaaS (11471SMAR)' 'Experienced' '11431' 'Senior Compiler Engineer (11431MCAR)' 'Experienced' '11373' 'Web User Interface  Senior Software Developer (11373GMAR)' 'Experienced' '11302' 'Senior Technical Lead  Statistics / Machine Learning (11302SMAR)' 'Experienced' '11254' 'Senior Software Developer  Graphical Language Editors (11254BWAL)' 'Experienced' '11087' 'Web MATLAB Desktop UI Developer (11087MCAR)' 'Experienced' '10415' 'Senior C++ Compiler Engineer (10415MCAR)' 'Experienced' '10089' 'Senior Software Test Infrastructure Developer (10089BHIL)' 'Experienced' '10038' 'MATLAB Software Engineer  Web GUIs (10038SMAR)' 'Experienced' '9537' 'Senior Compiler Engineer (9537MCAR)' 'Experienced' '8916' 'MATLAB on the Web  Backend Developer (8916JJUS)' 'Experienced' '8904' 'Senior Software Engineer  Web UI Platform (8904JJUS)' 'Experienced' '8892' 'Simulink Compiler Engineer (8892GMAR)' 'Experienced'
Perhaps after searching through this list of jobs, you realize that you would prefer an application engineering job. You can do more searches by placing the commands that are illustrated above into a jobsReader function. The jobsReader function takes a filename as the first input followed optionally by the location and department strings. You can specify 'all' for either location or department to get a complete listing. Save the content to a file and then use your jobsReader function. You can save the RSS feed by using websave to save the contents to a file.
filename = 'rss.xml';
websave(filename, rssURL);
Specify your new job search criteria and obtain the data from the file.
myLocation = 'USMANatick'; myDepartment = 'Application Engineering'; applicationEngineeringJobs = jobsReader(filename,myLocation,myDepartment)
applicationEngineeringJobs = Reference Title Type _________ ____________________________________________________________________ _____________ '14075' 'Application Engineer HardwareintheLoop (14075SMAR)' 'Experienced' '11969' 'Senior Application Engineer (11969SMAR)' 'Experienced' '11943' 'Principal Application Engineer  HardwareintheLoop (11943SMAR)' 'Experienced'
Rather than saving the content to a file, you can specify a function handle as a content reader for webread. When the ContentReader property is specified in weboptions, webread downloads the content to a temporary file and reads the file using the specified function.
Specify an anonymous function handle to return application engineering jobs in Natick.
contentReader = @(filename)jobsReader(filename,myLocation,myDepartment);
options = weboptions('ContentReader',contentReader);
applicationEngineeringJobs = webread(rssURL,options)
applicationEngineeringJobs = Reference Title Type _________ ____________________________________________________________________ _____________ '14075' 'Application Engineer HardwareintheLoop (14075SMAR)' 'Experienced' '11969' 'Senior Application Engineer (11969SMAR)' 'Experienced' '11943' 'Principal Application Engineer  HardwareintheLoop (11943SMAR)' 'Experienced'
Use the Facebook Graph API to obtain more information about MathWorks. The Graph API returns a JSON string. webread parses the string into a structure.
facebookAPI = 'https://graph.facebook.com/'; company = 'MathWorks'; facebookURL = [facebookAPI company]; mathworks = webread(facebookURL)
mathworks = id: '19564567448' about: 'Over one million people around the world speak M...' can_post: 0 category: 'Computers/technology' category_list: [1x1 struct] checkins: 2063 company_overview: 'MathWorks products are used throughout the autom...' cover: [1x1 struct] description: 'MATLAB® is a highlevel language and interactive...' founded: '1984' has_added_app: 0 is_community_page: 0 is_published: 1 likes: 211377 link: 'https://www.facebook.com/MathWorks' location: [1x1 struct] mission: 'Technology Our goal is to change the world by a...' name: 'MathWorks' parking: [1x1 struct] phone: '(508) 6477000' talking_about_count: 2268 username: 'MathWorks' website: 'http://www.mathworks.com' were_here_count: 2063
Display overview information about MathWorks. Use line wrapping to display the text; otherwise it will get trimmed when publishing the blog post. Use the function strsplit to split the string at each white space into a cell array of strings. Compute indices indicating which elements to print. With the indices, use the function strjoin to construct the string to print. Use the function strtrim to ensure that each line begins with a character.
splitData = strsplit(mathworks.company_overview,' '); maxLengthOfLine = 80; lengthOfString = cellfun(@length, splitData) + 1; while ~isempty(splitData) elementsToPrint = cumsum(lengthOfString) <= maxLengthOfLine; if ~any(elementsToPrint) elementsToPrint(1) = true; end lineout = strjoin(splitData(elementsToPrint)); fprintf('%s\n', strtrim(char(lineout))); splitData = splitData(~elementsToPrint); lengthOfString = lengthOfString(~elementsToPrint); end
MathWorks products are used throughout the automotive, aerospace, communications, electronics, and industrial automation industries as fundamental tools for research and development. They are also used for modeling and simulation in increasingly technical fields, such as financial services and computational biology. MathWorks software enables the design and development of a wide range of advanced products, including automotive systems, aerospace flight control and avionics, telecommunications and other electronics equipment, industrial machinery, and medical devices. More than 5000 colleges and universities around the world use MathWorks solutions for teaching and research in a broad range of technical disciplines.
Display the MathWorks cover source image from Facebook. webread returns an MbyNby3 image when reading image/jpeg content. Use imshow to display the image.
RGB = webread(mathworks.cover.source); figure imshow(RGB)
You can use webread to read specific fields from Facebook by specifying a 'fields' query parameter. Use the query parameter syntax to read the 'mission' and 'likes' fields.
fields = 'mission,likes'; data = webread(facebookURL,'fields',fields)
data = mission: 'Technology Our goal is to change the world by accelerating ...' likes: 211377 id: '19564567448'
To understand more about MathWorks, print the mission statement. The function strprint is a local function created from the code in the previous section to split and print the company overview. In this case, use strsplit to split the mission statement into sentences.
sentences = strsplit(data.mission, '.'); for k = 1:length(sentences) fprintf('\n') if ~isempty(sentences{k}) sentence = sprintf('%s.\n', sentences{k}); strprint(sentence) end end
Technology Our goal is to change the world by accelerating the pace of discovery, innovation, development, and learning in engineering and science. We work to provide the ultimate computing environment for technical computation, visualization, design, simulation, and implementation. We use this environment to provide innovative solutions in a wide range of application areas. Business We strive to be the leading worldwide developer and supplier of technical computing software. Our business activities are characterized by quality, innovation, and timeliness; competitive awareness; ethical business practices; and outstanding service to our customers. Human We cultivate an enjoyable, vibrant, participatory, and rational work environment that nurtures individual growth, empowerment, and responsibility; appreciates diversity; encourages initiative and creativity; values teamwork; shares success; and rewards excellence. Social We actively support our communities and promote social and environmental responsibility. Learn more about our Social Mission.
Display the location information to assist in travel planning.
location = mathworks.location
location = city: 'Natick' country: 'United States' latitude: 42.3 longitude: 71.349 state: 'MA' street: '3 Apple Hill Dr' zip: '01760'
You can use a RESTful web service from the Federal Aviation Administration (FAA) to find out the local weather and determine whether any flight delays exist for the Boston Logan International Airport. The web service allows you to specify JSON format for the output.
faaURL = 'http://services.faa.gov/airport/status/'; airportCode = 'BOS'; url = [faaURL airportCode]; format = 'application/json'; options = weboptions('Timeout',10); data = webread(url,'format',format,options)
data = delay: 'false' IATA: 'BOS' state: 'Massachusetts' name: 'General Edward Lawrence Logan International' weather: [1x1 struct] ICAO: 'KBOS' city: 'Boston' status: [1x1 struct]
Display the time.
data.weather.meta
ans = credit: 'NOAA's National Weather Service' updated: '8:54 AM Local' url: 'http://weather.gov/'
Display the weather information, including temperature and wind speed.
data.weather
ans = visibility: 2 weather: 'Light Rain Fog/Mist' meta: [1x1 struct] temp: '43.0 F (6.1 C)' wind: 'Northwest at 10.4mph'
Display the airport status.
data.status
ans = reason: 'No known delays for this airport.' closureBegin: '' endTime: '' minDelay: '' avgDelay: '' maxDelay: '' closureEnd: '' trend: '' type: ''
If you have access to the Mapping Toolbox™, you can display a map showing you a route from Boston Logan International Airport to MathWorks. Use the gpxread function to read a sample GPX file that ships with the Mapping Toolbox and contains a route from Boston Logan to MathWorks. Overlay that route onto an OpenSteetMap base layer using the webmap function. Display the route using the wmline function.
route = gpxread('sample_route'); webmap('openstreetmap') wmline(route.Latitude,route.Longitude)
Print out turnbyturn directions.
turns = route(~cellfun(@isempty, route.Description)); for k = 1:length(turns) fprintf('%s. %s\n',num2str(k),turns(k).Description) end
1. Head southeast 2. Keep left at the fork, follow signs for I90 W/I93 S/Williams Tunnel/Mass Pike and merge onto I90 W Partial toll road 3. Take exit 13 to merge onto MA30 E/Cochituate Rd toward Natick Partial toll road 4. Turn right onto Speen St 5. Merge onto MA9 E/Worcester St via the ramp on the left to Boston 6. Slight right onto Apple Hill Dr 7. Turn left toward Apple Hill Dr 8. Take the 1st right onto Apple Hill Dr Destination will be on the right 9. The MathWorks, Inc., Natick, MA.
Your final destination is the MathWorks campus. You can obtain a highresolution image of the area from the United States Geological Survey (USGS) National Map Server using the Web Map Service (WMS) protocol. The function webread returns an RGB image when reading 'image/jpeg' content type. Specify the WMS query parameters to obtain the image of interest.
usgsURL = 'http://raster.nationalmap.gov/arcgis/services/Orthoimagery/USGS_EROS_Ortho/ImageServer/WMSServer'; service = 'WMS'; layers = '0'; crs = 'CRS:84'; contentType = 'image/jpeg'; height = 512; width = 512; request = 'GetMap'; latlim = [42.298530964,42.301131366]; lonlim = [71.353251832,71.348120766]; format = '%12.9f'; bbox = [ ... sprintf(format,lonlim(1)) ',' sprintf(format,latlim(1)), ',', ... sprintf(format,lonlim(2)) ',' sprintf(format,latlim(2))]; version = '1.3.0'; RGB = webread(usgsURL, ... 'SERVICE', service, 'LAYERS', layers, 'CRS', crs, 'FORMAT', contentType, ... 'HEIGHT', height', 'WIDTH', width, 'REQUEST', request, 'BBOX', bbox);
Display the image.
figure imshow(RGB)
I am hopeful that this post gives you an introduction to accessing RESTful web services using MATLAB. There are many RESTful web services available on the Internet. How might you use these tools in your work? Let us know here.
Best of luck in your job search!
Get
the MATLAB code
Published with MATLAB® R2014b
Idin's pick for this week is Koch curve by Dimitrios Piretzidis.
This week being a holiday week in the US, and start of winter in the Northern Hemisphere where our corner of the world (Natick, Massachusetts) is located, I figured I would pick something a little more fun and winter related, like a snowflake!
The “Koch Snowflake” or “Koch Star” is a wellknown fractal curve (actually one of the earliest fractals to be described). There are many places on the web where you can read about the Koch Snowflake:
The mathematical formula for the Koch snowflake isn't very complicated, and there are at least a halfdozen Kochrelated submissions I found on File Exchange. I picked this particular one out of the bunch because it worked straight out of the box, and it has a simple user interface that makes generating a snowflake as easy as entering an integer and pressing a button! Be warned, though, for any iteration order above 6 or 7, it will take a while to generate your snowflake! A close runnerup for me was this submission:
http://www.mathworks.com/matlabcentral/fileexchange/27577fractalcurves
This gives you more fractal functions, and the code is perhaps a bit more polished, but it's not as easy as a graphical user interface with a push button!
Happy winter, everyone!
As always, your thoughts and comments here are greatly appreciated.
Get
the MATLAB code
Published with MATLAB® R2015a
help greetings
% greetings Seasonal holiday fractal. % greetings(phi) generates a seasonal holiday fractal that depends % upon the parameter phi. The default value of phi is the golden ratio. % What happens for other values of phi? Try both simple fractions % and floating point approximations to irrational values.
which greetings
http://www.mathworks.com/moler/ncm/greetings.m
phi = (1+sqrt(5))/2 greetings(phi)
greetings((1+sqrt(17))/2)
greetings((1+sqrt(13))/5)
greetings(22/7)
greetings(.469)
... and best wishes for the new year.
Get
the MATLAB code
Published with MATLAB® R2014b
]]>
Jiro's pick this week is "Distribute Figures" by Anders Simonsen.
I must be obsessed with figure management, MATLAB figure that is. I've previously Picked three File Exchange entries around this topic: cascade, togglefig, and Useful Figure Management Utilities. No wonder, I was drawn to this entry by Anders.
Have you ever generated lots of figures and wished you could quickly view all of them on one screen?
Anders's distFig lets you do that. Let's say that you have four figure windows.
logo figure surf(peaks) figure plot(rand(10,3)) figure scatter3(rand(50,1),rand(50,1),rand(50,1))
Just call distFig, and it distributes your figures on the screen. Pass in some parameters to control how the figures are distributed.
distFig('Rows',2,'Columns',2)
Anders provides detailed help with the function, as well as some sample code demonstrating the various uses. I also appreciate his responsiveness to user comments. You can see that he quickly updated his entry to fix some issues that came up in R2014b. Great stuff, Anders!
Comments
Give this a try and let us know what you think here or leave a comment for Anders.
Get
the MATLAB code
Published with MATLAB® R2014b
As you can imagine, to debug and validate my work, I need to compare a lot of signals in my model with the original data I received.
As a first step, I am trying to validate the static behavior of the model. This means that I simulate the model for only one time step and validate that the outputs are as expected for different inputs.
If the output does not match, then I need to inspect intermediary signals to figure out where the difference originates.
Display blocks
The first, most obvious, way to see the value of a signal in Simulink is the Display block. In my case, using Display blocks quickly made my model look like:
As you can see, the model quickly becomes a mess. And since those displays are only temporary, I will need to remove them later.
Value Label Display
Another option that I often find more convenient to quickly view the value of a signal for debugging is the Value Label Display.
In the Display menu, I can decide if I want the label to be displayed when hovering over the block, or when clicking on it.
There are also options to control how often the labels are updated, and their accuracy. In my case, I chose the long format to see all the digits available.
Once this is enabled, you can simply click on blocks to enable the displays and simulate the model to see the values.
Once your debugging is completed, you can remove all the labels in one click from the Display menu.
Now it's your turn
What is your preferred method to to view signal values during debugging? Let us know by leaving a comment here.
]]>I'm pleased to have Dave Foti back for a discussion of subclassing and class hierarchies in MATLAB. Dave manages the group responsible for objectoriented programming features in MATLAB.
In computer science, a class represents a set of objects that share a common definition, usually including data and functions or behaviors. A subclass represents a subset of objects that share a more specific definition usually by adding or specializing data and/or functions defined in a superclass. In practice subclassing has been used to reuse functionality for a new and more specialized purpose. In recent years, this usage has become largely discouraged due to the tendency for introducing errors in programs that reuse superclasses (often called base classes) for applications not anticipated by the original class. A class author hides the inner workings of a class to make the class useful as an abstract representation. That same hiding of details can also hinder the ability of a programmer to anticipate how the class will function in a new role. Let's explore some cases where subclassing does make sense in MATLAB and how to avoid some common programming pitfalls.
One of the most obvious uses for subclassing in MATLAB is a situation requiring a model of a real physical system that has already been classified. In such a situation, it might make sense to build a software model that closely resembles the real physical classification system. For example, one could create objects based on the Standard Model of particle physics [1]. In this way, the classification system already exists and can simply be represented in software. There is still design work to be done in how best to represent a physical system in software for a specific purpose, but all the classes can be built at roughly the same time by the same person or team and thus avoid the pitfalls of reusing an existing class later on after its details have been forgotten. Here is what a superclass might look like for elementary particles in the Standard Model:
type ElemParticle
classdef ElemParticle % Intrinsic properties of the particle properties(SetAccess = immutable) Spin Charge Mass end % State variables for a particular particle properties Position = [0 0 0]; Velocity = [0 0 0]; end methods function p = ElemParticle(spin, charge, mass) p.Spin = spin; p.Charge = charge; p.Mass = mass; end end end
We can then create subclasses for fermions and bosons:
type Fermion type Boson
classdef Fermion < ElemParticle methods function p = Fermion(charge, mass) p@ElemParticle(.5, charge, mass); end end end classdef Boson < ElemParticle methods function p = Boson(mass) p@ElemParticle(0, 0, mass); end end end
We can divide the Fermions into Leptons and Quarks:
type Lepton type Quark
classdef Lepton < Fermion methods function p = Lepton(charge, mass) p@Fermion(charge, mass); end end end classdef Quark < Fermion properties Color end methods function p = Quark(charge, mass, color) p@Fermion(charge, mass); p.Color = color; end end end
In the above example, the superclasses and subclasses were designed as a single collection of classes so there is no need for superclasses to anticipate the needs of unknown future subclasses. Any future evolution of the superclasses needs to consider the whole system, but it is a bounded list of classes.
Even if one is not modeling a physical system with an existing classification model that makes sense in software, there are still other examples where a system of classes can be designed at roughly the same time to fulfill some purpose such as managing a variety of data formats or anaysis techniques. For example, we might have a need to visualize some data as points. We might know that we have two kinds of data sets: one that is twodimensional and one that is threedimensional. We design a class that handles 2D data and a subclass that adds support for a third dimension. Here is what the 2D class looks like:
type PointList
classdef PointList properties X Y end methods(Sealed) function show(pc) figure; data = getData(pc); opts = {'+', 'MarkerSize', 3}; if numel(data) == 3 plot3(data{:},opts{:}); else plot(data{:},opts{:}); end end end methods(Access = protected) function data = getData(pc) data = {pc.X, pc.Y}; end end end
Because we know that we want to support 3D data, we design our plotting routine to be able to handle 2D or 3D data. This is an advantage of knowing about all the subclasses in the beginning. Then when we write the 3D subclass:
type PointList3D
classdef PointList3D < PointList properties Z end methods(Access = protected) function data = getData(pc) data = {pc.X, pc.Y, pc.Z}; end end end
it is very easy to support 3D plots without rewriting the whole show method. Here we produce a 3D point cloud for the 3rd harmonic of the Lshaped membrane:
L = membrane(3); [X, Y] = meshgrid(1:length(L), 1:length(L)); pts = PointList3D; pts.Z = L(:); pts.X = X(:); pts.Y = Y(:); show(pts)
Another good reason to use subclassing is when a superclass is designed specifically for reuse. A common design pattern [2] is called the mixin pattern [3]. In this pattern a class is often used to implement one specialized piece of functionality. This class represents all objects that have this common functionality even though such classes may do other things very differently. A mixin class should expect nothing of a subclass except that it uses the mixin functionality. A common example of a mixin class in MATLAB is the class called matlab.mixin.CustomDisplay which is used to customize the MATLAB default object display while retaining appropriate consistency with other objects. Here is an example of this mixin being used on our ElemParticle superclass. To separate this version from the previous example, I've put the new version of the classes in a package called particle.
type +particle/ElemParticle
classdef ElemParticle < matlab.mixin.CustomDisplay % Intrinsic properties of the particle properties(SetAccess = immutable) Spin Charge Mass end % State variables for a particular particle properties Position = [0 0 0]; Velocity = [0 0 0]; end methods function p = ElemParticle(spin, charge, mass) p.Spin = spin; p.Charge = charge; p.Mass = mass; end end methods(Access = protected, Sealed) function propgroups = getPropertyGroups(p) propgroups = matlab.mixin.util.PropertyGroup; propgroups(1).Title = 'Intrinsic'; propgroups(1).PropertyList = getIntrinsicPropertyNames(p); propgroups(2).Title = 'State'; propgroups(2).PropertyList = {'Position','Velocity'}; end end methods(Access = protected) function pnames = getIntrinsicPropertyNames(~) pnames = {'Spin','Charge','Mass'}; end end end
Note that whenever functionality is designed for a superclass such as this, it is important to implement the functionality for the whole class of objects including all instances of subclasses. Thus, the display should not just show what is common to ElemParticle, but should provide for specialization in subclasses. In this case, we add one method that can be specialized by subclasses to provide lists of intrinsic properties. Since we know about all the subclasses, we don't need to worry about what kinds of property groups a future subclass might need. We can see why it is important for superclasses to be designed with the knowledge that subclasses will be created and why it is very helpful if the subclasses are already known during this design. In this case, we know the family of subclasses of ElemParticle and thus we know that they will sometimes need to add new properties, but not whole new groups of properties. We choose the right amount of generality in the superclass so that each subclass can specialize display with minimal effort. In fact Quark is the only particle that needs additional code as shown here:
type +particle/Quark
classdef Quark < particle.Fermion properties Color end methods function p = Quark(charge, mass, color) p@particle.Fermion(charge, mass); p.Color = color; end end methods(Access = protected) function pnames = getIntrinsicPropertyNames(p) pnames = [getIntrinsicPropertyNames@particle.Fermion(p), 'Color']; end end end
We can now create and display an electron and a top quark:
e = particle.Lepton(1, 5.11e4)
e = Lepton with properties: Intrinsic Spin: 0.5 Charge: 1 Mass: 0.000511 State Position: [0 0 0] Velocity: [0 0 0]
t = particle.Quark(2/3, 173.07, 'red')
t = Quark with properties: Intrinsic Spin: 0.5 Charge: 0.66667 Mass: 173.07 Color: 'red' State Position: [0 0 0] Velocity: [0 0 0]
Other common patterns where classes are designed for unbounded subclassing include interfaces [4] and frameworks [5]. Interface classes define how two pieces of code can interact with each other without actually defining how those interactions are accomplished. Because MATLAB is a dynamically typed language, it is possible to write functions to an interface that isn't actually defined as a class. For example many functions in MATLAB can work with any object that implements the size, subsref, and subsasgn functions to perform reshaping operations on arrays. In some cases, it may be useful to define interfaces explicitly so that it is easier to determine what kinds of objects will work with a particular function and so that it is easier to see the capabilities of a class in the explicit interfaces it supports.
When designing a class for extension in the future, it is a good idea to think about how subclasses might restrict evolution of the base class. Some questions to consider are:
When designing methods in the base class, it is important to consider how subclasses will reuse and as necessary adapt these methods while maintaining the principle that subclass instances can be substituted for base class instances [6]. Some questions to ask are:
Once a class has been presented as a base class, its fundamental definition should not change or else subclasses could become invalid or produce erroneous results. The new or changed behaviors and capabilities may not be appropriate for all subclasses. The subclass authors could not possibly anticipate changes to the superclass when the subclasses were created. This is one of the reasons for adhering to the open/closed design principle [7].
Given the importance of considering how a class will be extended by subclasses, MATLAB provides ways to signal that a class wasn't designed for subclassing. Classes with the Sealed attribute set can't be used as superclasses. If an entire class hierarchy was designed together but not intended for extension with additional subclasses, then it can be sealed using the Sealed attribute on the leaf classes and the AllowedSubclasses attribute on the superclasses. MathWorks will at times seal classes or class hierarchies of widely used classes to indicate that these classes have not been designed to support subclasses.
These are some thoughts on subclassing in MATLAB, but I'd like to hear your thoughts as well. Have you used any of the mixin MATLAB classes like matlab.mixin.CustomDisplay or matlab.mixin.Heterogeneous? Do you have suggestions for mixin classes that you would like to see in MATLAB? You can post ideas and comments here.
Get
the MATLAB code
Published with MATLAB® R2014b
An incredible book, published in several editions from 1909 to 1933, by German mathematicians Eugene Jahnke and Fritz Emde, contains definitions and formulas for mathematical functions, handcalculated tables of function values, and meticulous handdrawn 2 and 3dimensional graphs. An English edition was published by Dover in 1945.
I will include only two of their illustrations in this blog. Here is one of their graphs of the gamma function in the complex plane. The surface is the magnitude of the function, with its poles at the integers on the negative real axis. The contour lines are the modulus and phase. I can't imagine how they made this plot in 1933.
Here is my attempt to have MATLAB recreate the Jahnke and Emde plot. It is only partially successful. I do have little ellipses to cut off the spikes  that was tricky. But I haven't begun to label the contours, or shade the edge faces.
I love this plot. It is the complexvalued Bessel function of order zero, also known as the Hankel function. There is a pole, a zero, and a saddle point. Again, imagine doing this by hand. Or, having some grad student do it.
Here we can use color, and our new parula color map, in place of 3d. But we can't get the nice feel of the hand drawn lettering.
[X,Y] = meshgrid(4:1/128:2,1.5:1/128:1.5); H = besselh(0,X+i*Y); s = 3.2; contourf(X,Y,abs(H),0:0.2:s) hold on c2 = s*(angle(H)/(2*pi)+0.5); r2 = s*((180:10:180)/360+0.5); contour(X,Y,c2,r2,'k') hold off axis equal axis tight grid on
Eugene Jahnke and Fritz Emde, Tables of Functions with Formulas and Curves, (4th ed.), Dover, 1945
Get
the MATLAB code
Published with MATLAB® R2014b
Sean's pick this week is Growth  Yet Another Obfuscated MATLAB Code by Hikaru.
This is awesome.
dbtype growth
1 figure(1) 2 clf,hold on,grid on,view(58,38),feval(@(mm,x)mm(@(O)mm(@(Y,U,V,W,v,nN,NN,Mm,mM,md,Nn)... 3 mm(@(G)mm(Y(G),x,0,0,v(mM(x))),mm(@(L)(@(M)(@(n,p,q,r)mm(mm(@(l,n)l{n},{@()1,@(){mm(Y... 4 (L), [25*... 5 v(1. /r*(... 6 n1) ) 25/ 2*(r... 7 1); md((n1),r)*5 ^2(... 8 r1* 1)*5^ 2/ 2;0] ,[0;... 9 0;1] ,[ 1;0 ;0 ], 6) ,M(n... 10 1,0 ,0,v(mM(x)))} }, ~(n>... 11 .0)+ 2*~(n <= 0))) )) ,(@(... 12 d)(( @( u,n ,t ,N )( (( mm (mm ((@(... 13 l,n) l{n}),{(@()mm (@ ()(mm(@()Nn(u (1),... 14 u(2) ,u(3) ,( '*g' ), U,5)) )) ),(@ ()mm... 15 ((@( w, b, p, m, mN ){( mm (@ (A )Nn(... 16 A(1, :) ,A(2,:),A(3,: ), '',... 17 W,N^ 2* 0.3,V ,( mN(N /6 ))),... 18 [u,( u+ n* w)] )) ,( mm (@ (f ,T, S){(... 19 d(u+ n*w,S,T,N1)) ,( d(u+n*w,f(S,n ,p),... 20 f(T, n,p), N 1)) ,( d(u+n *w ,f(S ,n,p... 21 *2), (f (T ,n ,( p)* 2)), N1 )) }, (@(r... 22 ,n,p )r *(Mm(p))nN(r ,n )*NN... 23 (p)+ (1 Mm(p )) *n*( n' *r))... 24 ,NN( m) *n t* Mm (m ), n* Mm(m... 25 )+NN (m)*t))}),9.8 /( 7.14... 26 N), 0.1,2 +1 /10.6 , .102... 27 +1,( @( r) [1 *1 *.32... 28 ,.22 ,1 *.08... 29 ]*r+ [0 ,1,(... 30 .5)] *( 1r)... 31 )))} ,1+(... 32 N>0))))))))))),@(f)mm(@(g)@(m,n,o,p)mm(f(g(g)),m,n,o,p),@(g)@(m,n,o,p)mm(f(g(g)),m,n,... 33 o,p)),O(0),O(1),O(2),O(3),O(4),O(5),O(6),O(7),O(8),O(9)),(@(I)mm(@(l,m)mm(l{1+(I>2)},... 34 char(m((m(58+I)46):(m(59+I)47))+18)),{@(J)J,@str2func},['[O`YS`aWhS1]Z]`ZW\SEWRbV',... 35 'TZ]]`Q`]aaaW\Q]aa_`b[]R^Z]b!TSdOZ/9>GLQTW[^c']))),@feval,9); %% (c) woodrush (2014) ...
growth
They look the same!
Hikaru, also provides the technical details about the programming features used. However, the detail I care about most was left out: How long did this take to make?
It is also a great example of pushing the MATLAB interepreter to its limits.
Maybe some day we can have an obfuscated MATLAB Code Contest, similar to the C one.
Have a great New Year everyone!
Give it a try and let us know what you think here or leave a comment for Hikaru.
Get
the MATLAB code
Published with MATLAB® R2015a
A block unlike any other
Simulink Functions are unlike any other Simulink block, even other userdefined functions like FunctionCall Subsystems. For starters, the block does not use traditional ports. Instead, like MATLAB functions, you specify arguments.
Inside the Simulink Function, arguments are represented by uniquelooking Argument port blocks.
Since they aren’t traditional ports, you won’t find these blocks in the Library Browser. To add or remove arguments, you type right on the block’s mask:
To use a Simulink Function in your model, you employ a Function Caller block.
Simulink associates each instance of this block with a Function by the function name; it requires that every function name in the model hierarchy be unique.
Functions in referenced models
Speaking of hierarchy, one very practical use of Simulink Functions are in a referenced model context. Functions defined in a referenced model are in the scope of its parent. But, the referenced model has to follow some very special design rules, the same as they do when they contain FunctionCall Subsystems.
The simplest way to adhere to these rules is to create a model that contains nothing but your collection of Simulink Function blocks, like so:
Now add a Model block to your parent model, point it to this collection, and you can call any of these functions from the parent, or from any other model further up the hierarchy.
Generating code
Code generated from Simulink Functions is as modular and straightforward as you would expect. Each function is defined in its own source file. The function prototype always takes the form of:
void foo(type_T in1, type_T in2, type_T *out1, type_T *out2)
Simulink is smart enough to recognize when you are trying to emulate passbyreference, and you name your output the same as your input (like I did in the timestwo example in the images above). When you do, the generated function prototype looks just like this:
void foo(type_T *x)
Now it’s your turn
This opens up a whole new way to model. What kind of design patterns can you create?
]]>In the summer of 2013 we were closing in a choice for the new MATLAB colormap. We were down to tweaking and finetuning.
But ... we needed a name!
For my many experiments, I had fallen back on an old graduate school habit of naming things after characters from Lord of the Rings. So I had filenames such as gandalf_20130623a.m and faramir_20130712b.m. (Faramir is my favorite LOTR minor character.) I certainly knew that wasn't going to work for the final name. How to choose one?
I looked over the names of the existing colormaps:
jet  
hsv  
hot  
cool  
spring  
summer  
autumn  
winter  
gray  
bone  
copper  
pink  
lines  
colorcube  
prism  
flag 
Well, we've got a few identifiable themes in there:
And there are some oddball names, such as jet, lines, prism, and flag. None of that seemed inspirational for naming a new colormap.
I decided to look for something descriptive. But descriptive of what?
To remind you, here's what the new colormap looks like:
showColormap(parula,'bar')
I picked the main colors (to my eye, these are blue, green, orange, and yellow) in the new colormap and started doing searches using these color names and different kinds of objects. Animals seemed obvious. I actually started with fish, but that got nowhere fast.
Then I tried birds, and up popped the tropical parula:
As Wikipedia describes it, the tropical parula "has mainly bluegrey upperparts, with a greenish back patch and two white wingbars. The underparts are yellow, becoming orange on the breast."
Perfect! I didn't really think the team would go for it, though. I sent around an email with the name parula and a picture of the bird to a small group of people working on the visual appearance changes for the new MATLAB graphics system. Somewhat to my surprise, everyone said they liked it. Later, a larger group of senior MATLAB designers reviewed it, and they also liked it. So the name stuck.
That left us with one big problem, however. How is parula pronounced?
I don't know. In trying to find a definitive answer, I have only managed to confuse myself. As a result, I don't even always pronounce it the same way myself. In looking at various references for pronunciation of bird names, I have seen all of these variations:
I know two amateur birders at MathWorkers who told me definitively how to pronounce it. Of course, they each gave me a different answer.
I guess that most American English speakers would choose the first variation, which has the accent on the second syllable. That has certainly been the case at MathWorks headquarters in Massachusetts.
So I'll tell you what I tell MathWorkers here: you can pronounce it however you like!
Get
the MATLAB code
Published with MATLAB® R2014b
Complementary Blocks
In Simulink, there are a few blocks that always come in pairs, for example Goto and From, and Data Store Write and Data Store Read.
In R2014b, if you have one of these blocks, you can simply hover the mouse on top of the block, and a light blue indicator will appear. Click on it and drag to make the complementary block appear.
Now it's your turn
Were you already aware of this feature? Do you find it useful? Let us know by leaving a comment here.
]]>
]]>
Sean's pick this week is the suite of naturalorder sorting tools by Stephen Cobeldick
If you work with data, there are many different naming schemes for files that you will likely encounter. For some of these the numerical sorting order might be the same as the ASCII order sorting order. But it is not always, and how to deal with this is a common question on MATLAB Answers.
I typically try to set it up the following way when naming my own files using the "%0xi" format in num2str. For example:
for ii = [1 2 17 495 3920] % %04i  up to four zeros in front of the number. filenameii = ['file' num2str(ii,'%04i') '.csv']; disp(filenameii) end
file0001.csv file0002.csv file0017.csv file0495.csv file3920.csv
Stephen's naturalorder sorting tools help sort files or names that do not necessarily have this setup. For example, let's look at the ASCII order of a few files using sortrows.
files = {'file1.csv','file111.csv','file21.csv','file211.csv'}.'; disp(sortrows(files))
'file1.csv' 'file111.csv' 'file21.csv' 'file211.csv'
Numerically file111.csv should not be before file21.csv. Now let's use natsort to do this for us:
disp(natsort(files))
'file1.csv' 'file21.csv' 'file111.csv' 'file211.csv'
The other utilies that Stephen has provided allow more control over this and for the extension of it to not just working on single elements of a cell but on a whole full file path. For example, how should I sort the following?
files = {'C:\Documents\Exp1\test1.csv','C:\Documents\Exp2\test1.csv','C:\Documents\Exp2\test2.csv','C:\Documents\Exp1\test2.csv'}.'; disp(files)
'C:\Documents\Exp1\test1.csv' 'C:\Documents\Exp2\test1.csv' 'C:\Documents\Exp2\test2.csv' 'C:\Documents\Exp1\test2.csv'
Sorted by experiment:
disp(natsortfiles(files))
'C:\Documents\Exp1\test1.csv' 'C:\Documents\Exp1\test2.csv' 'C:\Documents\Exp2\test1.csv' 'C:\Documents\Exp2\test2.csv'
To sort by test, we can split the file path into pieces and then use natsortrows on the pieces:
% Split on file separators filepieces = regexp(files, ['' filesep ''], 'split'); filepieces = vertcat(filepieces{:}); disp(filepieces)
'C:' 'Documents' 'Exp1' 'test1.csv' 'C:' 'Documents' 'Exp2' 'test1.csv' 'C:' 'Documents' 'Exp2' 'test2.csv' 'C:' 'Documents' 'Exp1' 'test2.csv'
% Sort them by fourth column (test) then third column (experiment)
[~, idx] = natsortrows(filepieces,[4 3]);
disp(files(idx))
'C:\Documents\Exp1\test1.csv' 'C:\Documents\Exp2\test1.csv' 'C:\Documents\Exp1\test2.csv' 'C:\Documents\Exp2\test2.csv'
These files provide excellent help and are well documented.
My only suggestion for Stephen would be to provide these files together in one File Exchange entry (or as a fourth entry). This is solely for the reason that I am lazy and downloading all of the separate zip files and unpacking them took an extra few minutes. However, since the second two files depend on natsort, they wouldn't work on their own without this process.
If given the choice, how do you choose to store your files? What challenges do you face when you receive files from others or from hardware? Let us know below.
Give it a try and let us know what you think here or leave a comment for Stephen.
Get
the MATLAB code
Published with MATLAB® R2014b
Today I’d like to introduce guest blogger Sarah Wait Zaranek who works for the MATLAB Marketing team here at MathWorks. Sarah previously has written about a variety of topics. Mostly recently, she cowrote a post with me about the new webcam capabilities in MATLAB. Today, Sarah will be discussing datastore, one of the new big data capabilities introduced in MATLAB R2014b.
datastore is used for reading data that is too large to fit in memory. For this example, we will be reading in data from the vehicle census of Massachusetts. It is a catalog of information about vehicle registered from 2008 to 2011. The dataset contains information about individual cars registered including vehicle type, location where the vehicle is housed, rated MPG, and measured CO_2 emissions. You can learn more about the data and even download it yourself, here. I have renamed the files in the demo for clarity's sake, but this is where they came from originally.
As mentioned, a datastore is an object useful for reading collections of data that are too large to fit in memory.
datastore can work with a single file or a collection of files. In this case, we will be reading from a single file. Our file does not include variable names at the top of the file. They are listed in separate header file as defined below.
% Define Data File and Header File datafile = 'vehiclebig.csv'; headerfile = 'varnames.txt'; % Read in Variable Names fileID = fopen(headerfile); varnames = textscan(fileID,'%s'); varnames = varnames{:}; fclose(fileID);
We can now create our datastore by giving the name of the data file as the import to the datastore function. We also specify our datastore not use the first row of our file as variable names. We will set those variable names explicitly using the names found in the 'varnames.txt' file.
ds = datastore(datafile,'ReadVariableNames',false); % Set Variable Names ds.VariableNames = varnames
ds = TabularTextDatastore with properties: Files: { 'H:\Documents\LOREN\MyJob\Art of MATLAB\SarahZ\datastore\vehiclebig.csv' } ReadVariableNames: false VariableNames: {'record_id', 'vin_id', 'plate_id' ... and 42 more} Text Format Properties: NumHeaderLines: 0 Delimiter: ',' RowDelimiter: '\r\n' TreatAsMissing: '' MissingValue: NaN Advanced Text Format Properties: TextscanFormats: {'%f', '%f', '%f' ... and 42 more} ExponentCharacters: 'eEdD' CommentStyle: '' Whitespace: ' \b\t' MultipleDelimitersAsOne: false Properties that control the table returned by preview, read, readall: SelectedVariableNames: {'record_id', 'vin_id', 'plate_id' ... and 42 more} SelectedFormats: {'%f', '%f', '%f' ... and 42 more} RowsPerRead: 20000
Note that we haven't read in our data yet. We have just provided an easy way to access it through ds, our datastore.
A really nice thing about a datastore is that you can preview your data without having to load it all into memory. datastore reads the data into a table which is a data type in MATLAB designed to work well with tabular data.
data = preview(ds); whos data data(:,1:7) % Look at first 7 variables
Name Size Bytes Class Attributes data 8x45 21426 table ans = record_id vin_id plate_id me_id owner_type start_odom start_date _________ ______ __________ __________ __________ __________ ____________ 2 1 5.3466e+06 1.2801e+07 1 NaN '20111107' 4 1 5.3466e+06 1.1499e+07 1 NaN '20090827' 7 2 6.6148e+06 1.2801e+07 1 NaN '20111119' 9 2 6.6148e+06 1.1499e+07 1 NaN '20080701' 10 3 6.4173e+06 1.2801e+07 1 NaN '20111206' 1 1 5.3466e+06 1 1 30490 '20090901' 3 1 5.3466e+06 2 1 55155 '20101002' 5 2 6.6148e+06 3 1 5 '20080702'
By default, datastore will read in every column of our dataset. datastore makes an educated guess for the appropriate format for each column (variable) of our data. We can, however, specify a subset of columns or different formats if we wish.
We can specify which variables (columns) by setting the SelectedVariableNames property of our datastore. In this case, we only want to bring in 5 columns out of the 45.
ds.SelectedVariableNames = {'model_year', 'veh_type', ... 'curbwt','mpg_adj','hybrid'}; preview(ds)
ans = model_year veh_type curbwt mpg_adj hybrid __________ ________ ______ _______ ______ 2008 'Car' 3500 21.65 0 2008 'Car' 3500 22.54 0 2008 'SUV' 4500 16 0 2008 'SUV' 4500 17 0 2005 'Truck' 5000 13.29 0 2008 'Car' 3500 22.09 0 2008 'Car' 3500 21.65 0 2008 'SUV' 4500 16.66 0
We can adjust the format of the data we wish to access by using the SelectedFormats property. We can specify to bring in the vehicle type as a categorical variable by using the %C specifier. You can learn more here about the benefits of using categorical arrays.
ds.SelectedFormats; ds.SelectedFormats{2} = '%C' % read in as a categorical
ds = TabularTextDatastore with properties: Files: { 'H:\Documents\LOREN\MyJob\Art of MATLAB\SarahZ\datastore\vehiclebig.csv' } ReadVariableNames: false VariableNames: {'record_id', 'vin_id', 'plate_id' ... and 42 more} Text Format Properties: NumHeaderLines: 0 Delimiter: ',' RowDelimiter: '\r\n' TreatAsMissing: '' MissingValue: NaN Advanced Text Format Properties: TextscanFormats: {'%*q', '%*q', '%*q' ... and 42 more} ExponentCharacters: 'eEdD' CommentStyle: '' Whitespace: ' \b\t' MultipleDelimitersAsOne: false Properties that control the table returned by preview, read, readall: SelectedVariableNames: {'model_year', 'veh_type', 'curbwt' ... and 2 more} SelectedFormats: {'%f', '%C', '%f' ... and 2 more} RowsPerRead: 20000
We can use the read function to read in a chunk of our data. By default, read reads in 20000 rows at a time. This value can be adjusted using the RowsPerRead property.
testdata = read(ds);
whos testdata
Name Size Bytes Class Attributes testdata 20000x5 683152 table
After you read in a chunk, you can use the hasdata function to see if there is still additional data available to read from the datastore.
hasdata(ds)
ans = 1
By using hasdata and read in a while loop with your datastore, you can read in your entire dataset a piece at a time. We will put in a counter just to track how many read operations took place in our loop.
counter = 0; while hasdata(ds) % Read in Chunk dataChunk = read(ds); counter = counter + 1; end counter
counter = 825
By using reset, we can reset our datastore and start reading at the beginning of the file.
reset(ds)
Now that we see how to get started using datastore, let's look at 3 different ways to use it to work with a dataset that does not entirely fit in the memory of your machine.
If you are interested in only processing certain columns of your text file and those columns can fit in memory, you can use datastore to bring in those particular columns from your text file. Then, you can work with that data directly in memory. In this example, we are only interested in the model year and vehicle type of the cars that were registered. We can use readall instead of read to import all the selected data instead of just a chunk of it at a time.
ds.SelectedVariableNames = {'model_year', 'veh_type'}; cardata = readall(ds); whos cardata
Name Size Bytes Class Attributes cardata 16145383x2 161456272 table
Now that you have the data read into MATLAB, you can work with it like you would normally work with your data in MATLAB. For this example, we will just use the new histogram function introduced in R2014b to look at the distribution of vehicle model years registered.
figure histogram(cardata.model_year) hold on histogram(cardata.model_year(cardata.veh_type == 'Car')) hold off xlabel('Model Year') legend({'all vehicles', 'only cars'},'Location','southwest')
Another way to subset your data is to filter the data down a chunk at time. Using datastore you can read in a chunk of data and keep only data you need from that chunk. You then continue this process, chunk by chunk, until you reach the end of the file and have only the subset of the data you want to use.
In this case, we want to extract a subset of data for cars that were registered in 2011. The new variables we are loading in (e.g., q1_2011), contain either a one or zero. Ones represent valid car registrations during that time. So we only save the rows which contain a valid registration sometime in 2011 and discard the rest.
reset(ds) ds.SelectedVariableNames = {'model_year','veh_type',... 'q1_2011','q2_2011','q3_2011','q4_2011'}; data2011 = table; while hasdata(ds) % Read in Chunk dataChunk = read(ds); % Find if Valid During Any Quarter in 2011 reg2011 = sum(dataChunk{:,3:end},2); % Extract Data to Keep (Cars Registered in 2011) idx = reg2011 >= 1 & dataChunk.veh_type == 'Car'; % Save to Final Table data2011 = [data2011; dataChunk(idx,1:2)]; end whos data2011 figure histogram(data2011.model_year) xlabel('Model Year') legend({'cars registered in 2011'})
Name Size Bytes Class Attributes data2011 3503434x2 35036782 table
But what if we can't hold the subset of the data we are interested in analyzing in memory? We could instead process the data a section at a time and then combine intermediate results to get a final result. In this case, let's look at the % of hybrid cars registered every quarter (in terms of total cars registered). So, we compute a running total of the number of cars registered per quarter as well as the number of hybrids registered per quarter. Then, we calculate the final % when we have read through the entire dataset.
% Reset Datastore reset(ds) % Select Data to Import quarterNames = varnames(end15:end); ds.SelectedVariableNames = [{'veh_type', 'hybrid'} quarterNames']; % Read in Vehicle Type as a Categorical Variable ds.SelectedFormats{1} = '%C'; totalCars = zeros(length(quarterNames),1); totalHybrids = zeros(length(quarterNames),1); while hasdata(ds) % Read in Chunk dataChunk = read(ds); for ii = 1:length(quarterNames) % Loop over car model years % Extract Data idx = dataChunk{:,quarterNames(ii)}== 1 & dataChunk.veh_type == 'Car'; idxHy = idx & dataChunk.hybrid == 1; % Perform Calculation totalCarsChunk = sum(idx); totalHybridsChunk = sum(idxHy); % Save Result totalCars(ii) = totalCarsChunk + totalCars(ii); totalHybrids(ii) = totalHybridsChunk + totalHybrids(ii); end end
percentHybrid = (totalHybrids./totalCars)*100; figure scatter(1:length(percentHybrid),percentHybrid,'filled') xlabel('Inspection Year') ylabel('% Hybrids') % Label tick axes ax = gca; ax.TickLabelInterpreter = 'none'; ax.XTick = 1:length(quarterNames); ax.XTickLabel = quarterNames; ax.XTickLabelRotation = 45;
You can also use datastore as the first step to creating and running your own MapReduce algorithms in MATLAB. Learn more about running MapReduce algorithms with MATLAB here. Perhaps this topic will be a blog post in the future.
Do you think you can use datastore with your big data? Let us know here.
Get
the MATLAB code
Published with MATLAB® R2014b
Our plots of the first eigenfunction of the Lshaped membrane have changed several times over the last fifty years.
Ph.D. thesis. Two dimensional contour plot. Calcomp plotter with ballpoint pen drawing on paper mounted on a rotating drum.
Calcomp plotter. Math department individual study student project. Original 3d perspective and hidden line algorithms. Unfortunately, I don't remember the student's name.
The first MathWorks documentation and the first surf plot. Apple Macintosh daisy wheel electric typewriter with perhaps 72 dotsperinch resolution.
Laser printer with much better resolution.
Sun workstation with CRT display and color. Hot colormap, proportional to height.
A hand made lighting model. I have to admit this looks pretty phony today.
Good lighting and shading. But at the time this was just a nice graphics demonstration. It was a few more years before it became the official company logo.
A prototype 3d printer. Pure granular silicon, a few dozen grains per inch. (Thanks to Robin Nelson and David Eells for curating photos.)
Cleve's Corner, "The Growth of MATLAB and The MathWorks over Two Decades", MathWorks News&Notes, Jan. 2006.
Get
the MATLAB code
Published with MATLAB® R2014b
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 handson MATLAB experience via the use of an integrated, webbased 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.
]]>Jiro's pick this week is roundsd by François Beauducel.
There are a number of entries, including this, this, this, and this, that perform this useful task of rounding numbers to specific significant digits. François provides a simple, wellimplemented function with good documentation. Here's how it works.
format long
roundsd(pi,3)
ans = 3.140000000000000
roundsd(pi/100,3)
ans = 0.031400000000000
It works with arrays as well.
roundsd([pi pi/100],3)
ans = 3.140000000000000 0.031400000000000
Did you know that in R2014b, the round function gives you the option to do the same thing? It provides the option of rounding to significant digits or to digits respect to the decimal point.
round([pi pi/100],3)
ans = 3.142000000000000 0.031000000000000
round([pi pi/100],3,'significant') % Reset format format short
ans = 3.140000000000000 0.031400000000000
Thanks, François (and others), for providing this functionality to folks using MATLAB versions prior to R2014b!
Comments
Give it a try and let us know what you think here or leave a comment for François.
Get
the MATLAB code
Published with MATLAB® R2014b
An email thread about a daterelated design question had me thinking about terminology this week. (Note: There's nothing about image processing in today's post. Sorry!)
Specifically, I was thinking about the term time zone. In the U.S. you might see a date written this way:
November 21, 2004 1:02 PM (EST)
EST means Eastern Standard Time in the U.S. Or you might see a date written like this:
11/21/2004 8:30 (04:00)
Many people refer to EST or 04:00 as the "time zone." However, because I work across the hall from the chief designer and implementer (Hi, P^2) of the new datetime suite of functions in MATLAB R2014b, I know that these are more properly referred to as time zone offsets.
A time zone is a much more complex concept than a simple offset from UTC. For example, a time zone encompasses both current and historical rules for using Daylight Saving Time.
The time zone offset for a given time depends not only on the time zone but also both the year and the dayofyear.
The dependence on the dayofyear is pretty obvious. In the fall season where I live, the time zone offset depends on whether the day is before or after the first Sunday in November.
d1 = datetime('October 15, 2014','TimeZone','America/New_York')
d1 = 15Oct2014
tzoffset(d1)
ans = 4:00
d2 = datetime('November 15, 2014','TimeZone','America/New_York')
d2 = 15Nov2014
tzoffset(d2)
ans = 5:00
Most of the state of Arizona does not observe Daylight Saving Time.
d3 = datetime('October 15, 2014','TimeZone','America/Phoenix')
d3 = 15Oct2014
tzoffset(d3)
ans = 7:00
d4 = datetime('November 15, 2014','TimeZone','America/Phoenix')
d4 = 15Nov2014
tzoffset(d4)
ans = 7:00
The dependence of time zone offset on the year is somewhat less obvious. It turns out that Daylight Saving Time rules change occasionally. In the U.S. the rules were last changed in 2007. Let's compare the time zone offset for noon on Halloween for 2007 and 2006.
d5 = datetime('October 31, 2007 12:00 PM','TimeZone','America/New_York')
d5 = 31Oct2007 12:00:00
tzoffset(d5)
ans = 4:00
d6 = datetime('October 31, 2006 12:00 PM','TimeZone','America/New_York')
d6 = 31Oct2006 12:00:00
tzoffset(d6)
ans = 5:00
Now for just a little fun with MATLAB R2014b date calculations. Most of the U.S. changes from Daylight Saving time to Standard Time at 2:00 AM on the first Sunday of November. Given the year, how can we compute that?
Here's just one way. Start with 2:00 AM on November 1 of the specified year.
year = 2014; d = datetime(year,11,1,2,0,0)
d = 01Nov2014 02:00:00
Now use dateshift to move to the first Sunday.
dateshift(d,'dayofweek','Sunday') % You can also use 1 instead of 'Sunday'.
ans = 02Nov2014 02:00:00
My horoscope for today said "Your next blog will include an example of an anonymous function for no good reason." So here it is. Let's make an anonymous function that returns the exact time for the transition from Daylight Saving Time to Standard time in the U.S.
f = @(year) dateshift(datetime(year,11,1,2,0,0),'dayofweek','Sunday')
f = @(year)dateshift(datetime(year,11,1,2,0,0),'dayofweek','Sunday')
f(2014)
ans = 02Nov2014 02:00:00
f(2015)
ans = 01Nov2015 02:00:00
f(2016)
ans = 06Nov2016 02:00:00
Of course, this function only works for years between 2007 and whenever the U.S. Congress changes its mind again!
Next week, I promise, I'll get back to writing about image processing. Or maybe colormaps. Unless I write about something else.
Get
the MATLAB code
Published with MATLAB® R2014b
Brett's Pick this week includes nods to two readers, and a submission called Fast Noise Estimation in Images by Tolga Birdal.
In my last post, back in September, I discussed several files for quantifying noise in images, without reference to a goldstandard image. I Picked and compared three submissions that simplified the quantification of image noise in isolated images. All of my tests, however, were conducted against image files that I had artificially blurred. Eric and Marco correctly took me to task for that, pointing out that 'blur' is not synonymous with 'noise' (though maybe we can agree that it can be a type of noise.)
Marco further pointed me to Tolga's submission, a fiendishly simple function that convolves the test image with a reference kernel to quantify noise.
In today's post, I would like to revisit noise quantification using my previous Picks, and Tolga's, and to look at their application to images containing different types of noise. I will again use the "rice.png" image that ships with the Image Processing Toolbox to consider how each does in quantifying the five types of noise "contant Gaussian," "locally varying Gaussian," Poisson, saltandpepper, and specklecreated by the imnoise function. To start, I'm going to use default input parameters (where possible) in the calls to imnoise.
img = cell(6,1); img{1} = imread('rice.png'); subplot(2,3,1) imshow(img{1});title('No Noise'); img{2} = imnoise(img{1},'gaussian'); subplot(2,3,2) imshow(img{2});title('Gaussian'); img{3} = imnoise(img{1},'localvar',0.05*rand(size(img{1}))); subplot(2,3,3) imshow(img{3});title('Local Gaussian'); img{4} = imnoise(img{1},'poisson'); subplot(2,3,4) imshow(img{4});title('Poisson'); img{5} = imnoise(img{1},'salt & pepper'); subplot(2,3,5) imshow(img{5});title('Salt & Pepper'); img{6} = imnoise(img{1},'speckle'); subplot(2,3,6) imshow(img{6});title('Speckle');
Now consider how each of the tools we've discussed:
quantifies noise in each of these images:
for ii = 1:6 imageBlur(ii) = blurMetric(img{ii}); blindImageQuality(ii) = blindimagequality(img{ii},2); noiseLevel(ii) = NLEstimate(img{ii}); fastNoiseEstimate(ii) = estimate_noise(img{ii}); end
NoiseTypes = {'NoNoise';'Gaussian';'LocalGaussian';'Poisson';'SaltAndPepper';'Speckle'}; Results = table(imageBlur,blindImageQuality,noiseLevelEstimage,fastNoise,... 'RowNames',NoiseTypes');
Interpreting this table is difficult and subjective. We can't necessarily compare metrics (which are on different scales), and we can't assess how each metric reflects different levels of noise. So let's do another series of analyses, increasing the noise successively, and recalculating the metrics. We can then plot some trends.
Let's increment four types of noise (omitting Poisson noise, for which imnoise is not adjustable) five times and recompute the metrics:
noNoise = imread('rice.png'); index = 0; for noiseType = 1:4 % Loop over noise types for noiseLevel = 1:5 index = index + 1; level = 0.01*noiseLevel; % Create noisy images and loop over metrics switch noiseType case 1 str = ['Gaussian (' num2str(level) ')']; img{noiseType} = imnoise(noNoise,'gaussian',0,level); case 2 str = [' Local Gaussian (' num2str(level) ')']; img{noiseType} = imnoise(noNoise,'localvar',level*rand(size(noNoise))); case 3 str = ['Salt & Pepper (' num2str(level) ')']; img{noiseType} = imnoise(noNoise,'salt & pepper',level); case 4 str = ['Speckle (' num2str(level) ')']; img{noiseType} = imnoise(noNoise,'speckle',level); end subplot(4,5,index); imshow(img{noiseType});title(str); imageBlur(noiseType,noiseLevel) = blurMetric(img{noiseType}); blindImageQuality(noiseType,noiseLevel) = blindimagequality(img{noiseType},2); noiseLevelEstimage(noiseType,noiseLevel) = NLEstimate(img{noiseType}); fastNoise(noiseType,noiseLevel) = estimate_noise(img{noiseType}); end end
f = figure; set(f,'DefaultAxesFontsize',8) NoiseTypes = {'Gaussian';'LocalGaussian';'SaltAndPepper';'Speckle'}; for noiseType = 1:4 ax = subplot(4,1,noiseType); noise = [Results{1,1},imageBlur(noiseType,:)]; plot(noise/max(noise)); if noiseType == 1 title('Normalized Noise Metrics','fontsize',12); end hold on noise = [Results{1,2},blindImageQuality(noiseType,:)]; plot(noise/max(noise)); noise = [Results{1,3},noiseLevelEstimage(noiseType,:)]; plot(noise/max(noise)); noise = [Results{1,4},fastNoise(noiseType,:)]; plot(noise/max(noise)); ax.XTick = 1:6; ax.XTickLabel = []; ylabel(NoiseTypes{noiseType}) end ax.XTickLabel = 0:5; xlabel('NOISE LEVEL') legend(MetricNames)
Conclusions are easier to draw now. For the most part, we can see metric values monotonically increasing or decreasing with the amount of noise. (The trend is less important than the monotonicity; we can easily think of the metric values as quantifying image quality.)
One final comment: it's very possiblemaybe even likelythat these metrics can be tweaked to give better results, and that they will perform differently for different image and noise types. I encourage you to play around with the inputs using your own images; hopefully this discussion will give you a framework for evaluating the different noise metrics.
Swag, and thanks, are due to Eric, Marco, and Tolga!
As always, I welcome your thoughts and comments.
Get
the MATLAB code
Published with MATLAB® R2014b
Reuse: Or how I learned to love the Bat Man Lego set…
Like many an engineer of my generation I grew up with the basic Lego® starter set: the 1x1, 2x2, 2x3, 1x8 and the massive 2x8. We considered the wheel a high tech bit of plastic. Years later when I first saw a custom set with its detailed instructions, I was outraged by what I assumed was the corporatesponsored crushing of my nephew’s creative possibilities. I scoffed, I scorned, I worried about the decline of western civilization if we fed these kids ships fully built!
My nephew had a whole raft of “complete” kits. From these kits he learned how to build houses and starships, bridges and mountains. What’s more he learned how to build in a way that was rugged and would survive the hardplay action of space ships crashing into erupting volcanos. One day I stopped and saw what he had built — a hybrid submarine from bits of the Bat Man street cruiser and the Millennium Falcon brought together with good old basic Legos. He christened it the Subbatfalc! It was awesome.
Figure 1 Not actual Subbatfalc
What's the link with Simulink?
So what does this have to do with MATLAB and Simulink? Well, everything. As engineers creating something new, we need to figure out where we start. Do we use the primitive blocks and create everything from first principals or do we find existing components and build from there?
MATLAB and Simulink gives you a wonderful mixture of basic and advanced blocks. Further, with objectoriented programing in MATLAB and Simulink Libraries and Model blocks in Simulink, we give you the ability to create your own advanced blocks (heck, if you use TLC we give you the basic plastic to craft anything).
However as a consultant in the Controls Design Area I have lost track of the number of times I have worked with a client who has built up a set of custom functions such as transfer functions, integrators, or table lookup algorithms. They’ve recreated basic blocks. Now, truth be told, in many cases these custom blocks had a slight performance edge over the builtin blocks. But customers also have to maintain and validate these blocks, taking time away from building something more important.
So how do you determine when to build from scratch and when to use builtin functionality? I would suggest the following 6 State chart.
Figure 2: Note remember to thank Sally!
My nephew took less than an hour to build the Subbatfalc. Having worked with the kits, he had the skills to build it from basic blocks but it would have taken him ten times longer. He made the rational choice to reuse; and so am I. So, here I am recanting my harangue against kits. It isn’t the end of civilization — it’s the start of more interesting time. Now, pardon me, I am going to go help him build the Mount Doom/Death Star mashup (Doom Star??) so we can see what happens when the Subbatfalc attacks!
Now it's your turn
How do you decide between using an existing complex block, versus building it from basic blocks? Let us know by leaving a comment here.
]]>The Method of Particular Solutions computes a highly accurate approximation to the eigenvalue we have been seeking, and guaranteed bounds on the accuracy. It also provides flexibility involving the boundary conditions that leads to the MathWorks logo.
In 196566, after finishing grad school, I spent a postdoctoral year at the ETH, the Swiss Federal Institute of Technology in Zurich. At the time, the faculty included three world renowned numerical analysts, Eduard Stiefel, Heinz Rutishauser, and Peter Henrici. Moreover, IBM had a branch of their research division in Ruschlikon, a Zurich suburb, where they held occasional seminars.
In October, Leslie Fox, from Oxford, visited the IBM center and gave a talk about the work two of his students, Joan Walsh and John Reid, were doing on the Lshaped membrane. Expansions involving Bessel functions that matched the singularity of the eigenfunction near the reentrant corner were coupled to finite difference grids in the rest of the domain. Shortly after Fox's talk, it occurred to me to use the Bessel function expansions over the entire domain and avoid finite differences altogether.
Within a matter of weeks I had done computations on the CDC 1604 mainframe at the ETH that had me very excited. To celebrate, I made an Lshaped ornament out of aluminum foil for the Christmas tree in our apartment.
Henrici contributed some theoretical background relating what I was doing to work of Stephan Bergman and a Russian mathematician I. N. Vekua. Together with Fox, we published a paper in the SIAM Journal on Numerical Analysis in 1967.
You don't need to know much about Bessel functions to appreciate the method of particular solutions. Using polar coordinates, a particular solution is any function of the form
$$ v_\alpha(\lambda,r,\theta) = B_\alpha(\sqrt{\lambda} r) \sin{(\alpha \theta}) $$
First of all, notice that the $\sin{(\alpha \theta)}$ portion insures that this solution satisfies the boundary condition
$$ v_\alpha(\lambda,r,\theta) = 0 $$
along the two rays with $\theta = 0$ and $\theta = {\pi}/{\alpha}$. For our Lshaped domain the reentrant corner has an angle of $\frac{3}{2}\pi$, so we will want to have
$$ \alpha = \frac{2}{3}j $$
where $j$ is an integer.
The Bessel function $B_\alpha$ is defined by an ordinary differential equation that implies any $v_\alpha$ satisfies the partial differential equation
$$ \Delta v_\alpha + \lambda v_\alpha = 0 $$
If $\alpha$ is not an integer, then $B_\alpha(r)$ has another important property. Its derivatives blow up as $r$ goes to 0. This asymptotic singularity is just what we need to match the singularity in the eigenfunction.
We can also take advantage of symmetries. The first eigenfunction is symmetric about the centerline of the L. We can also eliminate components that come from the eigenfunctions of the unit squares that are subdomains. This implies that for the first eigenfunction $\alpha = \frac{2}{3}j$ where $j$ is an odd integer that is not a multiple of 3, that is
$$ \alpha_j = \frac{2}{3}j, \ \ j = 1,5,7,11,13, ... $$
Form a linear combination of particular solutions.
$$ u(\lambda,r,\theta) = \sum_\alpha c_\alpha v_\alpha(\lambda,r,\theta) $$
This function satisfies our differential equation throughout the domain, satisfies the boundary condition on the two edges that meet at the reentrant corner, and has the desired singular behavior there. All we have to do is pick $\lambda$ and the coefficients $c_\alpha$ so that $u$ satisfies the boundary condition, $u = 0$, on the remaining edges of the L.
The general idea is to pick a number of points on the outer boundary and investigate the matrix obtained by evaluating the particular solutions in the linear combination at these points. Vary $\lambda$ until that matrix is nearly rank deficient. The resulting null vector provides the desired coefficients.
The details of how we carry out this general plan have varied over the years. The best approach is implemented in the function membranetx, available here, described in the PDE chapter of Numerical Computing with MATLAB. It uses the matrix computation workhorse, the SVD, the singular value decomposition.
The number of rows in the matrix involved is the number of points on the boundary and the number of columns is the number of terms in the sum. The figures below come from matrices that are 46by9. Let $(r_i,\theta_i)$ be the points chosen on the boundary. The matrix elements are
$$ a_{i,j}(\lambda) = B_{\alpha_j}(\sqrt{\lambda} r_i) \sin{(\alpha_j \theta_i)} $$
The function fmintx is used to vary $\lambda$ and find a minimum of the SVD of this matrix. The resulting matrix is essentially rank deficient and the resulting right singular vector is a null vector that provides the coefficients.
The minimum is found at
$$ \lambda = 9.639723844 $$
Here are the first four particular solutions. The vertical scale reflects the magnitude of the corresponding coefficient.
mps_figs
None of these four terms is remotely close to zero on the outer boundary, but when they are added together the sum is zero to graphical accuracy.
mps_figs(4)
So this is actually the first eigenfunction of the Lshaped membrane.
Now we can have some fun. Here is a surf plot of the sum of just the first two particular solutions.
mps_figs(2)
This function is an exact solution of the partial differential equation
$$ \Delta u + \lambda u = 0 $$
with the value of $\lambda$ reported above. It satisfies the boundary condition $u = 0$ on the two edges meeting at the reentrant corner, but it fails significantly on the outer edges. So it is only with some artistic license that we can refer to it as the first eigenfunction of Lshaped membrane.
When I first saw this surf plot, it didn't have this gorgeous color map, but I found the shape so attractive that I kept it around. The evolution of the shading, lighting and colors that turned this into our logo is the subject of my next blog.
The Fox, Henrici and Moler paper was not just about the method of particular solutions. We were also concerned about finding guaranteed upper and lower bounds for the exact eigenvalue of the continuous differential operator. A theorem proved in the paper relates the discrepancy in meeting the boundary condition to the width of a small interval containing the eigenvalue. In 1966 we proved
$$ 9.6397238_{05} < \lambda < 9.6397238^{84} $$
Today with membranetx on MATLAB, I am finding a value of $\lambda$ in the center of this interval.
[L,lambda] = membranetx(1);
format long
lambda
lambda = 9.639723844021997
In 2005 Timo Betcke and Nick Trefethen published a paper in SIAM Review that pointed out a serious shortcoming in the approach we described in the Fox, Henrici and Moler paper. When the method is attempted on other domains, especially domains with more than one singular corner, the matrix has a tendency to appear nearly rank deficient for all values of $\lambda$. It is difficult to identify local minima. The recommended cure is to add a number of interior points and require the approximate eigenfunction to be nonzero there.
Betcke and Trefethen also prove an approximation theoretic result that, when combined with their computations, adds four more digits to our rigorous bounds on the first eigenvalue.
$$ 9.63972384402_{16} < \lambda < 9.63972384402^{22} $$
Moreover, they are so confident in their numerical results that they are sure the exact eigenvalue lies in the center of this interval, and they are pretty sure they even know the next digit.
$$ \lambda \approx 9.63972384402195 $$
This agrees with what I get out of membranetx to almost full double precision accuracy. I think we finally have this fellow pinned down.
I recommend that anyone reading this blog who wants to know more about the method of particular solutions and computing the eigenvalues of the Lshaped membrane read the Betcke and Trefethen paper.
Leslie Fox, Peter Henrici, and Cleve Moler, Approximations and Bounds for Eigenvalues of Elliptic Operators, SIAM Journal Numerical Analysis 4, 89102, 1967, <http://blogs.mathworks.com/images/cleve/Fox_Henrici_Moler.pdf>
Timo Betcke and Lloyd N. Trefethen, Reviving the Method of Particular Solutions SIAM Review 47, 469491, 2005, https://people.maths.ox.ac.uk/trefethen/publication/PDF/2005_112.pdf
Cleve Moler, Numerical Computing with MATLAB, Chapter 11, Partial Differential Equations, 2004. <http://www.mathworks.com/moler/pdes.pdf>
Get
the MATLAB code
Published with MATLAB® R2014b
I hesitated for some time today about where to go next in my discussion about the new default colormap in MATLAB. I finally decided to spend a little time on the reactions so far. Those reactions have ranged from "Great, what took you so long?" to "I hate it; give me back jet!" In the Twitterverse, one scientist proclaimed "End of an era!" Another simply said "yuck." (Jet is still in MATLAB, by the way.)
Today I want to dig into some of those reactions a little bit. Along the way I'll show some experiments with a test pattern to explain some of my comments about jet.
To remind you, last time I summarized 15 years or so of criticisms about rainbow colormaps:
I have seen a claim that that rainbow incorporates a universally understood order based on associating blue with cold and red with hot. There's no doubt that this order makes sense for some people (although blue flame is hotter than red flame?). I quibble with the "universally understood" part, though. The dominant evidence is that individuals do not universally order the spectral hues the same way unless they refer to spectral order mnemonic, such as "Roy G. Biv." It is wellestablished, however, that individuals can easily and consistent order colors by relative brightness.
Another commenter said that the rainbow colormap has higher contrast than parula. Yes, this is true. Rainbow goes all the way from dark to bright in the bottom third of the data range, and it goes all the way from bright to dark in the top third of the data range. For some parts of the data range, therefore, there is higher visual contrast than with parula, which uses the entire data range to go from dark to bright. A price paid for this increased contrast, however, is the visual illusions and confusion caused by the reversed roles of dark and bright in different parts of the colormap.
Matteo Niccoli gives the following nice example of a rainbowbased visual illusion. The arrow in the righthand image below points to something that gives the strong but false impression of being a depression in the bone structure.
Image credit: Matteo Niccoli, "The Rainbow Is Dead ? Long Live the Rainbow!" Copyright 2014 MYCARTA. Used with permission.
And I showed another visual confusion example last time with these visualizations of two different sinusoidal patterns:
These visualizations give the impression that the two sinusoidal patterns have opposite phase, but in reality they have exactly the same phase.
Another commenter said that, with sufficiently familiarity, the bright "stripes" that are characteristic of the rainbow colormap can be an advantage for certain kinds of data. That might be true in some particular cases, but I don't think it's true in general. Let me show an example, using a test pattern, that illustrates how much difficulty the cyan and yellow stripes can cause in interpreting what should be easily recognizable data.
url = 'http://blogs.mathworks.com/images/steve/2014/eia2000cropped.png';
I = imread(url);
imshow(I)
colormap(jet)
Let me label some portions of the visualization and then ask some quetions about it.
text(35,90,'A','FontSize',25,'HorizontalAlignment','center',... 'BackgroundColor','white') text(125,90,'A','FontSize',25,'HorizontalAlignment','center',... 'BackgroundColor','white') text(215,90,'A','FontSize',25,'HorizontalAlignment','center',... 'BackgroundColor','white') text(190,320,'B','FontSize',25,'HorizontalAlignment','center') line([77 160],[225 225],'Color','white','LineWidth',5) text(65,225,'C','FontSize',25,'HorizontalAlignment','right',... 'BackgroundColor','white')
At periodic intervals, the dark blue vertical bars have a bright outline (see the points marked "A"). What is the data significance of these bright outlines? Especially when compared with the dark blue vertical bars that do not have a bright outline?
The vertical bars in between the "A" marks don't have bright outlines. Those bars alternate between dark blue and dark red, a lowcontrast combination that is a little difficult to focus clearly on (a redblue visual illusion called Chromostereopsis).
Above the "B", I see very thin vertical dark bars. When the bars are this thin, it's difficult to distinguish between the dark red and dark blue, which is unfortunate because those colors represent opposite ends of the data range. The bright outlines between the bars is broken up by a diagonal pattern of blue stripes.
The numbers around the border are oddly broken up and hard to decipher.
Let's investigate why some vertical bars are outlined and some not by looking at the data crosssection along the white path marked "C".
c = improfile(I,[77 160],[225 225]);
figure
plot(c,'*')
I'll superimpose two lines that show the data values corresponding to the bright cyan and bright yellow in the jet colormap.
hold on plot([1 84],[1 1]*0.375*255,'c','LineWidth',5) plot([1 84],[1 1]*0.625*255,'y','LineWidth',5) hold off
It is now apparent why some vertical bars are outlined but not others. In each transition from low to high and from high to low, there is one data value that is in between. When that one transitional data value happens to be close the data values mapping to cyan or yellow, the vertical bar is shown with a bright outline. Otherwise, there is no outline. Although the effect is visually dramatic and makes some regions of the data look quite different from others, the visual difference has no real data significance.
OK, now let's show this data using the parula colormap.
imshow(I) colormap(parula)
This visualization is trivially easy to interpret. The vertical bars in the middle are all the same. The digits are easy to read. The visual difference between low data values and high data values has high contrast and is easy on the eyes. And the data above the point marked "B" is seen to have a completely different form than the way it appeared with jet. The stripes have a lower spatial frequency and are seen to be slightly fanning out in direction instead of appearing vertical.
A valid complaint about this example is that it is completely artificial, with data rather unlike what people are likely to be looking at with jet or parula. I'd like to suggesting thinking about it this way, though: If jet can take familiar, easily interpretable data and render it almost unrecognizable, what can jet do to your own data?
Finally, I'll mention that a blog reader commented that using jet made it easy to distinguish positive and negative values from each other in a data set whose values are centered around zero. And yes, jet can be better than parula for that. Jet is bright and relatively unsaturated in the middle, and it is dark and saturated at both ends. These are characteristics often shared by colormaps designed specifically for highlighting a central, specific value and distinguishing the values above and below that. These are called divergent colormaps. In a future post I'll show where to find some good divergent colormaps and how to use them in MATLAB.
Get
the MATLAB code
Published with MATLAB® R2014b
Enabling Interface Display
In R2014b, you will notice a new entry at the top of the Display menu:
For a model that originally looks like this:
Enabling the interface display makes it look like this:
As you can see, this feature puts the emphasis on the inputs and outputs of the system.
Advantages of the Interface Display
For someone like me who receives new large Simulink models every day, this feature is very useful. When I explore new models, I can use the highlighting to trace sources and destinations of signals through multiple levels of nested subsystems. As you can see in the following animation, by zooming in, I can select one element of a bus signal and follow the highlighted signals to quickly find out where it is used.
Once the model has been simulated or updated once, the interface view displays for each port its data type and sample time, making it easier to understand the model.
One more thing to note. When in interface Display, you cannot edit the model. This is convenient when your goal is just to explore and understand the model without modifying it.
Now it's your turn
Give the new interface view a try and let us know what you think by leaving a comment here.
]]>Today, David Garrison, our guest blogger, will continue his series on the new graphics system in R2014b.
Here is Part 3 of the series.
In Part 1 of this series, I provided an introduction to the new MATLAB graphics system in R2014b with descriptions of a number of new features. In Part 2, I talked about one of the big changes in R2014b  graphics functions return MATLAB objects, not numeric handles. In this post, I will talk about compatibility considerations in using R2014b graphics.
This post is longer than my two previous posts. In the first section, I'll talk about the kind of changes we've made in R2014b graphics, why we made those changes, how they affect you, and where to go for more information. The second section describes visual changes in the new graphics system including a new colormap, new line colors in plots, and new sizes for plot titles and axis labels, and these affect everyone. The last section is about changes we've made that affect advanced graphics users and user interface builders. Most of you can skip the more advanced material in the last section.
When we talk about compatibility considerations we are talking about changes in R2014b that may not be compatible with code written before R2014b. These incompatibilities can affect your code in one of three ways.
First of all, let me say that here at MathWorks we take compatibility very seriously. We have a rigorous process for vetting potential changes in MATLAB that may affect existing MATLAB code. We evaluate them very thoroughly to determine the extent of their impact. We document those changes in our release notes and provide alternatives whenever possible. Those changes are necessary to allow us to continue to develop MATLAB.
One of our goals for the new graphics system was to introduce a new architecture that would allow us to build new graphics capabilities for the future. We couldn't continue with the old graphics architecture because it was too limiting. The new architecture has also allowed us to fix a large number of outstanding graphics bugs that have accumulated over the years.
Our second goal was to support existing graphics functionality and to minimize the disruption to MATLAB users. We did that by evaluating every change that would introduce an incompatibility in the graphics system. We tested those changes in house with over 100,000 MATLAB files and in person with a large number of advanced graphics users to determine the following:
The results of all this testing allowed us to do two things. First, it allowed to make changes in the new graphics system before the release to minimize the effect of certain incompatibilities. Second, it told us what kind of documentation and tools would be required to help people make the transition.
Anyone using graphics will notice a visual difference in R2014b. People who've seen R2014b tell us the new visual appearance is a big improvement over previous versions of MATLAB. Differences include changes to colormaps, line colors, font sizes, grids, and other visual properties. For most people that's it. That's all you need to know.
Those of you who do advanced graphics programming and build user interfaces may be affected by incompatibilities that will require changes to your code. I'll talk about the most common ones in the last section of this post. We've also provided plenty of documentation and some tools to help you make the transition.
In preparation for the release of the new graphics system, we created a lot of material to help people understand the changes and make any modifications to their code. Here is a list of resources for information on R2014b graphics:
Anyone who uses MATLAB graphics in R2014b will see that plots look quite different. In this section, I will show you some of the big visual differences. As we go through these examples, I think you'll agree that the changes in R2014b are a big improvement over previous MATLAB versions. Of course, there may be times when you may want to change some of these visual properties. I will either tell you how to do that or I will provide a link to the documentation for more information.
In MATLAB, the colormap is a matrix of RGB color values used to set the colors for images and surfaces. For example, when you call the surf function you get a plot where the surface color is proportional to the height of the plot. Colors are based on the values from the colormap.
In R2014b, we introduced a new colormap called parula which is now the default. In previous versions of MATLAB, the default colormap was jet. Here's a comparison of the colormaps between R2014a and R2014b.
We changed the default colormap to address problems with rainbow colormaps like jet. The new parula colormap is ordered from dark to light and is perceptually uniform. For more information about the issues with rainbow colormaps like jet, see the post on Steve Eddins blog called A New Colormap for MATLAB – Part 2 – Troubles with Rainbows. You can use the colormap function to change the colormap to whatever you want it to be. For example, to change the colormap to jet use the following command after you create your plot.
colormap(jet)
We've also changed the lines colors used for plots. These line colors were selected to make lines easier to distinguish and to help people with certain types of color blindness.
The line colors used in plots are controlled by the ColorOrder property of the Axes object. One other note regarding line colors in plots. In R2014b, when plot is called with hold on, MATLAB will use the next color in the color order for each call to plot. That gives plots with different line colors. In previous versions, each call to plot would start over in the color order so that each line would be the same color.
For information on controlling the colors in the ColorOrder see the documentation section Why Are Plot Lines Different Colors?
One last change regarding appearance of plots. In R2014b, plot titles and axis labels use a larger font size than in previous versions of MATLAB and plot titles are bold by default.
I intentionally used a very long title in this plot to make a point. In this case, the title is so long it extends beyond the boundaries of the axes. You can control your plot title and axis labels using the TitleFontWeight, TitleFontSizeMultiplier, and LabelFontSizeMultiplier properties. For more information see How Do I Make the Graph Title Smaller?
That covers the major visual changes for plots in R2014b. The next section covers changes that only affect advanced graphics users and GUI builders.
Most of you can stop reading here and get on with the rest of your day.
If you are still reading, I assume you do some advanced graphics programming in MATLAB or maybe you build MATLAB user interfaces. In this section I'm going to cover the four most common nonvisual changes in R2014b graphics that you are likely to encounter. If you do encounter one of these, you will probably have to make changes to your code.
I mentioned this briefly in Part 2 of this series. Graphics functions now return objects, not numeric handles. The R2014b documentation has detailed information about this subject in the section called Graphics Handles Are Now Objects, Not Doubles. I will give a couple of simple examples to illustrate what can happen with code written before R2014b.
Prior to R2014b, you could store a set of handles to graphics objects in an array and then add some numeric data to that array. In R2014b, that will cause an error.
x = pi:0.1:pi ;
y1 = sin(x);
y2 = cos(x);
myLines = plot(x,y1,x,y2) % plot returns an array of two Line objects
If you then try to set myLines(3) = 1.2, you get the following error.
Cannot convert double value 1.2 to a handle
MATLAB won't let you add numeric values to an array of graphics objects. A similar problem occurs if you try to use an object handle in a function where MATLAB expects a numeric value. A simple example of this happens with the sprintf function.
a = sprintf('You clicked on figure %d\n', gcf);
The %d specification in the sprintf format string expects an integer value. However, since gcf is a figure object, you get the following error.
Error using sprintf Function is not defined for 'matlab.ui.Figure' inputs.
Here is one final example. Because graphics handles used to be numbers, you could use them in logical expressions.
if (get(0, 'CurrentFigure')) disp(['Figure ' get(gcf, 'Name')']) % display the figure name for gcf else disp('No open figures') % there is no open figure end
This worked in earlier versions of MATLAB because get(0,'CurrentFigure') would return either an empty array or a numeric figure handle. Both of these values are valid in the logical test of the if statement above. In R2014b, this will cause an error.
Conversion to logical from matlab.ui.Figure is not possible.
We have tried to maintain compatibility with previous releases in some cases. For example, you can still use 0 to refer to the graphics root in functions like get and set. As a best practice, however, we now recommend using the groot function to get the graphics root. Similarly, we still support the use of literal integer values to refer to figures in functions like set, get, and figure. Again, the best practice is to use a variable which contains the object when using these functions.
If you find yourself really stuck, it is possible to cast object handles to numeric handles using the double function. You can then cast the number back to an object handle using the handle function. We don't recommend this as a long term solution. Be aware that we may choose to remove this feature in a future version of MATLAB. If we do, we'll let you know in advance.
In earlier versions of MATLAB, the legend and colorbar functions created objects whose type was 'axes'. Technically, legends and colorbars were subclasses of axes. In R2014b, the legend function creates a Legend object and the colorbar function creates a ColorBar object.
So what does that mean? First, the results of a findobj call in R2014b may be different than in previous versions. In R2014b, the command
findobj('Type', 'Axes')
will not return any legend or colorbar objects. To get those objects you will need to do the following:
findobj('Type', 'Legend') findobj('Type', 'ColorBar')
Second, legends and colorbars cannot become the current axes. Prior to R2014b, you could write code that looks like this:
plot(1:10) cb = colorbar; axes(cb) % Make the colorbar the current axes title('My Colorbar') % set the title of the colorbar
If you try to run this code in R2014b you will see an error.
Error using axes Handles of type ColorBar cannot be made the current axes
In R2014b, you can still use the title function. Just give it the colorbar handle as the first argument.
title(cb, 'My Colorbar')
Alternatively you can use the colorbar's Label property to get the text object for the label and set its string property.
cb.Label.String = 'My Colorbar'
Charting functions like bar, contour, stem and others return one or more graphics objects. In previous versions of MATLAB, those objects had names like barseries, contourgroup, and stemseries. These objects were each a special type of object called an hggroup. Each of these hggroup objects had a set of children. The children were low level objects like lines or patches. Here is an example from R2014a.
[~,h] = contour(peaks, 'LineLevels', 6:1:8) ; get(h, 'Type') % an hggroup object is returned by the contour function
ans =
hggroup
ch = get(h, 'Children') ; get(ch(1), 'Type') % the children of the hggroup are patch objects
ans =
patch
Prior to R2014b, the contour command returned an hggroup object that had children that were patch objects (one for each contour line). Using the patch objects, it was then possible to do interesting things with the contour lines. You could, for example, make the even numbered contour lines solid and make the odd numbered contour lines dashed as shown below.
As you've probably guessed, things are different in R2014b. First, the types of the objects returned by these functions are different. For example, the bar function returns a Bar object and the contour function returns a Contour object. Second, the objects returned by these functions no longer have any children. That means you can't get to the lowlevel objects.
So what do you do? You have to take a different approach. In the code below, I've create two contours  one with solid lines and one with dashed lines. This code will create a plot like the one shown above. It works in R2014b and in previous versions of MATLAB.
major = 6:2:8; minor = 5:2:7; [~,hmajor] = contour(peaks,'LevelList',major); % contour with evennumbered levels
hold on [~,hminor] = contour(peaks,'LevelList',minor); % contour with oddnumbered levels set(hminor,'LineStyle',':') % make the oddnumbered levels dotted hold off
There is one last change that I want to talk about. It can affect MATLAB GUIs created before R2014b. Suppose I have a simple GUI with a panel and a button. I might have code that looks something like this:
figure('NumberTitle', 'off', 'Name', version('release'), ... 'Position', [100 100 350 260], 'MenuBar', 'none', 'Toolbar', 'none') ; button = uicontrol('Style', 'pushbutton', 'String', 'My Button', ... 'Units', 'normalized', 'Position', [0.4 0.5 0.25 0.15]) ; panel = uipanel('Position', [0.10 0.10 0.8 0.8], 'Title', 'My UIPanel') ;
Seems OK right? Here's what it looks like in R2014a and R2014b.
So what happened here? I don't see my button in R2014b. Where did it go? It turns out that the button is still there but in R2014b it is drawn under the panel.
In previous versions of MATLAB, uicontrols were always drawn on top on uipanels regardless of the order in which they were created. In R2014b, the two components are drawn in creation order. Since the panel was created after the button, it is drawn on top. What I really want is to have the button be a child of the panel. In order to do that, I can change my code to look like this:
figure('NumberTitle', 'off', 'Name', version('release'), ... 'Position', [100 100 350 260], 'MenuBar', 'none', 'Toolbar', 'none') ; panel = uipanel('Position', [0.10 0.10 0.8 0.8], 'Title', 'My UIPanel') ; button = uicontrol('Parent', panel, 'Style', 'pushbutton', 'String', 'My Button', ... 'Units', 'normalized', 'Position', [0.4 0.5 0.25 0.15]) ;
You can also fix this problem in GUIDE by moving the button out of the panel and then moving it back in. This will automatically make the button a child of the panel.
You can check your existing GUIs by downloading a tool from Tools for transitioning to R2014b Graphics. For more information, see Why Are Some Components Missing or Partially Obscured?
Read the materials I've listed above to learn more about the changes to graphics in R2014b and how they might affect your MATLAB code. Download the tools and run them on your code and your GUIs. If you aren't able to find the information you need to work around one of the changes in R2014b, contact our Technical Support Team . They are very familiar with the changes in R2014b and can help you find a solution.
Have you encountered any of these changes in R2014b? Have you been able to find the information you need? Have you been able to use this information to make necessary changes to your code? Have you tried the transition tools? We'd love to hear your thoughts here.
I am hopeful that this series of posts has given you a good introduction to the R2014b graphics system. There's a lot of new stuff to digest in this release. Try these things out and let me know what you think.
Thanks to Loren for letting me take over her space for the last few weeks.
Get
the MATLAB code
Published with MATLAB® R2014b
The Partial Differential Equation Toolbox contains tools for the analysis of PDEs in two space dimensions and time. It is perhaps not surprising that one of the primary examples involves the Lshaped membrane.
If you have the PDE toolbox installed, bring up pdetool. Click on the first triangle in the toolbar. This initializes a tool for creating a planar domain with an unstructured mesh. The default domain is our L. The default mesh has equally spaced points on the boundary connected by an irregular grid in the interior with 150 nodes and 258 triangles.
The PDE tab allows the specification of a general elliptic operator with variable coefficients, but if we accept the defaults we have our desired Laplacian. Here is the first eigenfunction, obtained with the coarse grid, plotted with the default cool colormap. The first eigenvalue, reported in the title, is 9.9707.
Refine the grid three times to one with 8417 nodes and 16512 triangles, change to our new parula colormap, and add contour lines. The reported eigenvalue is 9.6525.
The pdetool has an option to plot the absolute value of the gradient. We see that the singularity at the origin acts like a black hole, sucking in all the color.
The mesh refinement is accomplished by adding grid points at the midpoints of the triangles, thereby quadrupling the number of triangles. It is possible to do this five times, producing a fine mesh with 258 * 4^5 = 264192 triangles. (Actually the mesh can be refined six times to one with over a million triangles, but the subsequent eigenvalue calculation runs out of memory.)
Here are the results for the first eigenvalue of the Lshaped membrane obtained with successive bisections of the unstructured triangular mesh.
format long load pdetool_results L
L = 9.970747465677787 9.745769128365856 9.675647712787020 9.652453140975609 9.644395207950412 9.641482807142362
These values are candidates for acceleration by Aitken's deltasquared process.
$$ x_n  \frac {(\Delta x_n)^2}{\Delta^2 x_n} $$
d1 = diff(L,1); L2 = L(3:end)  d1(2:end).^2./diff(L,2)
L2 = 9.643895738979518 9.640988739828559 9.640105597452624 9.639834371546435
We can even use deltasquared a second time, although this is not often justified.
t1 = diff(L2,1); L2(3:end)  t1(2:end).^2./diff(L2,2)
ans = 9.639720224107141 9.639714153353552
Compare this with the value I reported in my previous post, obtained by extrapolating from a regular square grid.
lambda1 = 9.639723753239963
We will see in my next post a much better way to compute the eigenvalue that will provide a result to nearly full double precision accuracy.
Get
the MATLAB code
Published with MATLAB® R2014b
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 highquality 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!
]]>If you're going to be there, come see my MATLAB Techniques for Image Processing talk on Wednesday, 12:45  13:45. Or stop by the MathWorks booth.
Now you'll have to excuse me while I finish preparing my talk.
]]>Scrolling
The Scope block now has a new scroll option to control the behavior when the duration of the simulation is longer than the Scope time range. The best way to describe that is to show it:
Signal Viewer and Floating Scope Look and Feel
In R2014b, the Signal Viewer and Floating Scope interfaces are now identical to the Simulink Scope block. Whether you prefer having a block or not in the model, you can experience the same visualization:
This means that the Viewer line properties, axes and figure colors can now be customized. Here is an example of how I like to configure it:
Signal Viewer supports multiple new features
For a few releases, the Signal Viewer and Floating Scope have lagged behind with many recent Simulink features. Not anymore! In R2014b, you will notice that the Signal Viewer and Floating Scope now offer features similar to the Scope block and support:
Now it's your turn
Do you prefer to use the Scope, Floating Scope, or Viewer? Why? Let us know by leaving a comment here.
]]>The first one I saw was Cleve's followup post on the eigenvalues and eigenfunctions of the Lshaped membrane, the basis of the MathWorks logo.
The last image from Cleve's post is a lovely contour plot of one of the eigenfunctions of the shape, rendered using the new parula colormap I've been discussing this month.
The second was the hilariously overthetop hightech Friday afternoon cookie detector created by Tom Lowell and posted on Loren's blog. Tom uses an IP camera on the ceiling and Image Processing Toolbox functions and the new webread function in R2014b and some File Exchange contributions to solve a difficult problem faced by MathWorkers in Natick, Massachusetts at the end of every week.
]]>
Today, I'd like to introduce a guest blogger, Tom Lowell, who is a program manager here at MathWorks. He works in our hardware connectivity group, and plays with Arduinos and Raspberry Pis in his spare time. He's relatively new to both MathWorks and MATLAB and decided to write a cookie detector as his first MATLAB project.
When I joined MathWorks, one of the first traditions I learned was that a large plate of home made cookies is delivered every Friday to every floor of every building here in our corporate headquarters in Natick. The problem is, they go fast! And you don't know when they'll be delivered. So if you're in a meeting or hunkered down at your desk, it's easy to miss out.
This seemed like a great opportunity to get some handson experience with MATLAB's image processing capabilities and solve an important problem. I'll describe the cookie detector and how it's built. And hopefully this will whet your appetite to go off and explore MATLAB's image processing capabilities on your own.
The cookie detector does one thing. It detects when a large circular plate of circular cookies is placed anywhere on the counter in the kitchen area on my floor. The system is very simple with only two components; an IP camera in the kitchen and MATLAB running on my computer.
IP cameras are similar to web cams except they run independently from any particular computer. They're attached to your network just like your laptop, your printer, or your smart phone. Most IP cameras are accessible via their IP address. Many, including mine, have a tiny web server running on their chipset. You just type their IP address into your web browser and you get back a JPG image or even video. The IP camera that I happen to be using is a Foscam FI8910W. But you can use any IP camera that supports URL addressing. Your MATLAB code will need to know the IP address of your IP camera. Therefore, it's desirable to configure your IP camera with a static IP address so it doesn't change over time.
You need a nearby outlet to power the IP camera, and a network connection. Most IP cameras support both WIFI or hardwired connections. We're using a WIFI connection because the computer that's running MATLAB is on the other side of the building. The ceiling in our kitchen happens to be metal. This made mounting the camera a breeze (using magnets).
All code for this project is contained within MATLAB. You will you also need the Image Processing Toolbox. The code is very simple. The main function loops forever, until a plate of cookies is detected. It downloads a new image from the camera every 30 seconds. MATLAB analyses the image using its builtin function imfindcircles. Once the cookies are detected, it sends a text message to several members on our team and then exits.
Getting the image from the camera into MATLAB was is the easiest part of the project. The camera has a tiny, builtin web server running on it that will return an image when requested via a URL. And MATLAB has several functions that will fetch data from web servers. I used webread which is new in R2014b. It reads content from a RESTful web service and returns it as Internet media types such as JSON, XML, images, or text.
Create a string with the camera's IP address plus specific API commands to pass to the camera. The format of this URL is camera specific.  snapshot.cgi tells the camera that we want a single image This camera also requires user and pwd as RESTful namevalue pairs.
IPcamURL = 'http://172.31.28.24:1026/snapshot.cgi'; user = 'cookies'; password = 'lovecookies'
Use a timeout so that network issues don't cause it to hang.
options = weboptions('Timeout', 10);
The IP camera's web interface returns content as a JPG image stream. webread converts that for you into an RGB image that is immediately usable by MATLAB functions.
RGB = webread(IPcamURL, 'user', user, 'pwd', password, options); imshow(RGB);
The analysis for cookies is also pretty straightforward. The code first checks the image for a large circle (the round cookie plate) and then checks the image for smaller circles (cookies) contained within the large circle. My code calls the Image Processing Toolbox function imfindcircles, passing it the image from the camera, along with a range of valid radii to check for. Any results are returned in a vector of circles.
Check for a plate. If one or more plates are found, check for cookies... detectPlateImfind() and detectCookiesImfind() are just customized wrappers around imfindcircles
[plateCenter, plateRadius] = detectPlateImfind(tempJpgFile, smallPlateRad, bigPlateRad, 'g', 1);
nCookies = 0; nPlates = size(plateRadius, 1); if nPlates > 0 nCookies = detectCookiesImfind(tempJpgFile, smallCookieRad, bigCookieRad, 'r', 0, plateCenter, plateRadius); end
The analysis takes about 1 second per image on my MacBook Air.
Writing the basic code was easy because MATLAB insulates the programmer from the complexities of the various imageprocessing algorithms. However, you can (and need to) tune your code by passing various parameters to these image processing algorithms. I'll briefly describe the type of tuning I did here. For more information, there's a great example of tuning imfindcircles in the MathWorks documentation center.
First, there are few things to notice about the real life image above. The plate is a perfect circle with welldefined edges, a known radius, and high contrast relative to the counter top. Also, the cookies are generally circular but they're a bit irregular. They may be different colors and sizes each week.
The eventual call to the builtin function imfindcircles takes advantage of all of these. Here's the call to find the plate. The range of acceptable plate radii is passed in (smallest, biggest). The "object polarity" of the silver plate relative to the dark blue counter is passed in as "Bright". And some other parameters are set that were determined by tuning with a positive test case.
Find the plate(s)
[cirCen, cirRad, cirMetric] = imfindcircles(grayimg,[smallest, biggest],... 'ObjectPolarity','Bright', 'Method', 'Twostage','Sensitivity',0.96);
Notice that the parameters to the second call to imfindcircles to find the cookies are different. The object polarity is "dark" because the cookies are dark relative to the silver plate. A different (default) detection method was used and an edge threshold was used.
Find any cookies (on the plates)
[cirCen, cirRad] = imfindcircles(grayimg,[smallest, biggest],... 'ObjectPolarity','dark','Sensitivity',0.96, 'EdgeThreshold',0.05);
There are tradeoffs to be made when you adjust these parameters. For example, if you want to detect every cookie on the plate, you might have to relax the parameters so much that the algorithm also incorrectly identifies other artifacts in your image as circles (other shapes, arcs, flashes of light). But if you tighten up your parameters to filter out those incorrect detections, you may lose detection of overlapping cookies, irregularly shaped cookies, etc.
The same logic applies to the plate or anything else you might want to detect. We don't want to miss the cookie plate, but we also don't want the tuning to be so permissive that we accidentally send 20 people rushing to the kitchen when someone sets their lunch plate on the counter.
The takeway here is that there is no one set of parameters to imfindcircles or any other image processing algorithm that will guarantee perfect results. You need to tune the parameters based on the real world characteristics of your image (such a the shape and color of the object, the lighting, angle, and other factors). This tuning is iterative and can take a while.
The output is a display of the image with the plate circled in green and the cookies in red. This is very easy to produce with MATLAB's basic plotting functions because imfindcircles returns a vector of circles (including the center points on the image and radii).
Although this is a pretty simple image, there are a few things to notice. It correctly detected plates of two different sizes. However, there are some false detects and it didn't get the exact number of cookies right. Lastly, those red crosses seen outside of the plate boundaries are detections that the loosely tuned algorithm thought were circles. However, some postprocessing code was able to filter those out by determining they were not on the plate. This code can detect the plates with very good accuracy because of their welldefined shape and contrast. The code does not do as well with the cookies, which are irregular and offer less contrast relative to their background.
Remember, the object here is not a perfect cookie detector. The goal is correct and timely notification of the cookies' arrival. And it's pretty good at that. These cookies were detected before they were unwrapped!
Once the code has determined that there is a plate of cookies, it sends a text message to interested team members. I'm not going to cover the text messaging in this blog. It works by sending an email through the unique mail gateway of the cellular carriers. For example, to send a text to a Verizon user, MATLAB can send an email to 'xxxxxxxxxx@vtext.com' (where 'xxxxxxxxxx' is the cell phone number). I used Ke Feng's send text message example that's posted on MATLAB Central. The cool thing is that it's all done in MATLAB.
Do you have a project that uses an IP camera? We'd love to hear how you might use MATLAB with your IP camera project at school, work, or home. Tell us about it here.
Get
the MATLAB code
Published with MATLAB® R2014b
After reviewing the state of affairs fifty years ago, I use classic finite difference methods, followed by extrapolation, to find the first eigenvalue of the region underlying the MathWorks logo.
My Ph. D. dissertation, submitted to the Stanford Math Department in 1965, was titled "Finite Difference Methods for Eigenvalues of Laplace's Operator". The thesis was not just about the Lshaped region, but the L was my primary example.
My thesis advisor was George Forsythe, who had investigated the L in some of his own research. At the time, many people were interested in methods for obtaining guaranteed upper and lower bounds for eigenvalues of differential operators. Forsythe had shown that certain kinds of finite difference methods could provide lower bounds for convex regions. But the L is not convex, so that made it a region to investigate.
Nobody knew the first eigenvalue very precisely. Forsythe and Wolfgang Wasow published a book in 1960, "Finite Difference Methods for Partial Differential Equations" (John Wiley & Sons), that features the L in several different sections. In section 24.3 Forsythe references some calculations he had done in 1954 that led him to estimate that
$$ \lambda_1 = 9.636 $$
He goes on to report an informal conjecture made in the late '50s from his colleague R. De Vogelaere that
$$ 9.63968 < \lambda_1 < 9.63972 $$
I repeat these numbers here to give you some idea of the accuracy we were working with back then. We will see in a later post that the right answer is just outside De Vogelaere's upper bound.
Here's a page from my thesis. The contour plot at the bottom of the page would eventually evolve into the MathWorks logo. The plot was made on a Calcomp plotter, which involved a computer controlled solenoid holding a ballpoint or liquid ink pen moving across a paper rolled on a rotating drum. This was a very effective device and was used in computer centers worldwide for a number of years.
The calculations shown on this page are the first eigenvalue of the finite difference Laplacian for ten values of the mesh size $h$. The finest mesh, $h$ = 1/100, has 29,601 interior points. As I recall, the calculation of the first eigenvalue, which I did with a kind of inverse power method, took about half an hour on Stanford's campus main frame, an IBM 7090.
Here is a plot of these values I computed fifty years ago. I've added a black line at what we now know to be the limiting value, the first eigenvalue of the continuous differential operator. You can see that the finite difference values are converging quite slowly. Furthermore, they are converging from above. Forsythe was not going to get his lower bound for this nonconvex region.
thesis_graph
The show rate of convergence of the finite difference eigenvalue with decreasing mesh size is caused by the fact that the corresponding eigenfunction is singular at the reentrant corner in the region. The gradient becomes infinite at that corner. You can see the contour lines crowding together near the corner. If you tried to make a drum head or tambourine by stretching a membrane over an Lshaped frame, it would rip.
To be more precise, use polar coordinates, $r$ and $\theta$, centered at the reentrant corner. Note that to cover the interior of the region $\theta$ will run from $0$ to $\frac{3}{2}\pi$. It turns out that approaching the corner the first eigenfunction has the asymptotic behavior
$$ u(r,\theta) \approx r^\alpha \sin{\alpha \theta} $$
with $\alpha = \frac{2}{3}$. This implies that the $k$ th derivatives of the eigenfunction have singularities that behave like
$$ \partial^k u \approx r^{\alphak} $$
When this fact is combined with the Taylor series analysis of the discretization error of the fivepoint finite difference operator we find that the first eigenvalue converges at the rate
$$ \lambda_{1,h}  \lambda_1 \approx O(h^{2 \alpha}) = O(h^\frac{4}{3}) $$
This is significantly slower than the $O(h^2)$ convergence rate obtained when the eigenfunction is not singular.
Now let's use my laptop and the sparse capabilities in MATLAB. The functions numgrid, delsq, spy, and eigs have been included in the sparfun directory since its introduction. We begin by generating a small Lshaped numbering scheme.
n = 10;
L = numgrid('L',n+1)
L = 0 0 0 0 0 0 0 0 0 0 0 0 1 5 9 13 17 21 30 39 48 0 0 2 6 10 14 18 22 31 40 49 0 0 3 7 11 15 19 23 32 41 50 0 0 4 8 12 16 20 24 33 42 51 0 0 0 0 0 0 0 25 34 43 52 0 0 0 0 0 0 0 26 35 44 53 0 0 0 0 0 0 0 27 36 45 54 0 0 0 0 0 0 0 28 37 46 55 0 0 0 0 0 0 0 29 38 47 56 0 0 0 0 0 0 0 0 0 0 0 0
To illustrate how this scheme works, superimpose it on a finite difference grid.
Lgrid(L)
The statement
A = delsq(L);
constructs matrix representation of the fivepoint discrete Laplacian for the grid L.
whos A
issym = isequal(A,A')
nnzA = nnz(A)
Name Size Bytes Class Attributes A 56x56 3156 double sparse issym = 1 nnzA = 244
So for this grid, the discrete Laplacian A is a 56by56 sparse double precision symmetric matrix with 244 nonzero elements. That's an average of
ratio = nnz(A)/size(A,1)
ratio = 4.3571
nonzero elements per row. Each interior grid point is connected to its four nearest neighbors. For example, point number 44 is connected to points 35, 43, 45, and 53. So the nonzero elements in column 44 are
A(:,44)
ans = (35,1) 1 (43,1) 1 (44,1) 4 (45,1) 1 (53,1) 1
There are 4's on the diagonal and 1's is positions on the offdiagonals corresponding to connections between neighbors. They can be seen in the spy plot which shows the location of the nonzeros.
spy(A)
Scale A by the square of the mesh size and then use eigs to find its five smallest eigenvalues.
h = 2/n
A = A/h^2;
eigvals = eigs(A,5,'sm')
h = 0.2000 eigvals = 30.3140 27.7177 19.0983 14.6708 9.6786
Run for decreasing mesh size, up to the memory capacity of my laptop.
type L_finite_diffs
% L_finite_diffs % Script for MathWorks Logo, Part Two % Finite Differences opts.tol = 1.0e12; eigvals = zeros(13,5); tic for n = 50:50:650 L = numgrid('L',n+1); h = 2/n; A = delsq(L)/h^2; e = eigs(A,5,'sm',opts); eigvals(n/50,:) = fliplr(e'); fprintf('%4.0f %10.6f %9.0f %10.0f %16.12f\n', ... n,h,size(A,1),nnz(A),e(5)) end toc save eigvals eigvals
type L_finite_diffs_results
n h size(A) nnz(A) lambda_1,h 50 0.040000 1776 8684 9.661331286788 100 0.020000 7301 36109 9.649547111146 150 0.013333 16576 82284 9.645736333124 200 0.010000 29601 147209 9.643932372382 250 0.008000 46376 230884 9.642903412412 300 0.006667 66901 333309 9.642247498956 350 0.005714 91176 454484 9.641797238078 400 0.005000 119201 594409 9.641471331746 450 0.004444 150976 753084 9.641225852818 500 0.004000 186501 930509 9.641035124947 550 0.003636 225776 1126684 9.640883203502 600 0.003333 268801 1341609 9.640759699387 650 0.003077 315576 1575284 9.640657572983 Elapsed time is 29.729023 seconds.
I know the Taylor series analysis of the discretization error involves terms with both $h^{4/3}$ and $h^2$, so let's use these as the basis for extrapolation. Use backslash to do a least squares fit.
load eigvals eigvals = eigvals(:,1); n = (50:50:650)'; h = 2./n; e = ones(size(h)); X = [e h.^(4/3) h.^2]; c = X(4:end,:)\eigvals(4:end); lambda1 = c(1); n = (50:10:650)'; h = 2./n; fit = c(1) + c(2)*h.^(4/3) + c(3)*h.^2; plot(n(1:5:end),eigvals,'o',n,fit,'',[0 700],[lambda1 lambda1],'k') set(gca,'xtick',100:100:600,'xticklabel',... num2str(sprintf('2/%3.0f\n',100:100:600))) xlabel('h')
Frankly, I was surprised when the extrapolation of the finite difference eigenvalues produced this result.
format long
lambda1
lambda1 = 9.639723753239963
I know from techniques that I will describe in later posts that this agrees with the exact eigenvalue of the continuous differential operator to seven decimal places. Before I wrote this blog I never tried this computation on computers with enough memory to handle meshes this fine and was never able to obtain this kind of accuracy from extrapolated finite difference eigenvalues.
Let's do a contour plot of the eigenfunction. The tricky part is using the grid numbering in L to insert the eigenvector back onto the grid.
close n = 100; h = 2/n; L = numgrid('L',n+1); A = delsq(L)/h^2; [V,E] = eigs(A,5,'sm'); L(L>0) = V(:,5); contourf(fliplr(L),12); axis square axis off
The plot shows off parula, the MATLAB default colormap with Handle Graphics II and Release 2014b. The colormap derives its name from the parula, which is a small warbler found throughout eastern North America that exhibits these colors. As our contour plot demonstrates, the color intensity varies smooth and uniformly as the contour levels increase. This is not true of jet, the old default. Steve Eddins has more about parula in his blog.
Get
the MATLAB code
Published with MATLAB® R2014b
Last week I showed you the new MATLAB colormap, parula. Parula has replaced jet as the default colormap in R2014b, which was released earlier this month.
This week I want to explain some of the motivations for replacing jet. Jet is an example of a rainbow colormap. A rainbow colormap is based on the order of colors in the spectrum of visible light. Some of you might know the mnemonic "Roy G. Biv" for remembering the order of colors: red, orange, yellow, green, blue, indigo, and violet.
Multiple methods exists for constructing rainbow colormaps, so they don't all look alike. Here is jet, the MATLAB version:
showColormap('jet','bar')
In my introductory post last week, I showed you a few visualization examples using jet, and I asked you some questions about them. Let's look at the questions again and figure out the answers.
Question 1: In the chart below, as you move from left to right along the line shown in yellow, how does the data change? Does it trend higher? Or lower? Or does it trend higher in some places and lower in others?
fluidJet hold on plot([100 160],[100 135],'y','LineWidth',3) hold off colormap(gca,'jet')
By my count there are six very distinct color changes as you move along the line from left to right. So what does that mean about the data?
I posed this question to someone who is very familiar with the jet colormap, and I was surprised by the response: "The data starts low, goes high, and then goes back low again."
"But what about all those other color stripes," I asked? "What do they mean?"
"That's just what jet does. I've learned to ignore all those stripes."
Oh! Well, it turns out that some of the color stripes are indeed mostly meaningless. But not all of them!
Here's what the data actually does along that path.
improfile(s.X,[100 160],[100 135])
Yes, the data starts low, goes high, and then goes back low. But if you ignore all the color stripes in between the blue and the red, then you miss the presence of the smallermagnitude peak on the left. Alternatively, if you believe that all the color stripes are meaningful, then the yellow stripes seem to falsely suggest distinct data regions to the left and the right of the main peak.
Now let's compare with parula.
fluidJet colormap(gca,parula) hold on plot([100 160],[100 135],'y','LineWidth',3) hold off
You can see that now there are only four distinct color changes as you follow the yellow line. These changes correspond directly to real features in the data. There are no extraneous color stripes that do not correspond to real data features.
Question 2: In the filled contour plot below, which regions are high and which regions are low?
filledContour15 colormap(gca,rgb2gray(jet))
In the plot above, I have done a crude simulation of what a jetbased visualization looks like if printed on a grayscale printer. Here's the full color version:
colormap(gca,jet)
If you are familiar with jet and know therefore that blue is low and red is high, then you can probably look at the fullcolor version and give a reasonable interpretation of it. When rendered as grayscale, however, all hope is lost. The red and blue colors of jet are just about equally dark.
Here are the fullcolor and grayscale results with parula.
filledContour15 colormap(gca,parula)
filledContour15 colormap(gca,rgb2gray(parula))
With parula, dark consistently means low, and bright consistently means high.
The last questions relate to the three plots below (A, B, and C) showing different horizontal oscillations.
Question 3: Which horizontal oscillation (A, B, or C) has the highest amplitude?
Question 4: Which horizontal oscillation (A, B, or C) is closest to a pure sinusoid?
Question 5: In comparing plots A and C, which one starts high and goes low, and which one starts low and goes high?
subplot(2,2,1) horizontalOscillation1 title('A') subplot(2,2,2) horizontalOscillation2 title('B') subplot(2,2,3) horizontalOscillation3 title('C') colormap(jet)
The answers are subjective and might depend on your monitor and ambient lighting. But here's what I see.
The data pattern in image A has significantly higher visual contrast than either image B or image C. That suggests that the data oscillation in A has the highest amplitude.
The data patterns in images B and C both exhibit vertical stripes within each period of the oscillation. That strongly suggests that horizontal oscillations shown in B and C have some constant or almostconstant regions and that they go up and down in something like a stairstep pattern. To my eye, image A looks the most like a pure sinusoid.
And finally, the pattern in image A starts with a bright stripe on the left, whereas the image B pattern starts with a dark stripe. That suggests that oscillation A starts high, while oscillation B starts low (or maybe vice versa).
You probably won't be surprised to find out that these are all trick questions. Actually, oscillations A, B, and C are all sinusoids with exactly the same amplitude and shape. They differ only in their constant offsets, which place them in different regions of the jet colormap. Note especially that oscillations A and C are actually in phase with each other; they both start high and go low. The differences we see are all visual artifacts caused by the jet colormap.
Compare the result using parula:
colormap(parula)
With these examples, I wanted to show several ways in which the jet colormap, and rainbow colormaps in general, can mislead the viewer about properties of the data.
It turns out that some people have been writing about these problems with rainbow colormaps for years. Here's a summary of the main criticisms:
For a much more detailed summary of the literature and online sources regarding rainbow colormaps, see the paper "Rainbow Color Map Critiques: An Overview and Annotated Bibliography" on mathworks.com.
Next time I'll starting getting into colormap construction principles.
Get
the MATLAB code
Published with MATLAB® R2014b
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:
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 commandline 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 HelloWorld 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. Rightclicking 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.
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 wikilike 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, rightclick 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 rightclicking 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.
]]>Coursera is a technology platform that kickstarted the current MOOCs boom. Even though there are more MOOCs players now, it still remains one of the leading companies in this space. But how are they doing these days for delivering higher education to the masses online?
Today's guest blogger, Toshi Takeuchi, would like to share an analysis using Courera's data.
I am a big fan of MOOCs and I benefited a lot from free online courses on Coursera, such as Stanford's Machine Learning course. Like many websites these days, Coursera offers its data through REST APIs. Coursera offers a number of APIs, but Catalog APIs are available without OAuth authentication. We can find out the details of courses offered by Coursera with these APIs.
We can try to answer questions like "how do STEM and nonSTEM courses break down among universities?"
JSON is a very common data format for REST APIs, and Coursera's APIs also returns results in JSON format. MATLAB now supports JSON out of the box in R2014b. You could always use JSON from within MATLAB by taking advantage of user contributed MATLAB programs on File Exchange, but builtin JSON support makes it easy for us to share scripts that use JSON, because we don't have to worry about dependencies.
Let's try the new feature using Coursera APIs. Calling a REST API is very simple with webread.
restApi='https://api.coursera.org/api/catalog.v1/courses'; params = 'sessions,universities,categories'; resp=webread(restApi,'includes',params,weboptions('Timeout',60));
webread returns the JSON response as a structure array. The data is further processed in a separate script processData.m  check out the details if interested.
We need to decide which categories represent STEM subjects. When there are multiple categories assigned to a given course, we treat it as a STEM course as long as one of them is included in STEM categories.
processData
STEM categories 'Computer Science: Theory' 'Economics & Finance' 'Medicine' 'Mathematics' 'Physical & Earth Sciences' 'Biology & Life Sciences' 'Computer Science: Systems & Security' 'Computer Science: Software Engineering' 'Engineering' 'Statistics and Data Analysis' 'Computer Science: Artificial Intelligence' 'Physics' 'Chemistry' 'Energy & Earth Sciences'
As a sanity check, let's plot the number of courses vs. number of sessions by university. A single course can be offered repeatedly in multiple sessions. Therefore you can determine the longevity or age of a given course by the count of sessions.
If it is a new course, or it was not repeated, then you only have one session per course. We can use this as the baseline, and check how universities scaled up their courses relative to this baseline.
R2014b comes with new MATLAB Graphics System, but you can still use the familiar commands for plotting.
% group by number of courses grouping = ones(height(universities),1)*2; grouping(universities.courses > 25) = 1; grouping(universities.courses <= 10) = 3; % plot figure gscatter(universities.courses,universities.sessions,grouping) h = refline(1,0); set(h,'Color','m','LineStyle',':') h = refline(2,0); set(h,'Color','m','LineStyle',':') h = refline(3,0); set(h,'Color','m','LineStyle',':') h = refline(6,0); set(h,'Color','m','LineStyle',':') xlabel('Number of Courses'); ylabel('Number of Sessions'); title('\fontsize{14} Courses by Sessions by University'); legend('Universities with 25+ courses','Universities with 10+ courses',... 'Universities with 110 courses','Ref line: 1 session per course',... 'Ref line: 2 sessions per course','Ref line: 3 sessions per course',... 'Ref line: 6 sessions per course','Location','NorthWest') % add university names for i = 1:height(universities) if universities.courses(i) > 10 && universities.sessions(i) > 20 text(universities.courses(i),universities.sessions(i),... universities.shortName{i},'FontSize',12) end end
You can see that Stanford, Penn (University of Pennsylvania), JHU (Johns Hopkins), and Duke are leading the pack. They are the early adopters, based on the number of sessions. It is interesting to see PKU (Peking University) leading international institutions. They offer a number of courses in Chinese. Coursera didn't start international partnership until recently, so it is quite remarkable the PKU has broadened their online content in relatively short time. More recent entrants are on the left with fewer courses and sessions.
Established players are trying to scale up by repeating the sessions. JHU seems to be particularly aggressive in terms of the number of courses they offer and how they are repeated as sessions.
Let's plot the number of courses by ratio of STEM courses by university. This will tell us which schools are making investments in online education content, and whether they focus on STEM or nonSTEM subjects. The size of the marker indicates the total number of sessions they are associated with, so it also gives us how long they have been involved in Coursera. Notice the parula colormap used in the colorbar, the new default colormap in R2014b.
% use sesion count for setting marker sizes markerSize = universities.sessions; % we need to scale the marker size markerSize = (markerSize  min(markerSize))/(max(markerSize)min(markerSize)); markerSize = markerSize * 1000; markerSize(markerSize == 0) = 1; % change the tick labels to reflect the original values barticks = num2cell(20:20:200); % create a scatter plot figure scatter(universities.courses,universities.stem_ratio,markerSize,markerSize,'fill') xlim([0 40]) h = colorbar('TickLabels',barticks); h.Label.String = '\fontsize{11}Number of Sessions'; title('\fontsize{14} Ratio of STEM courses by University on Coursera') xlabel('\fontsize{11}Number of Courses'); ylabel('\fontsize{11}Ratio of STEM Courses'); % add university names for i = 1:height(universities) if universities.stem_ratio(i) ~= 0 && universities.stem_ratio(i) ~= 1 && universities.courses(i) >= 5 text(universities.courses(i),universities.stem_ratio(i),universities.shortName{i},'FontSize',12) end end % add reference lines line([25 25],[0 1],'LineStyle',':') line([10 10],[0 1],'LineStyle',':') line([0 40],[0.5 0.5],'LineStyle',':')
Stanford is very heavy on STEM subjects, while others are more balanced. More recent entrants on the left have a wider variance in how STEM heavy their courses are. Perhaps rate of adaption is different among different academic disciplines?
We can plot the ratio of courses per category in order to see the relative representation of academic disciplines on Coursera. A course can belong to multiple categories, and in such cases a count is split equally across the included categories. Note that you can now rotate axis tick labels in R2014b.
% get the count of categories by university catByUniv = zeros(height(universities),height(categories)); for i = 1:length(T.categories) row = ismember(universities.id,T.universities(i)); col = ismember(categories.id,T.categories{i}); catByUniv(row,col) = catByUniv(row,col) + 1/length(T.categories{i}); end % segment the universities by number of courses catByTiers = [sum(catByUniv(grouping == 1,:));... sum(catByUniv(grouping == 2,:)); sum(catByUniv(grouping == 3,:))]; % get the ranking of categories by number of courses [~,ranking] = sort(sum(catByUniv(universities.courses > 25,:)),'descend'); % get the ratio of courses by category catByTiers = bsxfun(@rdivide,catByTiers,sum(catByTiers,2)); % plot a bar graph figure xticks = [{''};categories.name(ranking);{''}]; h = bar(catByTiers(:,ranking)'); xlim([0 26]); ax = gca; set(ax,'XTick',0:26);set(ax,'XTickLabel',xticks);set(ax,'XTickLabelRotation',270); title('\fontsize{14} Ratio of Courses Per Category') legend('Universites with 25+ courses','Universites with 10+ courses','Universites with 110 courses','Location','Best')
It looks like there was more STEM bias among the early adopters (universities with a lot of courses) but new entrants (universities with fewer courses) tend to have more nonSTEM courses. Categories like Social Sciences, Humanities, Business and Management, Education, Teacher Professioal Development, Music, Film and Audio are on the rise.
Why do we see this nonSTEM shift? There are a number of possible explanations.
There are questions we cannot answer with the data at hand. For example, is the shift driven by the convenience of supply side (universities) or by the demand for nonSTEM subjects by the public?
There are no strict prerequisites for Coursera courses, but the bar is still high for STEM courses. Therefore it is quite possible that the potential market size is larger for nonSTEM subjects.
You also saw how easy it is to use REST API with JSON response within R2014b, and got a quick look at some of the new features of updated MATLAB Graphics System. Download the new release, try those new features yourself and share what you find here!
Get
the MATLAB code
Published with MATLAB® R2014b
Differential Equation Editor
Try typing dee in MATLAB. This will launch an example model that looks like:
If you open one of the demo and doubleclick on the block, you will see a nice little user interface:
In this interface, you can type any equation you want, using the format of the Fcn block.
When I first saw this model and user interface, I thought "This thing looks really old, has it always been there, unnoticed?". So I looked into our source code history and found out this example ships with Simulink and has been untouched since the middle of the 90's!
How does that work?
I strongly encourage you to look under the mask and observe how this block is implemented. It is a pretty good example of how to create models programmatically. Once you click the rebuild button in the user interface, functions like delete_block, add_block and add_line draw the differential equation in a systematic manner.
For any set of equation, the resulting model will look like this, where one Fcn block is used for each state equation and for each output equation:
To help you visualize it better, I formatted the above model to make it easier to read:
Now it's your turn
Give this cute little example a try and let us know what you think by leaving a comment here.
]]>Today, David Garrison, our guest blogger, will continue his series on the new graphics system in R2014b.
Here is Part 2 of the series.
In Part 1 of this series, I provided an introduction to the new MATLAB graphics system in R2014b. I described a number of new features and alluded to one of the big changes in R2014b  graphics functions return MATLAB objects, not numeric handles.
When I use MATLAB, I usually don't think about the internal graphics system. Usually, I'm just trying to visualize some data. I might call one of the MATLAB charting functions, like plot or bar to understand something important about my data. I also create user interfaces using functions like uicontrol or uipanel. These are functions that are part of the MATLAB graphics system.
Let's suppose that I am creating a plot and I want the line in the plot to be red. I can do that in one of two ways. I can call the plot function with some extra arguments like this.
x = 0:0.1:10; y = sin(x); plot(x,y,'Color','red')
I can also use the output argument of the plot command to modify my plot after I've created it. Here I created the plot first, then used the set command to change the color of the line.
p = plot(x,y) ; set(p, 'Color', 'red')
Prior to R2014b, a call to a graphics function would return a number like this
p = plot(x,y)
p = 174.346
The value of the variable p was a special kind of number called a handle. The handle was a reference to a graphics object. The object itself was not available to the user but you could use the handle to retrieve or change properties of the object using the set and get commands. In the example above, we used the handle of the line object to set its color. In fact, when this capability first became available in MATLAB, this approach was so new and useful it was given a special name  Handle Graphics.
All graphics handles were just numeric values. Those values were integers for figures and fractional values for everything else. You could use a handle anywhere you could use a number. I've seen MATLAB code where people stored their handles with other data in a double array or used handles in functions where a number is expected as the argument (e.g. math functions).
As you probably know, graphics objects are stored as a tree structure with parents and children. There is a special object at the top of the tree called the graphics root. Prior to R2014b, the handle to the graphics root was always 0. For example to get the list of all open figures, you could type:
myFigures = get(0, 'Children');
In R2014b, graphics functions do not return numeric handles any more. They now return MATLAB objects. Now when you call a graphics function with an output argument, you see something like this:
p = plot(x,y)
p = Line with properties: Color: [0 0.4470 0.7410] LineStyle: '' LineWidth: 0.5000 Marker: 'none' MarkerSize: 6 MarkerFaceColor: 'none' XData: [1x101 double] YData: [1x101 double] ZData: [1x0 double] Use GET to show all properties
The plot function now returns a Line object. That line object has properties. When MATLAB displays the line object in the command window, it shows the object's most common properties. The all properties text is a hyperlink that you can click to see all the properties of the object.
If you use the whos function to see more information about p, you see something like this:
whos p
Name Size Bytes Class Attributes p 1x1 112 matlab.graphics.chart.primitive.Line
The Class column in the whos output indicates what kind of object if is. What you see here is the full class name of the object. Don't worry too much about the details of the class name. Most of you will never need to use it.
In R2014b, use the groot function to refer to the graphics root.
groot
ans = Graphics Root with properties: CurrentFigure: [1x1 Figure] ScreenPixelsPerInch: 96 ScreenSize: [1 1 1920 1200] MonitorPositions: [1 1 1920 1200] Units: 'pixels' Use GET to show all properties
One of the benefits of having graphic objects is that you can now get and set their values using dot notation. Dot notation has the form object.Property. It is the same syntax that you use when you want to refer to the field of a structure. For example, in R2014b, you can get the color of the line above like this:
p.Color
ans = 0 0.4470 0.7410
Similarly, you can set the line width of the line using dot notation.
p.LineWidth = 2;
Dot notation is also useful if you want to set properties of two objects to be the same. Prior to R2014b, you might write code that looks like this:
set(p2, 'Color', get(p1, 'Color'))
With dot notation, you can do the same thing with
p2.Color = p1.Color;
With dot notation you can also use tab completion when you are trying to access a graphics object property.
You can still use the set and get functions but I find dot notation a much more convenient way to refer to object properties. One note of caution, however. When you use dot notation you must use the correct capitalization of the property name. For example,
p.LineWidth
will return the value of the line width but
p.linewidth
will cause an error.
As I mentioned above, you can still use the set and get functions and there are cases where set and get are useful.
You can use the get command if you want to know all the properties of a graphics object. This is the same list of properties you see if you click the all properties link in the example above.
get(p)
AlignVertexCenters: 'off' Annotation: [1x1 matlab.graphics.eventdata.Annotation] BeingDeleted: 'off' BusyAction: 'queue' ButtonDownFcn: '' Children: [] Clipping: 'on' Color: [0 0.4470 0.7410] CreateFcn: '' DeleteFcn: '' DisplayName: '' HandleVisibility: 'on' HitTest: 'on' Interruptible: 'on' LineStyle: '' LineWidth: 2 Marker: 'none' MarkerEdgeColor: 'auto' MarkerFaceColor: 'none' MarkerSize: 6 Parent: [1x1 Axes] PickableParts: 'visible' Selected: 'off' SelectionHighlight: 'on' Tag: '' Type: 'line' UIContextMenu: [] UserData: [] Visible: 'on' XData: [1x101 double] XDataMode: 'manual' XDataSource: '' YData: [1x101 double] YDataSource: '' ZData: [1x0 double] ZDataSource: ''
The set function is useful if you want to know what options are available for a given property. In this case you can call set with the property name but no value for the property.
set(p,'LineStyle')
'' '' ':' '.' 'none'
Another use of set occurs when you need to reference a property from an array of graphics objects. Suppose you have an array of line objects created by the plot command:
x = 0:0.1:10; y1 = sin(x); y2 = cos(x); lineArray = plot(x,y1,x,y2)
lineArray = 2x1 Line array: Line Line
You can use the set function to change the colors of all the lines in the array in a single command.
set(lineArray, 'Color', 'blue')
What do you think about the change to graphics objects in R2014b? Have you tried using dot notation for getting and setting graphics object properties? We'd love to hear your thoughts here.
There are some other changes in R2014b that will impact some existing MATLAB code. These changes mostly impact advanced graphics users and people who build more complicated user interfaces. Next time, I will describe the important changes and what steps you should take to make your code compatible with R2014b. I'll also give you some guidance on how to write code that works in multiple releases.
Get
the MATLAB code
Published with MATLAB® R2014b
I believe it was almost four years ago that we started kicking around the idea of changing the default colormap in MATLAB. Now, with the major update of the MATLAB graphics system in R2014b, the colormap change has finally happened. Today I'd like to introduce you to parula, the new default MATLAB colormap:
showColormap('parula','bar')
(Note: I will make showColormap and other functions used below available on MATLAB Central File Exchange soon.)
I'm going spend the next several weeks writing about this change. In particular, I plan to discuss why we made this change and how we settled on parula as the new default.
To get started, I want to show you some visualizations using jet, the previous colormap, and ask you some questions about them.
Question 1: In the chart below, as you move from left to right along the line shown in yellow, how does the data change? Does it trend higher? Or lower? Or does it trend higher in some places and lower in others?
fluidJet hold on plot([100 160],[100 135],'y','LineWidth',3) hold off colormap(gca,'jet')
Question 2: In the filled contour plot below, which regions are high and which regions are low?
filledContour15 colormap(gca,rgb2gray(jet))
The next questions relate to the three plots below (A, B, and C) showing different horizontal oscillations.
Question 3: Which horizontal oscillation (A, B, or C) has the highest amplitude?
Question 4: Which horizontal oscillation (A, B, or C) is closest to a pure sinusoid?
Question 5: In comparing plots A and C, which one starts high and goes low, and which one starts low and goes high?
subplot(2,2,1) horizontalOscillation1 title('A') subplot(2,2,2) horizontalOscillation2 title('B') subplot(2,2,3) horizontalOscillation3 title('C') colormap(jet)
Next time I'll answer the questions above as a way to launch into consideration of the strengths and weaknesses of jet, the previous default colormap. Then, over the next few weeks, I'll explore issues in using color for data visualization, colormap construction principles, use of the L*a*b* color space, and quantitative and qualitative comparisons of parula and jet. Toward the end, I'll even discuss how the unusual name for the new colormap came about. I've created a new blog category (colormap) to gather the posts together in a series.
If you want a preview of some of the issues I'll be discussing, take a look at the technical paper "Rainbow Color Map Critiques: An Overview and Annotated Bibliography." It was published on mathworks.com just last week.
Get
the MATLAB code
Published with MATLAB® R2014b
MathWorks is the only company in the world whose logo satisfies a partial differential equation. Why is the region for this equation shaped like a capital letter L?
The wave equation describes how a disturbance travels through matter. If the units are chosen so that the wave propagation speed is equal to one, the amplitude of a wave satisfies
$$ {{\partial^2 u} \over {\partial t^2}} = \triangle u $$
The $\triangle$ denotes Laplace's operator
$$ \triangle = {\partial^2 \over {\partial x^2}} + {\partial^2 \over {\partial y^2}} $$
Geometry plays a crucial role here. Initial values of the amplitude and velocity of the wave are prescribed on a certain region. Values of the amplitude or its normal derivative are also prescribed on the boundary of the region. If the wave vanishes outside the region, these boundary values are zero.
Separating out periodic time behavior leads to solutions of the form
$$ u(x,y,t) = \cos{(\sqrt{\lambda}\,t)} v(x,y) $$
The functions $v(x,y)$ also depend upon $\lambda$ and the region. They satisfy the differential equation
$$ \triangle v + \lambda v = 0 $$
and are zero on the boundary of the region. The quantities $\lambda$ that lead to nonzero solutions are the eigenvalues, and the corresponding functions $v(x,y)$ are the eigenfunctions or modes. They are determined by the physical properties of the medium and the geometry of the region. The square roots of the eigenvalues are resonant frequencies. A periodic external driving force at one of these frequencies generates an unboundedly strong response in the medium.
Any solution of the wave equation can be expressed as a linear combination of these eigenfunctions. The coefficients in the linear combination are obtained from the initial conditions.
In one dimension, the eigenvalues and eigenfunctions are easily determined. The simplest example is a violin string, held fixed at the ends of an interval. The eigenfunctions are trig functions.
$$ v_k(x) = \sin{(k x)} $$
If the length of the interval is $\pi$, the eigenvalues are determined by the boundary condition, $v_k(k \pi) = 0$. Hence, $k$ must be an integer and
$$ \lambda_k = k^2 $$
If the initial condition is expanded in a Fourier sine series,
$$ u(x,0) = \sum_k a_k \sin{(k x)} $$
(And the initial velocity is zero), then the solution to the wave equation is
$$ u(x,t) = \sum_k a_k \cos{(\sqrt{\lambda_k}\,t)} v_k(x) $$
Here are graphs of the first nine eigenfunctions in one dimension. The corresponding eigenvalues are the squares of integers.
eigenvals = (1:9).^2
plot_modes('1d')
eigenvals = 1 4 9 16 25 36 49 64 81
The simplest region in two dimensions is a square. The eigenfunctions are again trig functions.
$$ v_{k,j}(x,y) = \sin{(k x)}\,\sin{(j y)} $$
If the sides have length $\pi$, the boundary conditions imply that $k$ and $j$ must be integers. Here are the first nine eigenvalues and eigenfunctions.
[k,j] = meshgrid(1:3); e = k.^2+j.^2; eigenvals = e(:)'
plot_modes('square')
eigenvals = 2 5 10 5 8 13 10 13 18
If the region is a circular disc, we switch to polar coordinates, $r$ and $\theta$. Trig functions are replaced by Bessel functions. The eigenfunctions become
$$ v_{k,j}(r,\theta) = B_j(\mu_k r) \,\sin{(j \theta)} $$
where $B_j$ is the $j$ th order Bessel function and $\mu_k = \sqrt{\lambda_k}$. To find the eigenvalues we need to have the eigenfunctions vanish on the boundary of the disc. If the radius is one, we require
$$ B_j(\mu_k) = 0 $$
In other words, we need to compute zeros of Bessel functions. Here are the first nine eigenvalues and eigenfunctions of the circular disc. The violin string has become a tambourine.
eigenvals = [bjzeros(0,3) bjzeros(1,3) bjzeros(2,3)].^2
plot_modes('disc')
eigenvals = Columns 1 through 7 5.7832 30.4713 74.8870 14.6820 49.2185 103.4995 26.3746 Columns 8 through 9 70.8500 135.0207
Replace the full circular disc by a threequarter circular sector. The angle at the origin is $270^\circ$ or $\frac{3}{2}\pi$ radians. We can make our eigenfunctions adapt to this angle. Take
$$ v_{k,j}(r,\theta) = B_{\alpha_j}(\mu_k r) \,\sin{(\alpha_j \theta)} $$
with
$$ \alpha_j = \frac{2}{3} j $$
and fractional order Bessel functions. The eigenfunctions satisfy their differential equation and also satisfy the boundary conditions on both sides of the angle.
$$ v_{k,j}(r,\theta) = 0 $$
at $\theta = 0$ and at $\theta = \frac{3}{2}\pi$.
By finding the zeros of the Bessel functions we can also have the eigenfunctions satisfy the boundary conditions on the outer circular portion of the boundary. Here are the first nine eigenvalues and eigenfunctions of the threequarter circular sector.
eigenvals = [bjzeros(2/3,3) bjzeros(4/3,3) bjzeros(6/3,3)].^2
plot_modes('sector')
eigenvals = Columns 1 through 7 11.3947 42.6442 93.6362 18.2785 56.1131 113.6860 26.3746 Columns 8 through 9 70.8500 135.0207
These eigenfunctions have another important property. Most of them are singular; the derivatives of the fractional order Bessel functions are unbounded at the origin. You can see the black concentration of grid lines in the plots. If you tried to make a tambourine with this sector shape, it would rip at the sharp corner. This singular behavior is needed to model the solution to the wave equation on this region.
For all the regions we have discussed so far it is possible to express the eigenvalues as zeros of analytic functions. For the interval and the square, the eigenvalues are integers, which are the zeros of $\sin{\pi x}$. For the circular disc and sector, the eigenvalues are zeros of Bessel functions. Once we have the eigenvalues, it is easy to compute the eigenfunctions using sines and Bessel functions.
So, an Lshaped region formed from three unit squares is interesting for at least two reasons. It is the simplest geometry for which solutions to the wave equation cannot be expressed analytically; numerical computation is necessary. Furthermore, the 270 degree nonconvex corner causes a singularity in the solution. Mathematically, the gradient of the first eigenfunction is unbounded near the corner. Physically, a membrane stretched over such a region would rip at the corner. This singularity limits the accuracy of finite difference methods with uniform grids.
I used the Lshaped region as the primary example in my doctoral thesis fifty years ago. MathWorks has adopted a modified surface plot of the first eigenfunction as the company logo. I am going to devote a series of blog posts to the L.
Here are the first nine eigenvalues and eigenfunctions, computed by a function from Numerical Computing with MATLAB, which I will discuss in a future posting. Compare these eigenfunctions with the ones from the circular sector, which has the same reentrant corner and resulting singularity.
for k = 1:9 [~,eigenvals(k)] = membranetx(k); end eigenvals plot_modes('L')
eigenvals = Columns 1 through 7 9.6397 15.1973 19.7392 29.5215 31.9126 41.4745 44.9485 Columns 8 through 9 49.3480 49.3480
Simple model problems involving waves on an Lshaped region include an Lshaped membrane, or Lshaped tambourine, and a beach towel blowing in the wind, constrained by a picnic basket on one fourth of the towel.
A more practical example involves ridged microwave waveguides. One such device is a waveguidetocoax adapter. The active region is the channel with the Hshaped cross section visible at the end of the adapter. The ridges increase the bandwidth of the guide at the expense of higher attenuation and lower powerhandling capability. Symmetry of the H about the dotted lines shown in the contour plot of the electric field implies that only one quarter of the domain needs to be considered and that the resulting geometry is our Lshaped region. The boundary conditions are different than our membrane problem, but the differential equation and the solution techniques are the same.
The photo is courtesy of Advanced Technical Materials, Inc. See their website, <http://atmmicrowave.com/>, for lots of devices like this.
waveguide
Get
the MATLAB code
Published with MATLAB® R2014b
We are very excited to announce that MATLAB r2014b introduces a new graphics system. Graphics developer, Mike Garrity, is launching a graphics specific blog. Here is the landing page for this release of MATLAB and for the new graphics system in general.]]>
I wanted to show you a glimpse of some of the new math functionality available in R2014b.
Recently on the MATLAB newsgroup, Christoph asked this question:
I have a vector A shown below, which has 6 elements. the elements are already sorted in descending order. now i want to create vector C by deleting elements from A, starting with element a1, until the sum of the vector equals or is smaller the value B
A= 26 23 20 19 15 14
B=70
So, the output should be
C= 20 19 15 14
Any idea how to do this?
I would do something like this. Get the numbers in increasing order and perform a cumulative sum.
Aflip = flipud(A); Asums = cumsum(Aflip);
Find the first sum that is GREATER than b. The next element is where you want to start. startIndex = find(Asums > b, 1, 'first');
But things are reversed now. So we need to find the starting index into the original array. Aind = length(A)  startIndex + 1;
Now delete the leading values you don’t want.
Afiltered = A((Aind+1):end);
In R2014b, the function cumsum has a new option, direction. (So do some other functions in the "cum" category.) You can specify to do the calculation in the default or 'forward' direction, or in 'reverse'. Here's what the new solution looks like.
A = [ 26 23 20 19 15 14]; B = 70; tmpA = cumsum(A,'reverse') ind = find(tmpA > B, 1,'last') C = A; C(1:ind) = []
tmpA = 117 91 68 48 29 14 ind = 2 C = 20 19 15 14
Have you had a chance to explore some of the math functionality in R2014b? Which features did you find helpful? Let me know here.
Get
the MATLAB code
Published with MATLAB® R2014b
Fast Restart
If you are dealing with large models, you are definitely aware that every time you click play to simulate a model, a good amount of time can be spent in initialization. In R2014b, if you know you will simulate your model multiple times without making structural changes, you can click the Fast Restart button:
At the end of the simulation, the model will remain in an Initialized state, and the toolbar will indicate this, telling you: "I am ready for fast restart".
At this point you can change the value of any tunable parameter and click play. The simulation will start instantly.
Quick Insert
Those who prefer typing over using the mouse will love the Quick Insert feature in R2014b. Type whatever you need in the canvas and the blocks magically appear! Once the block is there, a pop up offers to set the value of the most common parameter for this block:
Model Templates
You can now create models from templates:
This will launch the Simulink Template Gallery:
Of course, you can add your own templates to the gallery using the File menu of any model:
Now it's your turn
As you can see, R2014b contains many new features designed especially to speed up your workflow.
What is your favorite feature in R2014b? Let us know by leaving a comment here
]]>m = mobiledev;
m.AccelerationSensorEnabled = 1;Step 2: Send Data to MATLAB
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: '02Oct2014 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)Step 3: Stop Acquiring Data and Retrieve Logs
m.Logging = 0;To retrieve the data, use the accellog variable:
[a, t] = accellog(m);Step 4: Plot Raw Sensor Data
plot(t, a); legend('X', 'Y', 'Z'); xlabel('Relative time (s)'); ylabel('Acceleration (m/s^2)');
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)');
% Accounting for gravity. magNoG = mag  mean(mag); % Plot magnitude. plot(t, magNoG); xlabel('Time (s)'); ylabel('Acceleration (m/s^2)');
% 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;
m.AccelerationSensorEnabled = 0;
clear m;
Published with MATLAB® R2014b
]]>Mike is a computer graphics guru, well informed on the state of the art and kneedeep in the operational knowledge of How It All Works. He helped develop many of the features he’ll be describing, and he is a thoroughly entertaining writer. Beyond all this goodness, MATLAB graphics have taken a tremendous leap forward with the latest release, and Mike’s blog will give you an ideal perch to learn all about it. I recommend heading over there and seeing what he has to say.
]]>This is the third in a series of posts on the finite Fourier transform. The Fourier matrix produces an interesting graphic and has a surprising eigenvalue distribution.
The Fourier matrix of order $n$ is the $n$ by $n$ complex Vandermonde matrix $F$ whose elements $f_{k,j}$ are powers of the $n$ th root of unity
$$ \omega = e^{2 \pi i/n} $$
$$ f_{k,j} = \omega^{k j} $$
The matrix can be generated with the MATLAB statements
k = (0:n1)'; j = (0:n1); F = exp(2*pi*i*k*j/n);
Or, by taking the FFT of the identity matrix
F = fft(eye(n))
The statement
plot(F)
connects points in the complex plane whose coordinates are the real and imaginary parts of the elements of the columns of F, thereby generating a subgraph of the graph on n points. If n is prime, connecting the elements of all columns generates the complete graph on n points. If n is not prime, the sparsity of the graph of all columns is related to the computational complexity, and hence the speed, of the fast FFT algorithm. The graphs for n = 8, 9, 10, and 11 are generating and plotted by the following code.
for n = 8:11 subplot(2,2,n7) F = fft(eye(n)); plot(F) axis square axis off title(['n = ' int2str(n)]) end
Because n = 11 is prime, the corresponding graph shows all possible connections. But the other three values of $n$ are not prime. Some of the links in their graphs are missing, indicating that the FFT of a vector with that many points can be computed more quickly.
The remainder of this post examines F12 in more detail. Here is its graph.
clf n = 12; F = fft(eye(n)); plot(F) axis square axis off title(['n = ' int2str(n)])
The program fftmatrix, available here, or included with the NCM App, allows you to investigate these graphs. fftmatrix(n) plots all the columns of the Fourier matrix of order n. fftmatrix(n,j) plots just one column.
Let's plot the individual columns of F12. The first column of F12 is all ones, so its plot is just a single point.
for j = 1:12 p = mod(j1,4)+1; subplot(2,2,p); fftmatrix_mod(12,j1) title(['j = ' int2str(j)]) if p == 4, snapnow, end end
To see typical behavior, look at the third subplot, the red graph labeled j = 3, generated by the third column
F(:,3)
ans = 1.0000 + 0.0000i 0.5000  0.8660i 0.5000  0.8660i 1.0000 + 0.0000i 0.5000 + 0.8660i 0.5000 + 0.8660i 1.0000 + 0.0000i 0.5000  0.8660i 0.5000  0.8660i 1.0000 + 0.0000i 0.5000 + 0.8660i 0.5000 + 0.8660i
These are the powers of $\omega^2$. Because 2 divides 12 evenly these powers hit all the evennumbered points around the circle twice and miss all the oddnumbered points. Now look at the cyan graph labeled j = 11. These are the powers of $\omega^{10}$, which are the complex conjugates of the powers of $\omega^2$. So the two graphs lie on top of each other.
Six of the twelve columns in the graph of F12 connect only a subset of the nodes and ten of the columns lie on top of their complex conjugate columns. As a result, when all of the columns are combined to form the full graph, it is sparse. This sparsity, in turn, makes it possible to construct a fast finite Fourier transform algorithm for n = 12.
I've always been curious about the eigenvalues and vectors of the Fourier matrices. In 1979, three friends of mine at the University of New Mexico, Gus Efroymson, Art Steger, and Stan Steinberg, posed the question in the SIAM Review problem section. They didn't know that Jim McClellan and Tom Parks had actually solved their problem seven years earlier, when Jim was a grad student working under Tom at Rice. I didn't learn about the McClellanParks paper until fairly recently.
The Wikipedia page on the discrete Fourier transform says the facts about the eigenvalues are well known, but that the eigenvectors are the subject of ongoing research.
Let's rescale $F$ so that its columns have unit length.
F = F/sqrt(n);
This makes $F' \ F = I$
round(F'*F)
ans = 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1
Now $F$ is unitary, the complex generalization of orthogonal. This implies that all of its eigenvalues lie on the unit circle in the complex plane.
Furthermore, it turns out that $F^4 = I$.
round(F^4)
ans = 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1
This implies that any eigenvalue $\lambda$ must satisfy $\lambda^4$ = 1. Consequently the only possible eigenvalues are 1, 1, i, and i. You might guess that guess that if $n$ is divisible by 4, the eigenvalues would be equally distributed among these four values. But, surprisingly, to me at least, that doesn't happen.
lambda = eig(F)
lambda = 1.0000 + 0.0000i 1.0000 + 0.0000i 0.0000 + 1.0000i 1.0000 + 0.0000i 1.0000  0.0000i 0.0000 + 1.0000i 0.0000  1.0000i 1.0000 + 0.0000i 0.0000  1.0000i 1.0000  0.0000i 1.0000  0.0000i 0.0000  1.0000i
It's hard to pick through this unsorted output, but there are four +1's, three 1's, three i's, and only two +i's.
Here is a tricky piece of code that uses angle and the counting feature of sparse indexing to count the number of each of the four possible eigenvalues.
type eigfftmat
function c = eigfftmat(n) % EIGFFTMAT Count eigenvalues of the Fourier matrix. % c = eigfftmat(n) is a 4vector with counts for +1, 1, i, +i. % Compute the eigenvalues. e = eig(fft(eye(n))); % Sparse repeated indexing acts as a counter. c = full(sparse(mod(round(angle(e')/(pi/2)),4)+1,1,1))'; c([3 2]) = c([2 3]); end
When we run this code for a sequence of values of n, we see the pattern predicted by the McClellanParks analysis. The number of each of the four possible eigenvalues depends upon floor(n/4) and mod(n,4).
disp(' n +1 1 i +i') for n = 4:20 disp([n eigfftmat(n)]) end
n +1 1 i +i 4 2 1 1 5 2 1 1 1 6 2 2 1 1 7 2 2 2 1 8 3 2 2 1 9 3 2 2 2 10 3 3 2 2 11 3 3 3 2 12 4 3 3 2 13 4 3 3 3 14 4 4 3 3 15 4 4 4 3 16 5 4 4 3 17 5 4 4 4 18 5 5 4 4 19 5 5 5 4 20 6 5 5 4
The proofs in the McClellan and Parks paper involve the eigenvectors and are quite complicated. It turns out that this MATLAB expression
floor((n+[4 2 1 1])/4)
generates a 4vector of the multiplicities of the +1, 1, i, and +i eigenvalues for any given value of n.
J. H. McClellan and T. W. Parks, "Eigenvalues and eigenvectors of the discrete Fourier transformation", IEEE Trans. Audio Electroacoust 20, 6674, <http://dx.doi.org/10.1109/TAU.1972.1162342>, 1972.
G. Efroymson, A. Steger, and S.Steinberg, "A Matrix Eigenvalue Problem", SIAM Review, 21, 139139, <http://dx.doi.org/10.1137/1021013>, 1979.
Wikipedia, "Discrete Fourier Transform", <http://en.wikipedia.org/wiki/Discrete_Fourier_transform>, retrieved 09/03/2014.
Get
the MATLAB code
Published with MATLAB® R2014a
Let’s see how this works using a simple example!
The Problem
Let's say we have a setup where a DC motor is used to move a load between two positions. The position is measured using an encoder, and a PID is used to control the angle of the DC motor. The model looks like:
And when it simulates, the load tracks the desired position in that manner:
How could we detect if something breaks in the motor, or if the dynamics of the load change?
Online Identification
When starting the system we will spend a few seconds to identify its dynamics online (while the system is running). Once we have a model estimated, we will compare the encoder reading with the output of the model. If the error between the real system and the model becomes too large, it's a sign something undesired is happening... and we want to detect that!
I wrapped the control loop shown above in a subsystem and connected to it the Recursive Polynomial Model Estimator block from the System Identification Toolbox library.
Using the enable input port, we enable estimation for the first few seconds (only a few periods of the motor motion are necessary to get a good model). Once the estimation is completed, we monitor the estimation error, using very simple logic blocks to make a signal pass from zero to one as soon as the error level becomes larger than an acceptable level.
The Results
To test the logic in simulation, I faked the failure by doubling the friction in my DC motor model. As you can see, when the system dynamics change around 10 seconds, the error flag instantly signals the problem:
What's next?
The online estimation block also supports code generation. To see that in action, I recommend watching the video Online Fault Detection for a DC Motor by my colleague Karthik Srinivasan.
In this video, Karthik uses a model similar to the one I describe above, except that instead of a simulated DC Motor, he used driver blocks from the Arduino Simulink Suppport Package, to detect a failure using a real DC motor.
Now it's your turn
As you can imagine, fault detection is only one application of online parameters estimation. For example the identified model could be used in any sort of adaptive control.
Let us know what you think of the new online estimation features of the System Identification Toolbox by leaving a comment here.
]]>Marta is a senior application engineer at MathWorks UK, helping MATLAB users to design, develop and deploy productionquality code. Toshi is a Senior Marketing Manager at MathWorks.
by Marta Wilczkowiak and Toshi Takeuchi
In the previous post, you learned how we built an RFIDbased people tracking system from scratch in three months, a “technology stack” that spanned multiple hardware and software disciplines. We ended the post with how we managed to survive the moment of the truth: the first integration test, live at a conference with hundreds of delegates.
Now we had raw data to analyze to improve our system.
The Twist – Making Sense of RFID Data
Once we examined the real word data collected during the event more closely, we faced a new dilemma – we had originally assumed we could use signal strength to estimate distances. However, in reality there was so much fluctuation in signal strength that you could not make a confident estimate this way. This, incidentally, is why RFIDbased localization is a hot research topic today.
We were determined to work on an engineering solution, but it would take more time and effort. In the meantime, was there a way to get meaningful insights from the messy data at hand?
Web Analytics Meets Engineering
As we discussed the conference data internally at MathWorks, this issue came to the attention of Toshi Takeuchi, a USbased analyst. He immediately saw a similarity with the familiar problems he faced in analyzing massive volume of page views and clicks – always messy. What would happen if you applied web analytics techniques to this engineering problem? We embraced this suggestion and immediately shared our code and data with Toshi. Our collaboration was possible because we all spoke the same language – MATLAB.
Use of Proxy Metric
Nonphysical measurements are often used in web analytics. For example, the recommendations you see on online shopping or movie streaming websites are computed using distance measured in terms of similarity or affinity. Was there anything in the dataset, besides the signal strength, that we could turn into a proxy metric? The new approach Toshi came up with was to recast connection frequency as this virtual distance. If two tracking devices are connecting to one another more frequently than others in a given time interval, then we could think of them being “close” in the virtual space defined by connection frequency. This metric seemed far more reliable than signal strength. Our original goal was to understand the social interactions of the delegates, and physical distance metric was just a means to that end. If a virtual metric did the same job, then it seemed a good substitute.
The Result
With the support from Toshi, we were able to reconstruct the estimated movements of the delegates using the connection frequencies with stationary anchors such as demo stations (orange circles) and catering (green circles). Note that no personally identifiable information was collected or used for this analysis.
In the plot above, you see such an estimated movement for Delegate #120. This person appeared to have approached a nearby catering station first, then spent quite some time around BLOODHOUND Super Sonic Car in the middle, but in the afternoon stayed “closer” to software related demos than hardwarebased demos as measured in this “virtual” space.
We can use this information to map the frequencies of connections between delegates and demo stations, and visualize it as a heat map. The vertical bands with lighter red stripes indicates which demos had more connections than others, and horizontal stripes show to which demo stations each delegate connected frequently. You can see that people who connected with “RT” (a realtime testing exhibit) and “HW” (a hardware testing exhibit) had very little overlap with the people who connected with “Robo” (a robotics exhibit).
The Insight – the Three Factors of Delegate Behavior
Ultimately, the data is only as good as the insights you get from it. When we computed the estimated movements of the delegates with a technique called Principal Component Analysis, we could summarize the complex data by a mere three factors.
Furthermore, “Web” people and “Robo” people tend to share similar connection patterns as compared to “Conventional” people – what is common to both? They are both what we might call “convergence” applications where software and hardware are being integrated in an unconventional way.
While this might not be a novel insight, it was still very useful for us to know that the data backed up what we had been intuitively sensing all along. It was also interesting to see which delegates showed which particular inclinations on an individual basis (each dot in the plot identify a specific delegate).
Back to Reality – the Next Step
It is tempting to speculate if the virtual proximity on this plot correlates to actual shared interest among those delegates. If so, it is also interesting that two of the main factors turned out to be the web and robotics, which are examples of convergenceoriented applications that cross the traditional virtual/physical divide.
If there was an exit survey that asked the participants which demos they found interesting, we could use predictive analytics techniques to validate how well this virtual metric predicts the actual interest. Alternatively, we could use additional stationary anchors to determine which sessions delegates attended. If we placed more anchors on the floors in a more evenly spaced grid, including the areas populated by third party exhibitors, it would give more data points and more accurate measurements.
If we validated this correlation, it could open up opportunities for new ways for our delegates to find likeminded people.
The Journey Continues
Our journey is not over yet – we still have challenges with the system. But our project proved how such an interdisciplinary system could be built in three months using MATLAB, Simulink, and lots of team work.
]]>The story is broken into two separate blog posts. This week’s installment is told by Marta Wilczkowiak. Marta is a senior application engineer at MathWorks UK, helping MATLAB users to design, develop and deploy productionquality code.
by Marta Wilczkowiak
Technology can be the kiss of death to a conversation. But can we instead, through the thoughtful convergence of hardware, software, and realtime streaming data, actually encourage social interactions?
The barrier that existed between physical and virtual in engineering is disappearing in the world where internet companies buy robotics companies, fly drones, test selfdriving cars and wearable devices, while the industrial players are moving into sensor networks and Internet of Things. While those new approaches are promising, such interdisciplinary engineering collaboration is very difficult to achieve.
This is a story about our attempt to build such a convergent system in just three months, through crossdisciplinary collaboration. This bluesky project is called Expo Conversations.
EXPO Conversations: The Idea
MATLAB EXPO is the pinnacle of MathWorks UK events. It is a unique opportunity for the MATLAB community to network with peers across a range of industries and academia. Planning starts early in the year and involves engineers from all technical groups at MathWorks. During one of those meetings, we set ourselves a challenge: to create a system bringing together all the new technologies we tell our customers about, including low cost hardware, data analytics and cloud computing. At the same time, we saw an opportunity to understand how ideas and influences flow throughout EXPO.
Our idea was to use low cost hardware to detect the interactions between the EXPO participants, and provide a web interface displaying the live analysis of the social activity. To encourage social interactions we would run a contest and give a prize to the most prolific networker at the event.
Figure 1: Simulated social interactions in the EXPO Exhibition Hall
This conversation took place in June. EXPO was scheduled for October 8th so we had little time to develop something from nothing!
To summarize, the objective was to:
July: First tag
14 weeks to go
July started with discussion about the best approach to detect the attendee interactions: Cameras? Apps on mobile phones? RF tags? We wanted participants to remain anonymous so we opted for RF tags, with attendees optingin simply by picking up a tag on their way into the exhibition hall. Conor, our 2nd year summer intern, with support from Mark (embedded systems) and Graham (communications) developed the first tags using Arduino Uno microcontrollers (ATmega328P) and RF12 433MHz transceivers. They showed that two early tag prototypes can communicate with each other and the received signal strength is somewhat related to the distance between tags. Characterizing this relationship remained a holy grail throughout the project! You can read more about this part of the project in Conor’s MakerZone blog post.
August: Algorithms
9.5 weeks to go
With promising early results, we kicked off the data analysis part of the project. We started by brainstorming the elements we would like to understand about the social interaction patterns. Some of our favorite suggestions: Can we quantify the relationship between caffeine intake and social activity? Can we detect wallflowers, conversation boosters… or conversation killers? Can we see if visitors to a demo station recommend it to others? We quickly realized that most of those questions needed a definition of a conversation. We decided that two tag holders are in a conversation if they stay near to each other for 30 seconds.
But how do we know whether two tags are close to each other? We decided to use Received Signal Strength Indication (RSSI) as a proxy for distance between two tags. On the simplicityrobustness scale, it is closest to simplicity and that’s what we needed given the project timeframe.
By detecting conversations between tags, we can learn the popularity of demo and catering stations (see Figure 2) by monitoring tags located at those places of interest. All the tag data would be transmitted to Raspberry Pi base stations, which would send it in batches to the analytics server every 10 seconds.
Figure 2: EXPO Exhibition Hall floor plan with marked demos and catering stations (click to enlarge)
We had 10 engineers helping on the data analysis challenge in August. However this posed its own challenges:
We needed to design a system that would allow numerous people with different backgrounds and levels of experience to:
We put in place a 5point plan to address these challenges:
Figure 3: Highlevel system architecture: uniform access to real and simulated data; analysis agents representing system logic and producing aggregate data about the event to be displayed in MATLAB or via web interface
In the meantime, the hardware team was busy creating in reality what we had simulated: the networking system that would relay the tag data to a server.
By the end of August we had early versions of all the analysis agents, tested on data from the simulator. We also had a long list of known issues, including the fact that our prototype code could not keep up to the predicted data rate.
September: “Production”
5 weeks to go
September saw the transformation of our prototype system into an early version of what could be called a production system.
First, we harnessed our code in a series of unit tests (better late than never!) to allow faster iterations when fixing bugs and optimizing the code.
Second, we optimized the prototype implementations for speed. A couple of the biggest wins were:
Vectorization of code of calculation of the social network graph Implementation of some of the tag data filtering in SQL instead of MATLAB to reduce time overhead of data transfersBoth operations reduced the time needed to analyze a batch of data from 10s of seconds to fractions of seconds.
Third, we reimplemented our analysis engine to run each analysis agent on a separate thread. Out of many choices for doing so, we decided to recode the thin analysis engine skeleton in Java and deploy each analysis agent as a separate function on MATLAB Production Server. Thanks to using the loosely coupled agents in the first place, this required only changing one line in the agent’s code.
As a final step we used a machine on EC2 with an elastic IP to set up a Tomcat webserver with our webpage and a directory that would be regularly synchronized with our local directory containing all the results.
At the end of September we had a system analyzing the tags every 10 seconds and updating the webpage every minute but still working only with the simulator. Throughout September attempts to connect the two failed – due to an Endian disagreement, too slow a processing rate in the early version of the analysis engine, data overflow in the prototype data gathering stage, and more.
But we were not ready to cancel the project. In the end, we would not often have an opportunity to test it in a hall with 500 people, 8 catering stations, 9 demo stations and a BLOODHOUND SSC lifesize model in the middle (yes, a big pile of metal in the middle of our RF experiment).
October: EXPO or The First Integration Test
1 week to go
First week of October has gone between integration attempts and bug fixing. You can see that the moods were mixed:
Figure 4: Rainbow of emotions; Note as moods are improving when we move from bug chasing to handling some real hardware
The day before EXPO (October 7th) three of us set off to Silverstone, the conference venue. Two hours after arriving everything was connected. The moment the webpage started to show the two test tags in “conversation” with each other was worth all those sleepless nights. We stared in awe. By 7pm all 8 base stations were connected and receiving. At that point we put the system to sleep and left the site.
The morning of EXPO passed in a flash. We started with a pile of 200 tags (see Figure 5). By 9.30 am all of them has been distributed. Attendees were curious about the system and wanted to participate. The webpage started to update with plausible results (see sample in Figure 7). All team members, most of whom so far have seen only a small part of the system, were hypnotized by the updating web page. Throughout the day the demo drew a big crowd who asked about everything from working with low cost hardware, through machine learning and “big data” to web deployment. We also highlighted parts of the system during all relevant talks: Machine Learning, Accelerating MATLAB algorithms, Deployment, Modelling, Simulation, and RealTime Testing for ModelBased Design. Finally, in the wrap up talk we announced the tag number that won the networking competition prize. We knew the system was working when we discovered that the winner was our excolleague Michael, known for being a social butterfly.
So through this project we achieved:
But as EXPO was also the first time the system ran together, we still had a lot of work to understand what it was really doing.
Learn more as the story continues: Expo Conversations – Part 2
Figure 5 Tags waiting to be handed to attendees
Figure 6: Queue to pick up the tags
Figure 7: Screenshot of the live social networking analysis