We’ve added a direct connection between File Exchange and GitHub. In practice, this means that if your code already lives on GitHub, you can leave it there. The File Exchange is now smart enough to distribute the code directly from GitHub to anyone who asks for it from MATLAB Central. You get all the benefits of collaborative development in GitHub while community members get access to the latest version of your projects. No need to maintain the files in two different locations. Push your changes up to GitHub, and seconds later that’s what gets served from the File Exchange.
We’re very excited about this feature, because we know there are already a lot of excellent MATLAB-related repos on GitHub. With this feature, we join two strong communities and open a new era on MATLAB Central.
]]>Sean's pick this week is kml2struct by James Slegers.
Earlier this week, my friend sent me a Google Maps link containing our hiking tracks recorded with the GPS on his smart phone. From Google Maps you can download the data as a KML (Keyhole Markup Language) file.
I wanted to plot it and experiment with the data in MATLAB. Once again, the File Exchange was there for me!
% Read the KML file into a struct: kmlS = kml2struct('2014-04-12PresidentialTraverse.kml'); % Convert to table (new datatype in R2013b) to make manipulations easier kmlT = struct2table(kmlS);
Now looking at the table, we can see the four important pieces:
First, I'll get the bounding box of the whole trip. To do this, we'll stack each segment's bounding box into the third dimension and then pick the min and the max:
boxes = kmlT.BoundingBox; % Extract Bounding box from table boxes3d = cat(3,boxes{:}); % Stack along third dimension bbox = [min(boxes3d(1,:,:),[],3); max(boxes3d(2,:,:),[],3)].'; % Min and max along third dimension give limits latlim = bbox(2,:)+[-0.01 0.01]; % Buffer them lonlim = bbox(1,:)+[-0.01 0.01];
Next, I only want to work with the lines, i.e. the actual tracks. The points represent termini, which I don't need right now. Using the new categorical data type and logical indexing, we can extract the latitude and longitude from the table.
% Make Geometry categorical kmlT.Geometry = categorical(kmlT.Geometry); % Extract the latitude and longitude for the lines latlon = kmlT{kmlT.Geometry=='Line', {'Lat','Lon'}};
I'll get the elevation data from NASA using the Web Map Service in the Mapping Toolbox.
nasaLayers = wmsfind('nasa*elev', 'SearchField', 'serverurl'); ned = refine(nasaLayers, 'usgs_ned'); [Z, refmatZ] = wmsread(ned, 'Latlim', latlim, 'Lonlim', lonlim); Z = double(Z);
And finally, plot a contour map with the tracks overlaid on it.
figure ax = usamap(latlim, lonlim); geoshow(Z, refmatZ, 'DisplayType', 'texturemap') contourm(Z, refmatZ, 20, 'Color', 'k') demcmap(Z) title('Presidential Traverse 04/12:13/2014','FontSize',16) % Add each segment for ii = 1:length(latlon) geoshow(latlon{ii,:}, 'LineWidth', 2) end
Have you ever recorded a trip and then tried to analyze it in MATLAB? The tasks above would be more straight-forward if you had access to the original GPX files. These typically come with the elevation data and time stamps allowing you to get even more statistics with more accuracy.
Give it a try and let us know what you think here or leave a comment for James.
Get
the MATLAB code
Published with MATLAB® R2014a
plot(r287(ind:ev),hhGrams)
or
startIndex = ind;
stopIndex = ev;
xVec = r287;
yVec = hhGrams;
plot(xVec(startIndex, stopIndex), yVec)
All I did was change the names of those variables for the plotting step. Ideally, I would have had well named variables throughout my code, but at least this little corner of the code it more understandable now. This kind of simplification (even though it is more lines of code), when applied many, many times will tame complexity.
Simulink Projects Git Integration
Simulink Projects was introduced in R2011b to support effective collaboration across teams. Until now, the focus has been teams within an organization. With the introduction of Git integration in R2014a you can expand your team using social collaboration sites such as GitHub.
This new functionality is ideal for crowd-sourced development projects. With that in mind, I published our Simulink model of the Hyperloop onto GitHub. Here’s how I did it.
Creating a new Simulink Project on GitHub
To create a new Simulink Project to be shared on GitHub, I had to start by creating the new repository on GitHub. The URL of that repository is needed when creating the Simulink Project. Starting at https://github.com/new I went ahead as below:
Once I had the empty repository on GitHub, I went to MATLAB and set up a New > Simulink Project > From Source Control. This opened the Project Retriever window below. I set Git as the Source control integration and pasted the URL as the Repository path.
Once I had the local Simulink Project, I could begin adding content in my sandbox (or “working tree” if you’re fluent in git-ish). To publish the content, I first had to Commit the modified files to my local repository and then Push that repository out to GitHub.
Contributing to an existing GitHub repository
Although GitHub is designed for open, social collaboration, it’s not completely without access restrictions. To modify a repository on GitHub, you must be the owner or defined as a collaborator. Collaborators are defined in the Settings of the repository.
If you’re interested in making contributions to an interesting repository (a revolutionary transportation concept, for example), you must first create your own Fork, or copy.
If you’re confused about Branches and Forks, take a look at the graphic below. If you’re like me, you’ll need to print it out and sleep with it under your pillow.
Once you’ve created your own Fork, follow the same steps as a new repository to load the project into Simulink. You can freely Push your changes out to your own Forked repository.
If you think your changes are really great and should be included in the original repository, create a Pull Request in GitHub. The owner of the original repository will be alerted via e-mail about your desire to contribute.
Accepting contributions
I was lucky enough to have a contributor (some Guy who knows a lot about SimMechanics) improve my calculation for passenger g-forces in their local frame. You can see the e-mail that I received when he submitted a Pull Request.
After accepting the Pull Request in GitHub, I opened my Simulink Project and Fetched the latest archive. But, when I opened up the model, I didn’t see Guy’s updates. That’s because the files in the project window are my Sandbox, which is the Current Branch. I need to merge the updates using Manage Branches. In the Branches pull-down, I found a branch called refs/remotes/origin/master, which is the remote repository that I just Fetched. After clicking Merge I was working with Guy’s contribution.
Now it's your turn
Retreive or fork from Matt's Hyperloop repository, or create your own Simulink Project on GitHub and let us know by leaving a comment here.
]]>The Singular Value Decomposition of the digram frequency matrix of a text coded with a simple substitution cypher can reveal information about the vowels and consonants in the text.
I have mentioned Don Morrison before in Cleve's Corner, in my post on Pythagorean Addition. He was the mathematician and computer scientist who recruited me to the University of New Mexico in 1972. He was remarkably imaginative and today's post is about another one of his creative ideas. We wrote a paper about it, cited below, in the American Math Monthly thirty years ago.
In English, and in many other languages, vowels are usually followed by consonants and consonants are usually followed by vowels. This fact is reflected by the principal component analysis of what I will call the digram frequency matrix for a sample of text. (I am using the term digram to refer to any pair of letters, although in the discipline of orthography the term more precisely refers to a pair of characters that represents one phoneme.)
English text uses 26 letters, so the digram frequency matrix is a 26-by-26 matrix, $A$, with counts of pairs of letters. Blanks and all other punctuation are removed from the text and the entire sample is thought of as circular or periodic, so the first letter follows the last letter. The matrix entry $a_{i,j}$ is the number of times the $i$-th letter is followed by the $j$-th letter. The row and column sums of $A$ are the same; they count the number of times individual letters occur in the sample. The fifth row and fifth column usually have the largest sums because the fifth letter, which is "E," is usually the most frequent.
A principal component analysis of $A$ produces a first component,
$$ A \approx \sigma_1 u_1 v_1^T $$
that reflects the individual letter frequencies. The first right- and left-singular vectors, $u_1$ and $v_1$, have elements that are all of the same sign and that are roughly proportional to the corresponding frequencies.
We are primarily interested in the second principal component,
$$ A \approx \sigma_1 u_1 v_1^T + \sigma_2 u_2 v_2^T $$
The second term has positive entries in vowel-consonant and consonant-vowel positions and negative entries in vowel-vowel and consonant-consonant positions. The sign patterns adjust the frequency counts to reflect the vowel-consonant properties of the language.
My Numerical Computing with MATLAB toolbox ncm includes a function, digraph.m, that carries out this analysis. Let's run the function on Lincoln's Gettysburg Address, gettysburg.txt.
fid = fopen('gettysburg.txt');
txt = fread(fid);
fclose(fid);
Convert to integers between 1 and 26.
k = upper(char(txt)) - 'A' + 1;
k(k < 1 | k > 26) = [];
Use the counting feature of sparse to generate the digram frequency matrix.
j = k([2:length(k) 1]); A = full(sparse(k,j,1,26,26)); cnt = sum(A);
Compute the SVD.
[U,S,V] = svd(A);
Plot the second left and right singular vectors. The $i$-th letter of the alphabet is plotted at coordinates $(u_{i,2},v_{i,2})$. The distance of each letter from the origin is roughly proportional to its frequency, and the sign patterns cause the vowels to be plotted in one quadrant and the consonants to be plotted in the opposite quadrant.
There is more detail. The letter "N" is usually preceded by a vowel and often followed by another consonant, like "D" or "G," and so it shows up in a quadrant pretty much by itself, sort of a "left vowel, right consonant" . On the other hand, "H" is often preceded by another consonant, namely "T," and followed by a vowel, "E," so it also gets its own quadrant, a "right vowel, left consonant".
for i = 1:26 if cnt(i) > 0 text(U(i,2),V(i,2),char('A'+i-1)) end end s = 4/3*max(max(abs(U(:,2))),max(abs(V(:,2)))); axis(s*[-1 1 -1 1]) axis square line([0 0],[-s s],'color','b') line([-s s],[0 0],'color','b') box title(sprintf('%d characters',length(k)))
A simple substitution cipher is obtained by permuting the letters of the alphabet. The digram matrix $A$ is replaced by $PAP^T$ where $P$ is a permutation matrix. The singular vectors are simply multiplied by $P$. Letters in the plot are permuted, but the locations and quadrants remain the same. So it is possible to identify the vowels and consonants, and perhaps -- if the sample is large enough -- the letters standing for "N" and "H" in the cryptogram.
I admit that simple substitution ciphers are not the world's most difficult codes to crack, but I found this to be an unexpected use of the SVD.
Cleve Moler and Donald Morrison, "Singular Value Analysis of Cryptograms", The American Mathematical Monthly, Vol. 90, No. 2, 1983, pp. 78-87, <http://www.jstor.org/stable/2975804>
Get
the MATLAB code
Published with MATLAB® R2014a
Idin's pick this week is Communications System Toolbox Support Package for RTL-SDR Radio by MathWorks Communications System Toolbox Team.
There has been a proliferation in the availability and usage of low-cost hardware platforms such as Raspberry Pi, BeagleBoard, and USRP boards. These boards can be used for quickly realizing ideas in the physical world, as well as great educational tools. MathWorks has invested heavily over the last few years to work with as many of these platforms as possible, which makes my life as a MathWorks application engineer much more exciting (see list of supported hardware here).
The USB-based software-defined radio (SDR) card from NooElec selling for under $25 is a nifty piece of hardware for any wireless communications enthusiast. There are a few different versions of this radio, but I will focus on the one pictured below (officially: NooElec NESDR Mini SDR and DVB-T USB Stick (R820T)).
What makes this radio more accessible is the MATLAB support package, which allows off-the-air signals to be captured directly into MATLAB and/or Simulink. This means you can make simple receivers to demodulate/decode signals that are present around us.
What can you do with this radio?
The radio front-end on this card supports frequencies from 25MHz to 1.7GHz, which means you can receive FM radio, some navigation signals, and most GSM and LTE bands in use in the United States. The bandwidth is limited to about 3MHz, so you can easily cover FM radio or GSM, but you can't do a full LTE receiver (most LTE deployments go up to 10MHz). But even on LTE, you can still capture and decode the LTE broadcast information (MIB), which is transmitted within the middle 1.4MHz of the band (the middle 72 OFDM subcarriers). If you have access to LTE System Toolbox, this should be a relatively simple exercise.
What's already there?
When you download the MathWorks support package, you will get a few example models including an FM receiver (for both MATLAB and Simulink), a spectral analysis app, and an FRS/GMRS walkie-talkie receiver. So you can listen to your favorite FM radio stations, or receive signals from an off-the-shelf walkie-talkie.
How does it work?
My favorite piece about this support package is its simplicity. With the support package, you simply get a "receiver block" in Simulink, and a corresponding System object in MATLAB. In Simulink, a simple spectrum analyzer looks like this:
where you set three parameters for the radio (as shown below): its center frequency, sampling rate, and tuner gain.
The MATLAB interface is just as straight forward.
% Create receiver System object hSDRrRx = comm.SDRRTLReceiver(... 'CenterFrequency',102.5e6, ... 'EnableTunerAGC',true, ... 'SampleRate',1e6, ... 'SamplesPerFrame',256*20, ... 'OutputDataType','double');
% View spectrum; assuming "hSpectrum" is already initialized for count = 1:5000 [data,~] = step(hSDRrRx); % no 'len' output needed for blocking operation data = data-mean(data); % remove DC component step(hSpectrum,data); end
For the curious, the spectrum analyzer object, hSpectrum, can be initialized as follows:
hSpectrum = dsp.SpectrumAnalyzer(... 'SpectrumType','Power density',... 'SampleRate',1e6, ... 'SpectralAverages',50);
So What?
I would argue that every electrical engineering student taking a digital communications course must own one of these devices. Even for the non-student, and the generally curious person, this will be the best $25 you will spend this year. And the MATLAB support package for the radio makes this one of the cheapest and easiest learning tools for wireless communications.
Comments
Let us know what you think here or leave a comment for MathWorks Communications System Toolbox Team.
Get
the MATLAB code
Published with MATLAB® R2014a
Our colleagues Jay Abraham and Corey Lagunowich saw the recent news stories about Apple's iOS security vulnerabilities, and looked a little closer at how to avoid these problems using Polyspace.
goto fail
Back in February, technology blogs and news outlets (see this article from The Guardian for an example) were abuzz about a newly discovered vulnerability in Apple’s iOS and Mac OS X. There was a problem in the Transport Layer Security (TLS) and Secure Sockets Layer (SSL) code that could be exploited by what is known as a Man in the Middle attack (MitM). The vulnerability was dubbed “goto fail”.
The Anatomy of the Defect
The defect originated in the TLS/SSL library in source file sslKeyExchange.c (accessible on www.opensource.apple.com). The stripped code snippet that represents the defect is shown in the screenshot below.
The problem with this code is that there are two repeated goto statements, and the second is not surrounded by an if statement. This means that any code after the second “goto fail” will never be executed (and is referred to as “unreachable code” or “dead code”). Therefore the function SSLHashSHA1.final() was never called. This scenario provided a set of conditions that could be exploited for an attack.
Presence, Not Absence
How do software bugs like this escape into production? In order to answer this question, let’s turn to the lessons taught by a founding member of the computer scientist community, Edsger W. Dijkstra.
Dijkstra is considered to be a leading contributor to the bedrock principles of computer science in the areas of programming languages and formal verification. Of all his tenets, perhaps the most famous is:
In theory, performing code review and executing the right set of test cases can catch most code defects. But this is hard in practice. Even for the simplest operations, such as adding two 32-bit integers, one would have to spend significant time to complete exhaustive testing. And if you don’t test for it – you won’t find it.
So how do we apply Dijkstra’s principle in order to develop safe and secure embedded systems? How do we show the absence of bugs? We use Static Code Analysis, like what is performed by Polyspace Bug Finder and Polyspace Code Prover.
Detecting Defects
How would Polyspace have stopped defects like these from being released? The screenshot below shows the results of Polyspace Bug Finder’s analysis on the code snippet. It flags the final if statement as dead code, because it cannot be reached.
The “goto fail” bug is a fairly straightforward example of this. What about something more complicated, something that requires being able to follow the data flow to detect? For that we can use Polyspace Code Prover.
A More Complex Example
The example below shows code where Polyspace Code Prover has flagged three items as not green (i.e. not proven reliable):
Because Polyspace can follow the data flow, it knows that the if-statement on line 16 will always evaluate to false. If count ever did equal zero, a run-time error would have occurred on line 13 before the program reaches line 16. And Polyspace Code Prover can detect all this without requiring program execution, code instrumentation, or test cases.
Takeaways
The Apple “goto fail” is a teaching moment for embedded software engineers that develop and maintain critical automotive ECU, avionics, and medical devices. It highlights the need to augment code review and testing with static code analysis. You can look to Polyspace to scrutinize critical software components, find bugs, and ensure that no new bugs creep into your code.
Give look at what Polyspace can do for you and let us know what you think by leaving a comment here.
I've seen two questions recently about the speed of the fft function in MATLAB. First, a tech support question was forwarded to development. The user wanted to know how to predict the computation time for an FFT of a given length, N. This user was interested in values of N in the neighborhood of 4096 (2^12).
The second was a post in the MATLAB newsgroup comp.soft-sys.matlab. This user wondered why padding to the next power of two wasn't always the fastest way to compute the FFT.
Inspired by these questions, I want to show you today how to do some FFT benchmarking in MATLAB.
It turns out that, in general, the time required to execute an N-point FFT is proportional to N*log(N). For any particular value of N, though, the execution time can be hard to predict and depends on the number of prime factors of N (very roughly speaking). The variation in time between two close values of N can be as much as an order of magnitude.
Whenever I do FFT benchmarking, I find it very helpful to look at three sets of numbers:
Also, I have learned to look at plots that are log scale in N and that have roughly the same number of test values within each octave (or doubling) of N.
Constructing sets of N values along these lines takes a little thought. Here's some code.
First, how many powers of 2 do we want to examine? Based on the customer questions I saw, I want to examine the range from 1024 (2^10) to 131072 (2^17).
low2 = 10; high2 = 17; n2 = 2.^(low2:high2);
Next, I want to pick 10 composite numbers and 10 prime numbers in between successive powers of 2. I'd like to pick the numbers "randomly," but I also want my experiment to be repeatable. To satisfy these seemingly contradictory constraints, I'll reset the MATLAB random number generator before beginning.
rng('default'); % Initialize the vectors holding the prime N's and composite N's. np = []; nc = []; for m = low2:high2 k = (2^m):(2^(m+1)); kp = k(2:end-1); isp = isprime(kp); primes = kp(isp); composites = kp(~isp); % Use randperm to pick out 10 values from the vector of primes and 10 % values from the vector of composites. new_np = primes(randperm(length(primes),10)); new_nc = composites(randperm(length(composites),10)); np = [np new_np]; nc = [nc new_nc]; end
Now let's use the function timeit to measure the execution time required to compute FFTs for all these values of N. (If you don't have a recent version of MATLAB that has timeit, you can get a version of it from the File Exchange.)
t2 = zeros(size(n2)); for k = 1:length(n2) x = rand(n2(k),1); t2(k) = timeit(@() fft(x)); end
tp = zeros(size(np)); for k = 1:length(np) x = rand(np(k),1); tp(k) = timeit(@() fft(x)); end
tc = zeros(size(np)); for k = 1:length(nc) x = rand(nc(k),1); tc(k) = timeit(@() fft(x)); end
Now do a loglog plot of all these times.
loglog(n2,t2,'o') set(gca,'xtick',2.^(10:17)) xlim([2^10 2^17]) hold on loglog(nc,tc,'+') loglog(np,tp,'*') hold off legend({'Powers of 2','Composite numbers','Prime numbers'}, ... 'Location','NorthWest') xlabel('N') ylabel('Execution time (s)') title('FFT execution time as a function of N')
You can see that there's a wide spread of execution times for the values of N that are not powers of 2.
One thing I'm not seeing is what the MATLAB Newsgroup poster reported. That is, I don't see a non-power-of-2 time that's faster than the next highest power of 2.
So let's look a little harder for composite numbers that are faster than what we've seen so far. Specifically, I'm going to look for values of N with prime factors no bigger than 3.
nn = []; for m = low2:high2 k = (2^m):(2^(m+1)); kp = k(2:end-1); kp = kp(randperm(length(kp))); nn_m = []; for q = 1:length(kp) if max(factor(kp(q))) <= 3 nn_m = [nn_m kp(q)]; if length(nn_m) >= 4 % We've found enough in this part of the range. break end end end nn = [nn nn_m]; end
Measure execution times for these "highly composite" numbers.
tn = zeros(length(nn),1); for k = 1:length(nn) x = rand(nn(k),1); tn(k) = timeit(@() fft(x)); end
Add the times to the plot.
hold on loglog(nn,tn,'d') hold off legend({'Powers of 2','Composite numbers','Prime numbers', ... 'Highly composite numbers'},'Location','NorthWest')
You can see that sometimes a non-power-of-2 can be computed very fast, faster than the next higher power of 2. In this experiment we found one such value of N: 39366. This number has 10 prime factors:
factor(39366)
ans = 2 3 3 3 3 3 3 3 3 3
I hope you enjoyed these experiments with FFT benchmarking. I can tell you from personal experience that it can turn into almost a full-time hobby!
Get
the MATLAB code
Published with MATLAB® R2014a
With spring comes the tax filing deadline. This post is also about taxes. I'd like to introduce this week's guest blogger, Toshi Takeuchi. Toshi analyzes web data and runs online ad campaigns here at MathWorks.
Hi, I am Toshi. I am a big fan of Nate Silver who made analyzing data very cool and mainstream. Because I analyze data a lot, it bugs me when I see questionable analyses passed around in the news media.
So when I saw this CNBC post on Google+, my “bogus data analysis” radar started sending high alerts.
This map shows you the ranking of states based on average tax amount, adjusted to the cost of living index. Let me pretend some data journalism here.
Well, I happen to think that the tax amount is correlated more directly to income, rather than cost of living. The average tax amount should be higher if you live in a state with high median income. Cost of living may be also higher in those states, but that's a secondary effect.
In order to understand the true picture, you actually need to think in terms of tax to income ratio instead. This is what you get when you use this metric.
You can see that the color shifted in a number of states if you compare the first map and the second map. Massachusetts, where I live, actually looks pretty good; so are states in Mid-Atlantic region, while they were red in the first map. On the other hand, the color flips in the other direction in some Gulf Coast states in the South.
If you believed in the original analysis and moved from Massachusetts to one of those states, then your taxes may go down, but your income may also go down even more. Not a good move, IMHO.
% Disclaimer: don't trust my analysis, either - I only did it for % debunking the original story; it is not meant as a robust analysis. % Please don't plan your relocation based on this analysis, just in case.
If you are interested in playing a data journalist, this type of analysis is fairly easy with MATLAB.
All I had to do was to download the median income dataset from the Census Bureau website and merge two datasets with the newly introduced table (available in MATLAB since R2013b). I also used Mapping Toolbox to visualize the data.
First, I went to the data sources to get the data. You can use Excel to import HTML tables into a spreadsheet directly. Census data is also available in Excel format. To match the time period with the original analysis, I used
Historical (1984 to 2012): Median Household Income by State - Single-Year Estimates [XLS - 98k]
Data sources
% load data from files into tables Tax = readtable('bestWorstStateTax.csv'); Income = readtable('medianIncome.csv'); % inspect the content disp Tax disp(Tax(1:5,:)) disp Income disp(Income(1:5,:))
Tax Rank State AvgAnnualStateLocalTaxes ____ ______________ ________________________ 1 'Wyoming' 2365 2 'Alaska' 2791 3 'Nevada' 3370 4 'Florida' 3648 5 'South Dakota' 3766 PercentDiffFromNationalAvg AdjRankCostOfLivingIdx __________________________ ______________________ -0.66 1 -0.66 4 -0.52 2 -0.48 3 -0.46 5 Income State MedianIncome StandardError ____________ ____________ _____________ 'Alabama' 43464 2529.4 'Alaska' 63648 2839.1 'Arizona' 47044 2921.7 'Arkansas' 39018 2811.5 'California' 57020 1237.5
Now we have two tables in the workspace, and you can also see that each column has a header and can contain a different data type. Both tables contains the same column called "State" containing the text string of state names. We can use that as the key to join those two tables. We don't need all the columns for this analysis, so I will join just the columns I need.
% |table| is smart - it automatically uses that "State" column as the key. % Just using |State| and |AvgAnnualStteLocalTaxes| and |State| and % |MedianIncome|. T1 = join(Tax(:,2:3),Income(:,1:2)); % rename columns T1.Properties.VariableNames = {'State','Tax','Income'}; % compute tax to income ratio T1.Ratio = T1.Tax./T1.Income; % create a new table ranked by tax to income ratio T2 = sortrows(T1,{'Ratio'}); % inspect the new table disp T2 disp(T2(1:5,:))
T2 State Tax Income Ratio ______________ ____ ______ ________ 'Wyoming' 2365 57512 0.041122 'Alaska' 2791 63648 0.043851 'Washington' 3823 62187 0.061476 'Nevada' 3370 47333 0.071197 'South Dakota' 3766 49415 0.076212
Check whether the new metric produced any meaningful differences.
disp('Top 20') disp('By Avg. Tax By Avg. Ratio') disp([T1.State(1:20) T2.State(1:20)])
Top 20 By Avg. Tax By Avg. Ratio 'Wyoming' 'Wyoming' 'Alaska' 'Alaska' 'Nevada' 'Washington' 'Florida' 'Nevada' 'South Dakota' 'South Dakota' 'Washington' 'Florida' 'Texas' 'Colorado' 'Delaware' 'Texas' 'North Dakota' 'North Dakota' 'Colorado' 'Utah' 'New Mexico' 'Delaware' 'Alabama' 'Massachusetts' 'Arizona' 'New Hampshire' 'Utah' 'Virginia' 'Mississippi' 'Maryland' 'Indiana' 'District of Columbia' 'Louisiana' 'Rhode Island' 'West Virginia' 'Arizona' 'Montana' 'Hawaii' 'Oklahoma' 'New Jersey'
Now we will start using the functions from Mapping Toolbox. First we will assemble the required pieces of data to prepare a map.
Note: I also used Bioinformatic Toolbox function redgreencmap to create the colormap to go from green to red to mirror the scheme in the original map. If you don't have this toolbox, you can easily create a custom map in MATLAB. Colormaps are arrays of RGB values (triplets) in the range from 0 to 1.
% get the US geography data as a structure array states = shaperead('usastatelo', 'UseGeoCoords', true); % Get the state names as a cell array of strings. names = {states.Name}; % This is a vector that the stores ranking of each state. ranking = zeros(length(names),1); for i=1:length(names) ranking(i)=find(strcmpi(names(i),T2.State)); end % Create a colormap that goes from green to red in 51 steps. colors = redgreencmap(length(ranking)); % Sort colors by state ranking. stateColors = colors(ranking,:); % Separate Hawaii and Alaska from the Continental US. indexHawaii = strcmp('Hawaii',names); indexAlaska = strcmp('Alaska',names); indexConus = 1:numel(states); indexConus(indexHawaii|indexAlaska) = [];
Now we are ready to draw the map.
% This creates a figure with axes of US geography. % It contains three axes - Continental US, Alaska and Hawaii. figure; ax = usamap('all'); % We don't need the axes, so turn them off. set(ax, 'Visible', 'off') % Draw the states with specified color within the Continental US. for j = 1:length(indexConus) geoshow(ax(1), states(indexConus(j)),'FaceColor',stateColors(indexConus(j),:)) end % Now do the same for Alaska and Hawaii. geoshow(ax(2), states(indexAlaska),'FaceColor',stateColors(indexAlaska,:)) geoshow(ax(3), states(indexHawaii),'FaceColor',stateColors(indexHawaii,:)) % We don't need geographical details, so turn them off for each axes. for k = 1:3 setm(ax(k), 'Frame', 'off', 'Grid', 'off',... 'ParallelLabel', 'off', 'MeridianLabel', 'off') end % Add a colorbar. colormap(flipud(colors)) c= colorbar('YTickLabel',... {'51','41',... '31','21','11','1'}); ylabel(c,'Ranking')
You can download the data files from these links:
Did you enjoy seeing how to use MATLAB for debunking some bad news analysis? Would you like to try? Perhaps you already do this yourself. Tell us know about it here.
Get
the MATLAB code
Published with MATLAB® R2014a
Since this is the first major overhaul of the Blogs site, we’re all pleased to see the new changes released. In fact, we’d like to give you a tour of the Blogs site to show off some of the new features that we’ve added. So sit back and enjoy the view as we take a trip through the new Blogs site.
If you’re a regular reader of the blogs, chances are you haven’t been to the landing page in a while. The last time you went, it probably looked like this:
Now, the landing page looks like this:
In addition to sporting big, colorful images of our bloggers, there are a few other things we’d like to point out:
Go to the home page of Loren Shure’s blog, and you’ll see several new features:
Okay, let’s get to the good stuff and see Loren’s post on Intersecting Lines.
In addition to this terrific post by Loren, you’ll notice a couple of other things:
What do you think of the new Blogs site? What do you like about it? What’s still missing? Leave us a comment below.
]]>Jiro's pick this week is 2048 MATLAB.
"2048"
Need I say more? I'm certain that many of you know what I'm talking about. 2048 is an online and mobile game created by Gabriele Cirulli recently. All I can say is that it is an addictive game that has already resulted in a few MATLAB implementations.
A couple of weekends ago, I also spent a day implementing a MATLAB version. I probably spent more time than I needed to on trying to make it look as close as possible to the original version, but it was a lot of fun.
But then, one of my colleagues, Seth (who was also just featured in this blog last week), suggested including a mechanism for feeding in different algorithms to solve the game. There must be an AI algorithm that would solve the game or produce a high score consistently.
In the app, I've included the ability to select a MATLAB function as the solving algorithm. Then you can see the algorithm in action as it tries to solve the game. The algorithm must be defined as a MATLAB function that takes in a 4x4 matrix representing the game board (NaN for empty spots) and returns a character array representing the direction to move ('up', 'down', 'right', or 'left'). Here's an example of a simple algorithm that randomly chooses a direction.
function direction = myAI(board) d = {'up', 'down', 'right', 'left'}; direction = d{randi(4)};
I've also included a command line simulator, which you can use to do a montecarlo simulation of a particular algorithm.
game = Game2048Simulator(@myAI); simulate(game, 1000) % 1000 simulations viewResult(game, 25) % 25 histogram bins
We can display the highest score and compute what percentages of the simulation resulted in each of the highest block.
% Display highest score disp(['Highest score: ', num2str(max(game.Result.Score))]) disp(' ') % Get the unique values for high blocks highblocks = game.Result.HighestBlock; blocks = unique(highblocks); % Compute the percentage of occurrance for each high block value nSims = length(highblocks); for id = 1:length(blocks) percentage = nnz(highblocks == blocks(id)) / nSims; disp([num2str(blocks(id)), ' ', num2str(percentage*100), '%']) end
Highest score: 3172 16 0.2% 32 6.8% 64 34.5% 128 49.7% 256 8.8%
Submit your algorithm!
So here's a challenge for you. Submit your algorithm and win some MathWorks swag! I will give prizes for each of the following criteria:
It was brought to my attention that since the simulator I've included with the File Exchange entry would stop the game once you reach 2048, I am limiting the highest score you could get. For my next update, I will change it so that the game will continue until no more blocks can move.
If your algorithm is fairly short, please submit the algorithm through the comments for this blog post. Be sure to use code formatting. For lengthy algorithms or multi-function algorithms, feel free to email them to me (as ASCII text).
In the next blog post I write, I will go over some of the winning algorithms you have submitted.
Enjoy!!
Get
the MATLAB code
Published with MATLAB® R2014a
The Cure for Tuning Headaches
Simulink makes it easy to model and simulate feedback control systems. But how do you pick the gains of your controller to get adequate performance and robustness? Simple enough when tuning a single PI loop, harder for control systems with multiple loops, configurations, and operating conditions. You probably use a combination of know-how, experience, trial-and-error, and home-grown tools. Wouldn’t it be nice if you could just enter your specifications and let the computer figure out the gain values? Welcome to systune!
Tuning Workflow in Simulink
To see how this works, let’s tune a cascade controller for setting and regulating the speed of a DC motor. This controller consists of two feedback loops: an inner loop for controlling the current in the armature, and an outer loop for controlling the motor speed. Both loops use digital proportional-integral (PI) controllers, so there is a total of four gains to tune.
The first step is to create a slTuner object for interacting with your Simulink model and to specify which blocks you want to tune:
ST0 = slTuner('DCMotor',{'SpeedController','CurrentController'});
Next, you list the signals and points of interest for tuning and validation. For example, the reference signal Ref and the output Speed will come handy for specifying how fast the control system should respond and for checking the actual response of the tuned system.
addPoint(ST0,{'Ref' , 'Speed' , 'SpeedMeas' , 'CurrentMeas'});
If these two steps look familiar, it’s because slTuner is just an extension of the slLinearizer interface we discussed in this earlier blog post.
Third and last step, you specify the tuning goals, that is, how the control system is supposed to perform. There is a lot to choose from, with goals ranging from tracking and disturbance rejection to loop shape, stability margins, and minimum closed-loop damping. For the DC motor application, I use the following goals:
I use the StepResp and LoopShape objects from the TuningGoal library to express these goals:
Goal1 = TuningGoal.StepResp('Ref','Speed',0.05); BandWidth = 2*pi*200; % 200 Hz in rad/s Goal2 = TuningGoal.LoopShape('CurrentMeas',BandWidth); Goal2.Openings = 'SpeedMeas';
Note that the inner-loop bandwidth should be evaluated with the outer loop open, so I specify a loop opening at the location SpeedMeas in the Simulink model.
I can now launch the tuning algorithm:ST1 = systune(ST0,[Goal1,Goal2]);
Final: Soft = 0.802, Hard = -Inf, Iterations = 60
Each tuning goal receives a normalized score and systune works on improving the overall score. A final score ≤ 1 means “pass” and a final score > 1 means “fail”. Here systune did well with a score of 0.8, and plotting the tuning goals confirms that the tuned responses are dead on spec:
viewSpec([Goal1,Goal2],ST1)
The slTuner object ST1 represents the tuned control system and I can use it to access other system response for further validation. This is entirely similar to the slLinearizer workflow. For example, I can use getLoopTransfer to compute the stability margins of the current loop:
L = getLoopTransfer(ST1,'CurrentMeas',-1); margin(L), grid on
When I am happy with the linear analysis results, I use writeBlockValue to push the tuned values of the PI gains to the Simulink model and perform additional validation in Simulink.
writeBlockValue(ST1)
Control System Tuner App
With the new Control System Tuner app, you can do all the above without a single line of code! You launch this app from the “Analysis” menu in the Simulink model. Here is what it looks like:
Now It's Your Turn
There is a lot more you can do with systune and slTuner:
Check out the many examples and applications and try it on your control system. And as usual, we'd love to hear your thoughts!
]]>datetick
explicitly.
]]>Employing a factorization based on the least significant singular values provides a matrix approximation with many surprisingly useful properties. This Reverse Singular Value Decomposition, RSVD, is also referred to as Subordinate Component Analysis, SCA, to distinguish it from Principal Component Analysis.
The Singular Value Decomposition of a matrix $A$ is computed by
[U,S,V] = svd(A);
This generates two orthogonal matrices $U$ and $V$ and a diagonal matrix $S$ with diagonal elements $\sigma_k$ that provide the expansion
$$ A = \sigma_1 E_1 + \sigma_2 E_2 + ... + \sigma_n E_n $$
where $E_k$ is the rank one outer product
$$ E_k = U(:,k) V(:,k)' $$
Traditionally, the singular values are arranged in descending order. In contrast, the Reverse Singular Value Decomposition Approximation of rank $r$ is obtained by arranging the singular values in ascending order,
$$ 0 \le \sigma_1 \le \sigma_2 \le ... $$
and then using the first $r$ terms from the expansion.
Here is a function that computes the RSVD for square or rectangular matrices with at least as many rows as columns.
type rsvd
function X = rsvd(A,r) % RSVD Approximation by the Reverse Singular Value Decomposition % rsvd(A,r) approximates A by a matrix of rank r obtained from % the r least significant singular values in ascending order. [m,n] = size(A); [U,S,V] = svd(A,'econ'); k = n:-1:n-r+1; X = U(:,k)*S(k,k)*V(:,k)';
In certain situations, the RSVD can reduce or even eliminate roundoff error. For example, according to its help entry the elmat function hilb attempts to compute
hilb(N) is the N by N matrix with elements 1/(i+j-1)
But the function can only succeed to within roundoff error. Its results are binary floating point numbers approximating the reciprocals of integers described in the help entry. Here is the printed output with n = 5 and the default format short.
format short
H = hilb(5)
H = 1.0000 0.5000 0.3333 0.2500 0.2000 0.5000 0.3333 0.2500 0.2000 0.1667 0.3333 0.2500 0.2000 0.1667 0.1429 0.2500 0.2000 0.1667 0.1429 0.1250 0.2000 0.1667 0.1429 0.1250 0.1111
We are seeing the effects of the output conversion as well as the underlying binary approximation. Perhaps surprisingly, the reverse singular value decomposition can eliminate both sources of error and produce the exact rational entries of the theoretical Hilbert matrix.
format rat
H = rsvd(H,5)
H = 1 1/2 1/3 1/4 1/5 1/2 1/3 1/4 1/5 1/6 1/3 1/4 1/5 1/6 1/7 1/4 1/5 1/6 1/7 1/8 1/5 1/6 1/7 1/8 1/9
The RSVD is capable of uncovering spelling and typographical errors in text. My web site for Numerical Computing with MATLAB has a file with the text of Lincoln's Gettysburg Address.
I've used this file for years whenever I want to experiment with text processing. You can download the file and then use the MATLAB data import wizard with column delimeters set to spaces, commas and periods to produce a cell array, gettysburg, of individual words. There are 278 words. The longest word has 11 characters. So
A = cell2mat(gettysburg)
converts the cell array of strings to a 278-by-11 matrix of doubles.
It turns out that an RSVD approximation of rank nine finds three spelling errors in the original text.
B = rsvd(A,9); k = find(sum(A,2)~=sum(B,2)) disp([char(A(k,:)) char(B(k,:))])
weather whether consicrated consecrated govenment government
In all the years that I have been using this data, I have never noticed these errors.
We have also found that the RSVD is capable of softening the appearance of aggression in photographs. A matrix is obtained from the JPEG image by stacking the RGB components vertically. Roughly 90% of the small singular values then produce a pleasant result.
RGB = imread('creature.jpg'); A = [RGB(:,:,1); RGB(:,:,2); RGB(:,:,3)]; [m,n] = size(A); B = rsvd(A,ceil(0.90*n)); m = m/3; C = cat(3,B(1:m,:),B(m+1:2*m,:),B(2*m+1:3*m,:)) image(C)
Get
the MATLAB code
Published with MATLAB® R2014a
Today we are very happy to give you a sneak peek of a feature we have been working on for a long time:
Now it's your turn
If you wonder when this feature will be released, just look at your calendar.
Let us know what you think by leaving a comment here.
]]>Now for something a little different. As some of you may know, most of the free time that I spend in front of a computer is absorbed by MATLAB Answers. So this week I am going to highlight one of my favorite questions and some of the answers to it.
Sean's pick this week is Seth's Answer to Sven's Question about numerical packing.
I started working at MathWorks in late November, 2011. When Sven asked this question a few weeks later, my wife had not yet moved to Massachusetts. It was December in my temporary not-so-energy efficient house so I was saving costs by keeping the house at 50F (10C) and staying late at work.
This question came in at the end of the work day and immediately caught my attention. I spent the next five hours coming up with an answer.
A year or so later, Seth, our Optimization Toolbox Product Marketing Manager, was showcasing some new technology to the engineers where he debuted the Sudoku Solver that Alan wrote about a few weeks ago and the solution to the Traveling Sales Man problem. This actually interested me more since it brought back memories of Transportation Engineering classes in college. This new technology was still a few releases away from being ready but the potential was already obvious. I sent Seth this question to see if it would be a good candidate and he replied with the optimal answer in about an hour.
When R2014a was release earlier this Month, he agreed to post his answer.
So what is the question? The goal of the question is to optimally group some elements into bins to minimize both the number of bins and the size of the biggest bin. Each bin has a maximum size and the elements can not be rearranged from their original order. This is a type of optimization problem known as a knapsack problem (thanks Walter!).
The first answer to roll in came from Elige Grant. He posted a way to do this with brute force. For small problems, this will be quick enough and the optimal solution will obviously be found and guaranteed.
As the problem size grows, however, brute force starts to become cumbersome in terms of computational resources. My solution was to use the only solver at the time that we had for mixed integer optimization, ga, the genetic algorithm solver. This solver has an option, 'IntCon', to constrain variables to integers, a requirement. I'd also only been able to have access to the Global Optimization Toolbox for the last few weeks, so this was a chance to play with something new. It did work too (!), for small problems. Once again, as the problem size grew, this approach fell apart and there's no way to prove optimality (or lack thereof) from it.
Enter intlinprog.
In R2014a we introduced intlinprog, a mixed-integer linear programming solver. This solver will solve any mixed-integer linear programming problem, be it a minimization or a feasibility study. If requested it will also prove that the solution is optimal for a minimization problem.
Let's examine how it's done:
% The call to |intlinprog| looks like this:
[xopt,fval,eflag,output] = intlinprog(f, intIdxs, A, b, Aeq, beq, lb, ub);
LP: Optimal objective value is 765.000000. Cut Generation: Applied 10 implication cuts. Lower bound is 765.000000. Heuristics: Found 1 solution using rss. Upper bound is 939.000000. Relative gap is 18.50%. Cut Generation: Applied 2 gomory cuts, and 7 implication cuts. Lower bound is 785.000000. Relative gap is 16.37%. Branch and Bound: nodes total num int integer relative explored time (s) solution fval gap (%) 9 0.01 2 9.220000e+02 1.044609e+01 12 0.01 2 9.220000e+02 0.000000e+00 Optimal solution found. Intlinprog stopped because the objective value is within a gap tolerance of the optimal value; options.TolGapAbs = 0 (the default value). The intcon variables are integer within tolerance, options.TolInteger = 1e-05 (the default value).
We can see the sparsity patterns of A and Aeq the linear and equality constraints using spy()
subplot(2,1,1) spy(A) title('Linear Constraints') subplot(2,1,2) spy(Aeq) title('Equality Constraints')
Now let's look at output
output
output = relativegap: 0 absolutegap: 0 numfeaspoints: 2 numnodes: 12 constrviolation: 1.7764e-15 message: 'Optimal solution found.
Since absolutegap is 0, this is the best solution.
If we were curious on time, that took about 0.037s to run!
Do you have a problem that could be solved with mixed-integer linear programming, or one that you're not sure about and would like our thoughts on? Let us know about it here.
Also, if you liked this post, and you're still reading at this point (which probably biases the data): Would you like me to write another out of model Pick Of The Week post that highlights an intlinprog practical example?
Get
the MATLAB code
Published with MATLAB® R2014a
Before R2014a
As your models are getting larger and larger, you might be interested to convert some of the subsystems to Referenced Models.
If you are not familiar with Model Referencing and its advantages, I recommend getting started with the documentation page Component-Based Modeling with Model Reference.
To help converting subsystems to referenced models, we have been providing the function Simulink.SubSystem.convertToModelReference and a menu option which you can see when right-clicking on a subsystem:
For all those who tried this approach, you very likely ran into a series of errors that had to be fixed one-by-one before successfully converting your subsystem into a referenced model.
Improved Workflow
In R2014a, the Model Reference Conversion Advisor is introduced to improve this workflow. When you right-click on the subsystem, this new tool based on the Model Advisor shows up, listing all the steps you will need to go through to convert the subsystem:
Depending on your preferences, you can run the tasks one by one to understand each part of the process, or you can configure the advisor to run all the checks and automatically do all the modifications required by model referencing.
If needed, the Model Reference Conversion Advisor will create and save the bus objects necessary for the referenced model interface. At the end, it will generate a report so you can review the changes that have been automatically applied.
Now it's your turn
Try the improved workflow for converting a subsystem to a referenced model, and let us know what you think by leaving a comment here
]]>MATLAB user Meshooo asked a question on MATLAB Answers about a problem with the createMask function associated with impoly. Meshooo observed a discrepancy between the output of bwboundaries and the mask created by createMask.
I want to describe the issue in more general terms here as a conflict between the geometry of the bwboundaries function and the geometry of the poly2mask function (which is used by createMask).
Here's a simple example that illustrates the discrepancy. Start by creating a small binary image.
BW = [ ...
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 0 0
0 0 1 1 1 1 1 0 0
0 0 1 1 1 1 1 0 0
0 0 1 1 1 1 1 0 0
0 0 1 1 1 1 1 0 0
0 0 1 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 ];
Call bwboundaries, which traces the perimeter pixels of all the objects (and holes) in the image.
B = bwboundaries(BW);
There's only one boundary in this case. Extract and plot it.
B = B{1}; Bx = B(:,2); By = B(:,1); plot(Bx,By) axis ij axis equal axis([.5 9.5 .5 10.5])
If we now pass Bx and By to poly2mask, we don't get exactly the same binary mask image that we started with.
BW2 = poly2mask(Bx,By,10,9)
BW2 = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 1 1 1 1 0 0 0 0 0 1 1 1 1 0 0 0 0 0 1 1 1 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
To understand the reason for the discrepancy, it helps to understand that bwboundaries treats the foreground pixels of the input image as points in space. Here's a plot to illustrate:
plot(Bx,By) axis ij axis equal axis([.5 9.5 .5 10.5]) hold on [yy,xx] = find(BW); plot(xx,yy,'*') hold off legend('Boundary polygon','Foreground pixels')
Now let's do a plot that shows the pixels as squares of unit area. I'll include some code to overlay the pixel edges as gray lines.
imshow(BW,'InitialMagnification','fit') hold on x = [.5 9.5]; for k = .5:10.5 y = [k k]; plot(x,y,'Color',[.7 .7 .7]); end y = [.5 10.5]; for k = .5:9.5 x = [k k]; plot(x,y,'Color',[.7 .7 .7]); end plot(Bx,By,'r') plot(xx,yy,'*') hold off
Now you can see that the polygon produced by bwboundaries does not completely contain the pixels along the border of the object. In fact, most of those pixels are only half inside the polygon (or less).
That's the clue needed to explain the discrepancy with poly2mask. That function treats images pixels not as points, but as squares having unit area. Its algorithm is carefully designed to treat partially covered pixels in a geometrically consistent way. I wrote several blog posts (POLY2MASK and ROIPOLY Part 1, Part 2, and Part 3) about this back in 2006. Here's a diagram from Part 3 that illustrates a bit of the algorithm for handling partially covered pixels.
It turns out that there is a way to get boundary polygons from bwboundaries that are consistent with the poly2mask geometry. The idea is to upsample the binary image so that the polygon produced by bwboundaries is outside the pixel centers instead of running directly through the centers.
BW3 = imresize(BW,3,'nearest');
B3 = bwboundaries(BW3);
B3 = B3{1};
Now shift and scale the polygon coordinates back into the coordinate system of the original image.
Bx = (B3(:,2) + 1)/3; By = (B3(:,1) + 1)/3; imshow(BW,'InitialMagnification','fit') hold on x = [.5 9.5]; for k = .5:10.5 y = [k k]; plot(x,y,'Color',[.7 .7 .7]); end y = [.5 10.5]; for k = .5:9.5 x = [k k]; plot(x,y,'Color',[.7 .7 .7]); end plot(Bx,By,'r') plot(xx,yy,'*') hold off
You can see that the pixel centers are now clearly inside the modified polygon (Bx,By). That means that when we try poly2mask again, we'll get the same mask as the original image.
BWout = poly2mask(Bx,By,10,9)
BWout = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 1 1 1 1 1 0 0 0 0 1 1 1 1 1 0 0 0 0 1 1 1 1 1 0 0 0 0 1 1 1 1 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
isequal(BW,BWout)
ans = 1
To everyone in the Northern Hemisphere: Happy Spring!
Get
the MATLAB code
Published with MATLAB® R2014a
Kristen's pick this week is Packing Santa's Sleigh by Alfonso Nieto-Castanon.
Hi, my name is Kristen, and this is my first post for the Pick of the Week as a guest blogger. My background is in statistics, machine learning, and data science, which are the areas I support as an Education Technical Evangelist at MathWorks. I use MATLAB every day to support faculty and students at universities around the world as they use MATLAB to teach and conduct research. The most exciting part of my job is seeing the diverse problems that are solved with MATLAB! Today, I’m going to talk about one example...
Over the holidays, MathWorks sponsored a Packing Santa's Sleigh competition on Kaggle, a platform for predictive modeling and data analytics competitions. Congratulations to Alfonso, who had the top MATLAB submission and scored second place overall! (Alfonso is no stranger to winning MATLAB contests – see his MATLAB Programming Contest wins here and number one rank on Cody here.)
Participants received the $x$-, $y$-, and $z$-dimensions of an ordered list of one million presents. To solve the problem, they were required to pack the presents onto a 1000 x 1000 sleigh, minimizing the total height of the sleigh and keeping the presents as close to the original order as possible. Presents could be rotated as long as they remained parallel to the $x$-, $y$-, and $z$-axes. (You might recognize this as a holiday version of the three-dimensional bin packing problem.)
The solutions were judged on the compactness of packing (minimized height) and present order (delta from original order). The actual evaluation metric $M$ was calculated by
$M = 2 \cdot max_{i}(z_{i}) + \sigma(\Gamma)$, where
$z_{i}$ = $z$-coordinate of the $i$th present,
$\sigma(\Gamma) = \sum_{i=1}^{N_{p}} |i-\Gamma_{i}|$,
$\Gamma$ = order the presents appear, and
$\Gamma_{i}$ = the Present ID in the $i$th position.
The ideal present order is $\Gamma_{ideal} = 1,2,3,\dots,N_{p}$, so that $\sigma(\Gamma_{ideal}) = 0$, where $N_p$ is the total number of presents.
A naïve solution might be to simply pack the presents in reverse order from the bottom to the top. Each layer would have a height $H$, the next layer would start at that height, and the process would be repeated until all of the presents are in the sleigh. This would optimize the present order, but leaves much room to reduce the height of the sleigh.
Alfonso's solution implements a similar layer-based approach, but he adds a creative twist to how individual layers are defined: He divided the presents in each layer into two sets based on their heights, $H_{1}$ (Set 1) and $H_{2}$ (Set 2). His algorithm keeps the original order of the presents by packing from the bottom up and reverses the solution at the end. Figure 1 shows a naïve layering approach (on the left) and a generalization of Alfonso's solution (on the right).
After the first layer, each layer starts partially occupied by presents from the layer below it. The process for each layer looks like this:
After all of the layers are packed, the algorithm compacts them further by attempting to lower each present vertically, if there is room. Finally, the sleigh is flipped upside down to improve the optimal volume density and achieve the optimal order. This got Alfonso his winning score.
$^{1,2,3}$ Alfonso provides the definitions of these terms, as well as a more detailed description with additional tricks and nice visualization of his solution, on the competition Forum here. You can see other Kaggle competitions that he's won here.
Congratulations to Alfonso, the other competition winners, and all 362 competing teams!
Visit the competition Data page to download the data and sample MATLAB code to get you started. Good luck!
Give it a try and let us know about your solution here or leave a comment for Alfonso!
Get
the MATLAB code
Published with MATLAB® R2014a
Flappy Bird is currently the most popular file on the File Exchange, unseating Oliver Woodford’s venerable champion export_fig for the first time in years. How much longer until export_fig is back on top? Or does the ‘Bird have real staying power?
]]>A Real-Life Problem
Last year I bought a house. For the first time in my life, I have an oil furnace, and a oil-fired storage water heating system. Last summer I was traumatized to hear the furnace starting in the middle of the day when it is 100 degrees Fahrenheit outside and I am not using any hot water.
I thought it would be a good idea to use the Simscape Thermal Liquid domain to analyze if it would be worth changing my current system to something else, for example a tankless system. To help me with that, my colleague Zack Peters offered to put together a few Simulink models.
The Model
I like using Simscape because you create a model by assembling components the way they are assembled in real life. In the model shown below, the blocks through which the water flows are colored in blue. You can find:
To model the heat going in and out of the tank, we used blocks from the Simscape Thermal Domain, highlighted in magenta in the model. You can find:
To test the simulation, we used the U.S. Department of energy residential water heater test procedures and standards for which consumer products should be tested against. One test, the 24-Hour Simulated Use Test determines the energy consumed by a water heater over a 24-hour period using a prescribed flow pattern. It requires that 64.3 gallons be drawn over the course of 6 separate draws separated 1 hour apart, followed by approximately 19 hours of no activity. It also specifies that a temperature set point of 135 F should be used.
When I compared the outlet temperatures of both systems, I could see that:
If we integrate the energy consumption, we can see that, obviously, the tankless water heater uses less energy. It transfers heat more efficiently and turns on only when hot water is needed. However there is more to tell to this story. The storage water heater is an oil fuel system while the tankless water heater runs off of propane, each of which have different costs per gallon in the market and different energy content stored per gallon.
Taking those numbers into account, we found out that the two systems use nearly the same number of gallons of fuel per day. When multiplied by the current cost of fuel here in New England, we find that the storage water heater costs $\$$2.18 per day while the tankless is only $\$$1.86 per day. That’s a difference of about $\$$116 per year.
Conclusion
I am pretty impressed to see that with just a few Simscape blocks we have been able to obtain numbers very close to the bills I received from my oil provider last summer.
If needed, we could increase the complexity of each component, but for now I think I have enough information to choose the system I need for my house. Now that I have this framework, I will be able to analyze different systems and compare those results to see if the cost/benefits of changing my system is worth it.
Now it's your turn
Download Zack's models for the Storage and Tankless systems, and let us know what you think of the Simscape Thermal Liquid domain by leaving a comment here.
]]>"Is there anything else you would like to tell us?"
This was the question whose answer led me to revisit an old friend of mine, the roifill function. This is a function that smoothly fills a specified region of interest in a gray-scale image. I implemented it in the early 1990s following a hallway conversation with another MathWorks developer, who said to me something along the lines of "I bet you could fill an image region based on the equation for a soap film surface along a wire rim."
(I don't know when the term inpainting became common, but I certainly had not heard it way back then.)
Here's an example of roifill in action.
I = imread('rice.png');
imshow(I)
If you now call roifill, you can select a polyline around any object you like, and the function will remove the object and smoothly fill in the resulting hole. Here's a sample output (can you find the missing grain of rice?):
Anyway, fast-forward to just two weeks ago. A customer that I know was helping us with a usability test for a planned MATLAB feature. I was sitting in the observation room. After the test tasks were completed, the user experience expert asked a common question: Is there anything else you'd like to tell us about?
It turns out that this customer knows a great deal about the Image Processing Toolbox, and he had some interesting feedback about the roifill function and the mask image.
roifill, you see, takes a "mask" image that specifies the set of pixel values that are to be filled. At least, that's the way I should have designed it. Unfortunately, I didn't, and my design decision was causing a lot of confusion. That's because roifill actually only replaces the interior pixels of the mask.
Let me show you a diagram to explain what I mean. Suppose your fill mask looked like this:
It would be a very reasonable expectation to think that roifill would replace all the pixels corresponding to this mask. And that's what our intrepid customer thought should happen. But it actually replaces only the interior pixels of the mask. That is, it replaces the pixels marked in gray:
The pixels on the edge of the mask are used to establish boundary conditions for the fill equation. These pixels stay the same. Hence the confusion.
As I walked back to my office after the usability test, I tried to figure out why I would have originally designed the function this way. I think there were two reasons. First, I was much less experienced with designing image processing functions then. Second, I'm pretty sure I was only thinking about an interactive workflow. That is, a user interactively draws a region around an object using a mouse, and then roifill processes the result.
In the interactive workflow, it doesn't matter (very much) how the mask border pixels are handled. However, in a noninteractive workflow, I now think it really only makes sense for the mask to represent the set of "bad pixels" that all need to be replaced by the function.
Ah well. Lessons learned. I will suggest that the Image Processing Toolbox development team think about how to address this behavior quirk in a future release.
In the meantime, dear reader ... Is there anything else you would like to tell us?
Let us know by leaving a comment.
Get
the MATLAB code
Published with MATLAB® R2014a
We are hosting another worldwide conference this month, on-line. Its coming Marcn 26, and has special tracks including
You can register on-line here.
There are two timezones for the conference, one for the Americas (with four sessions in Spanish), and one for India/Europe.
Please be sure to attend Roy Lurie's keynote entitled "Embracing Technical Computing Trends with MATLAB: Accelerating the Pace of Engineering and Science".
I will be attending at least part of the conference myself. Hope to see you there.
Get
the MATLAB code
Published with MATLAB® R2014a
High resolution measurements of the depth of the mold for the United States one cent coin provide an interesting data set.
During a visit around 1985 to what was then called the National Bureau of Standards, I was shown an instrument that made high resolution measurements of the thickness of small objects. The instrument would pass a stylus across the object multiple times, recording its height or depth, something like a square phonograph player. My hosts demonstrated the instrument by measuring the depth of a mold obtained from the United States Mint to make a one cent piece, a penny.
The result is a 512-by-512 array of integers in the range from 0 to 255. Back then, graphics on an array that large was beyond the capacity of MATLAB, so I subsampled the data to produce the 128-by-128 array called penny.mat that has been distributed in our demos directory for many years.
Now we can handle the full 512-by-512 array. Here it is, saved as uint8's. If you want to follow along with this blog on your own machine, click on this link, download the array, and convert it to doubles.
load penny512.mat
P = double(flipud(P));
The copper colormap is just a linear ramp of rgb values up to a copper color.
set(gcf,'renderer','zbuffer') colormap(copper) rgbplot(copper) title('copper')
Our first view of the data is a surface plot, with lighting and shading. You have to realize this is not a photograph of a penny. It is a top down view of a three-dimensional surface plot of the depth data obtained from the scanning stylus. The placement of the light and calculation of its reflections are critical.
surf(P); shading interp material metal lighting gouraud daspect([1,1,20]) axis tight axis off set(gca,'zlimmode','auto','climmode','manual'); light('pos',[1,1,2000],'style','inf'); view(2)
Zoom in on the date. You can see the striations resulting from the scanning action.
surf(P(50:130,360:490)); shading interp material metal lighting gouraud daspect([1,1,20]) axis tight axis off set(gca,'zlimmode','auto','climmode','manual'); light('pos',[1,1,2000],'style','inf'); view(2)
Now display a contour plot with 20 copper colored contour levels.
contour(P,20) axis square axis off
Here is a closer three dimensional view with tiny points on the surface.
surf(P,'facecolor','none','edgecolor','none',... 'marker','.','markeredgecolor','flat','markersize',1); daspect([1,1,20]) axis tight axis off set(gca,'zlimmode','auto'); camzoom(3);
Thanks to Eric Ludlam and Mike Garrity for help with this post.
Get
the MATLAB code
Published with MATLAB® R2014a
I'd like to introduce this week's guest blogger Alan Weiss. Alan writes documentation for mathematical toolboxes here at MathWorks.
Hi, folks. While I have not written a blog post for Loren before, if you use Optimization Toolbox™ or Global Optimization Toolbox then you might have read my work.
I am excited to describe how to use a new solver. Beginning with Release 2014a, Optimization Toolbox has mixed-integer linear programming. The intlinprog solver attempts to solve problems with linear objective functions, linear constraints, and (this is the new part) the constraint that some variables must be integer-valued.
The addition of integer constraints opens a surprisingly wide field of problems that are now solvable with Optimization Toolbox. There are examples in the documentation of a complicated factory production problem and solution, a travelling salesman problem, and a Sudoku puzzle solver.
This article starts with an explanation of how the Sudoku puzzle solver works. Then it extends the method to solving Hyper Sudoku puzzles. You might not be familiar with Hyper Sudoku; I will describe it presently.
Let's look at the Sudoku puzzle solver. As you probably know, a Sudoku puzzle is to fill a 9-by-9 grid with integers from 1 through 9 so that each integer appears only once in each row, column, and major 3-by-3 square, as in the following diagram, which includes puzzle clues.
The clues are encoded as follows: each row of a clue matrix has the form [i,j,k]. This means row i, column j, has clue k. If the first row is [1,2,2], then the clue at row 1, column 2 is equal to 2. Both the sudokuEngine and drawSudoku functions accept clues in this form, or the more conventional 9-by-9 matrix form.
B = [1,2,2; 1,5,3; 1,8,4; 2,1,6; 2,9,3; 3,3,4; 3,7,5; 4,4,8; 4,6,6; 5,1,8; 5,5,1; 5,9,6; 6,4,7; 6,6,5; 7,3,7; 7,7,6; 8,1,4; 8,9,8; 9,2,3; 9,5,4; 9,8,2]; drawSudoku(B)
The documentation shows how to solve this Sudoku puzzle using the sudokuEngine function:
drawSudoku(sudokuEngine(B));
LP: Optimal objective value is 29565.000000. Cut Generation: Applied 1 strong CG cut, and 2 zero-half cuts. Lower bound is 29565.000000. Relative gap is 0.00%. Optimal solution found. Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value; options.TolGapAbs = 0 (the default value). The intcon variables are integer within tolerance, options.TolInteger = 1e-05 (the default value).
How does sudokuEngine solve Sudoku puzzles? It simply encodes the rules and clues in a binary integer programming problem. The intlinprog function does the work. Isn't that great? You don't have to figure out a solution algorithm. Compare this with several Sudoku solvers in the File Exchange and Cleve's article about how to program a Sudoku solver. Those all use ingenuity to develop a Sudoku-solving algorithm. But we don't have to figure out an algorithm. We just have to figure out how to represent the rules of Sudoku as an integer linear programming problem. sudokuEngine uses a binary integer programming formulation.
The key trick in turning the rules of Sudoku into a binary integer program is to turn the 9-by-9 grid of clues and answers into a 9-by-9-by-9 cube of binary numbers. Think of the cube as a stack of 9-by-9 grids, one at level 1, one at level 2, up to a grid at level 9. A clue of, say, 5 in the original (3,4) grid becomes a 1 at level 5 in the cube at grid coordinate (3,4). In other words, the (i,j) grid can be viewed as a stack of 9 levels, and a value k in the (i,j) grid becomes a 1 at stack level k.
RGB = imread('sudokuGrid3D.jpg'); image(RGB); axis image; axis off;
Of course, for this to make sense, the solution needs to have exactly one 1 in the stack at each grid coordinate. So if x(i,j,k) represents the value of the (i,j) grid at level k, then we must have
$$\sum_{k=1}^9 x(i,j,k) = 1$$
at a solution.
What are the other rules of Sudoku? Each row has exactly one number of each value from 1 to 9. In equations, for each $j$ and $k$, we have
$$\sum_{i=1}^9 x(i,j,k) = 1$$
Each column has exactly one number of each value, too. So for each $i$ and $k$,
$$\sum_{j=1}^9 x(i,j,k) = 1$$
Each major 3-by-3 grid has exactly one number of each type (1 through 9). For the grid elements $1\le i\le 3$ and $1\le j\le 3$, and for each $1\le k\le 9$,
$$\sum_{i=1}^3\sum_{j=1}^3 x(i,j,k) = 1.$$
To represent all nine major grids, just add 3 or 6 to each $i$ and $j$ index:
$$\sum_{i=1}^3\sum_{j=1}^3 x(i+U,j+V,k) = 1,$$
where $U,V~\epsilon~\{0,3,6\}.$
The intlinprog function has a similar syntax to most Optimization Toolbox functions. The constraints on the solution x are bounds and linear equalities. All the entries in x are binary numbers, so all the lower bounds are 0 and all the upper bounds are 1.
lb = zeros(9,9,9); ub = 1 + lb;
The other constraints are all linear equalities. Optimization Toolbox solvers use the syntax
$$Aeq\cdot x = beq$$
to represent linear equalities, where x is regarded as a column vector.
x is naturally a 9-by-9-by-9 matrix. However, you can transform it to a column vector using the colon operator xcolumn = x(:). You can recover the 9-by-9-by-9 matrix representation using the reshape function: x = reshape(xcolumn,9,9,9).
Let's see how to express the first equation,
$$\sum_{k=1}^9 x(i,j,k) = 1,$$
as a matrix Aeq * x = beq. Each row of the Aeq matrix is matrix multiplied by the x matrix. This is the same as saying that entry i of beq is equal to the dot product of row i of Aeq with the vector x. So to determine the entries of row i of Aeq, make a 9-by-9-by-9 array of the multipliers of x in the equation, then turn the array into a row vector. For i and j ranging from 1 to 9, we have
temp = lb; % A zero array, 9-by-9-by-9 temp(i,j,:) = 1; % Set the coefficients of x to 1 % Then take temp(:)' as a row of Aeq
N = 9^3; % number of independent variables in x, a 9-by-9-by-9 array M = 4*9^2; % number of constraints, see the construction of Aeq Aeq = zeros(M,N); % allocate equality constraint matrix Aeq*x = beq beq = ones(M,1); % allocate constant vector of ones, beq counter = 1; for i = 1:9 % one in each depth for j = 1:9 temp = lb; % Clear temp temp(i,j,1:end) = 1; % Set each depth entry to 1 Aeq(counter,:) = temp(:)'; % Change the temp array to a row, put it in Aeq counter = counter + 1; end end
Now for the equation representing each row,
$$\sum_{i=1}^9 x(i,j,k) = 1$$
do a similar mapping from array to row of Aeq
for j = 1:9 % one in each row for k = 1:9 temp = lb; % clear temp temp(1:end,j,k) = 1; % one row in Aeq*x = beq Aeq(counter,:) = temp(:)'; % put temp in a row of Aeq counter = counter + 1; end end
Similarly, put in the column constraints.
for i = 1:9 % one in each column for k = 1:9 temp = lb; temp(i,1:end,k) = 1; Aeq(counter,:) = temp(:)'; counter = counter + 1; end end
Finally, set the coefficients for the major squares.
for U = 0:3:6 % one in each square for V = 0:3:6 for k = 1:9 temp = lb; temp(U+(1:3),V+(1:3),k) = 1; Aeq(counter,:) = temp(:)'; counter = counter + 1; end end end
Include the initial clues in the lb array by setting corresponding entries to 1. This forces the solution to have x(i,j,k) = 1.
for i = 1:size(B,1) lb(B(i,1),B(i,2),B(i,3)) = 1; end
The Aeq matrix, beq vector, and bounds lb and ub represent the rules of Sudoku and the clues. The objective function vector f is not important in this case; we just want a solution to the puzzle. We have to tell the solver that all variables are integers.
opts = optimoptions('intlinprog','Display','off'); % no exit messages intcon = 1:N; % all variables are integers f = []; % no objective function, we just want a feasible solution
Now call intlinprog to solve the puzzle.
[x,~,eflag] = intlinprog(f,intcon,[],[],Aeq,beq,lb,ub,opts);
To convert the solution to usable form, simply add up the numbers at each $(i,j)$ entry, multiplied by the depth at which the numbers appear:
if eflag > 0 % good solution x = reshape(x,9,9,9); % change back to a 9-by-9-by-9 array x = round(x); % clean up non-integer solutions y = ones(size(x)); for k = 2:9 y(:,:,k) = k; % multiplier for each depth k end S = x.*y; % multiply each entry by its depth S = sum(S,3); % S is 9-by-9 and holds the solved puzzle else S = []; end drawSudoku(S)
Hyper Sudoku differs from Sudoku by having more constraints. There are four 3-by-3 squares in addition to the major 3-by-3 squares that also require exactly one entry of each numeral from 1 through 9. The following image shows a Hyper Sudoku puzzle taken from Wikipedia
RGB = imread('hyperBoard.png'); % Colored squares are the new constraint regions image(RGB); axis image; axis off;
It is very easy to write the rules for Hyper Sudoku based on our work to date. There are four new constraint regions, each requiring 9 equations. So allocate a new matrix Anew, and fill it in.
Anew = zeros(4*9,N); % allocate 4*9 new rows beq = [beq;ones(4*9,1)]; % extended beq matrix
The constraints involve the 3-by-3 matrices with row and column indices 2 through 4, and also indices 6 through 8.
counter = 1; % reset counter lb = zeros(9,9,9); % reset lb to all zeros for U = [2,6] % start at 2 and 6 for V = [2,6] % for both rows and columns for k = 1:9 temp = lb; temp(U+(0:2),V+(0:2),k) = 1; Anew(counter,:) = temp(:)'; counter = counter + 1; end end end Aeq = [Aeq;Anew]; % Append Anew to Aeq
The Hyper Sudoku puzzle pictured above can be represented by the following matrix.
B2 = [1 8 1 2 3 2 2 8 3 2 9 4 3 5 5 3 6 1 4 6 6 4 7 5 5 2 7 5 4 3 5 8 8 6 3 3 7 5 8 8 1 5 8 2 8 8 7 9 9 1 6 9 2 9];
Include the clues in the lower bound.
for i = 1:size(B2,1) lb(B2(i,1),B2(i,2),B2(i,3)) = 1; end
Solve the puzzle.
[x,~,eflag] = intlinprog(f,intcon,[],[],Aeq,beq,lb,ub,opts); if eflag > 0 % good solution x = reshape(x,9,9,9); % change back to a 9-by-9-by-9 array x = round(x); % clean up non-integer solutions y = ones(size(x)); for k = 2:9 y(:,:,k) = k; % multiplier for each depth k end S = x.*y; % multiply each entry by its depth S = sum(S,3); % S is 9-by-9 and holds the solved puzzle else S = []; end % drawSudoku(S) % this would draw the board. But for a Hyper Sudoku figure: RGB = imread('hyperBoardSolved.png'); image(RGB); axis image; axis off;
The drawSudoku and sudokuEngine utilities ship with R2014a, so you can solve your own puzzles easily. However, perhaps not as easily as you might like. I also submitted an entry to the File Exchange called visualSudoku that is a visual editor for Sudoku puzzles. visualSudoku has no real code, it relies on sudokuEngine to solve puzzles, but it makes it easy to solve, save, and import Sudoku puzzles.
I also uploaded a visual Hyper Sudoku solver, and an engine for solving Hyper Sudoku puzzles.
The constraint matrices for Sudoku and Hyper Sudoku are quite sparse. I did not use that sparsity in the formulation. For examples using sparse constraint matrices, see the factory example or the travelling salesman problem in the documentation.
Did you enjoy seeing how to solve Sudoku or Hyper Sudoku puzzles with integer programming? Do you have any other types of puzzles you solve this way? Tell us about it here.
Get
the MATLAB code
Published with MATLAB® R2014a
Brett's Pick this week is more of a challenge than a "Pick".
If I were to search the MATLAB Central File Exchange for "face detection" (with the quotation marks) I would get a dazzling--and somewhat overwhelming--array of 44 hits. Trying to detect faces (or anything else*) in images seems to me a reasonable thing to want to do, and in my mind typifies the challenges that the Computer Vision System Toolbox was made to address. In fact, I shudder to think of what classical "image processing" code designed to detect faces might look like.
But how do we navigate the 44 files on the File Exchange that facilitate that detection process? Where do we start? I can sort the results by number of downloads, by rating, by relevancy, ....
Some of those downloads seem to have great potential. But where does one start?
For this post, I'm going to go "out-of-model" and focus not on a submitted File Exchange file, but on core capabilities of the Computer Vision System Toolbox. This Toolbox has a "Cascade Object Detector" that facilitates out-of-the-box detection of faces, eyes, noses, mouths, and upper bodies. In fact, there are two in-product face detectors, and "FrontalFaceCART" (frontal face detection using "Classification And Regression Tree" analysis) is the default tool behavior.
Using the detectors in the Computer Vision System Toolbox is stunningly simple; we simply instantiate an instance of the cascade object detector and, for visualization purposes, a vision.ShapeInserter:
faceDetector = vision.CascadeObjectDetector; shapeInserter = vision.ShapeInserter('BorderColor','Custom','CustomBorderColor',[0 255 255]);
We then read input image:
I = imread('visionteam.jpg');
imshow(I);shg;
Then we simply "step" the instantiated objects on the image of interest. The CascadeObjectDetector returns the bounding boxes of the detected objects, and the ShapeInserter delimits them in the image:
bbox = step(faceDetector, I); % Draw boxes around detected faces and display results I_faces = step(shapeInserter, I, int32(bbox)); imshow(I_faces), title('Detected faces');
Detecting faces using MathWorks' functionality is trivially easy. It requires MATLAB, the Image Processing Toolbox, and the Computer Vision System Toolbox, and a few readily reproducible lines of code.
Swag to the first person who shows me how to detect faces in "visionteam.jpg" with a File Exchange submission!
By the way, you can readily train MATLAB to find objects that are not in the "pre-trained" list above. For that, you might find these tools useful:
Get
the MATLAB code
Published with MATLAB® R2014a
Diagnostic Viewer
In R2014a, the Diagnostic Viewer has been completely re-engineered. This new interface can store the diagnostics from multiple past simulations, and offers filtering capabilities. Also, most warnings that were previously reported at the MATLAB prompt are now all routed to this central location.
Data Dictionary
Block parameters can now be resolved from a Data Dictionary file instead of the base workspace. This offers a lot of advantages like data grouping, change tracking, comparing and merging. Look at the documentation page What is a Data Dictionary? for more info.
Annotations with rich text, graphics, and hyperlinks
To help you document you model, we added the ability to include rich text, images, hyperlinks and tables to annotations.
Sliders, dials, and spinboxes available as parameter controls in masks
In the last few releases, a lot of functionality has been added to the Mask Editor. Here is an example of the kind of masks you can now create:
Git Source Control Support in Simulink Project
Git Source control is now available in Simulink Project. It allows you to create branches and merge them using the Model Comparison tool from Simulink Report Generator. It also allows you to push and fetch changes from a remote repository like github.
Now it's your turn
Let us know what is your favorite R2014a feature by leaving a comment here.
]]>At MathWorks, the software developers' sense of calendar time is SO out of sync with the real world. In the real world, MathWorks shipped R2014a, the first of our two annual releases, last week. In the developers' world, R2014a is old news. Right now all the developers are heads-down and sprinting toward the feature deadline for the next release, R2014b.
(On a personal, but somewhat related, side note, I recently reached my 20-year milestone with MathWorks. In my first decade with the company, MATLAB releases happened several years apart in unpredictable and somewhat chaotic fashion. The transition in 2004-2005 to a regular, repeatable, and completely reliable twice-per-year release schedule is one of the most impressive changes I've seen here.)
For image processing, the big news in MATLAB R2014a is webcam support. You can use MATLAB to bring in live images from USB Video Class webcams, which includes those built into laptops and those that plug into your computer via a USB port. Here's a video that explains how it all works.
Winning the award for "why didn't you do that years ago?" is multidimensional array support for flipud, fliplr, and rot90 functions. For image processing applications, that means you can use these functions directly on truecolor (M-by-N-by-3) images, like this:
rgb = imread('peppers.png'); subplot(1,3,1) imshow(rgb) title('Original') subplot(1,3,2) rgb2 = fliplr(rgb); imshow(rgb2) title('fliplr result') subplot(1,3,3) rgb3 = rot90(rgb); imshow(rgb3) title('rot90 result')
Those of you building small robots to take over the world might want to check out the new Raspberry Pi hardware support.
The MATLAB Math Team was busy making a collection of functions that make it easier to write linear algebra algorithms, including a bunch of is* functions (ishermitian, istril, etc.), an efficient Sylvester equation solver, and an option for computing left eigenvectors. The control system experts here seemed very excited about the Sylvester equation solver. I have no idea why, because in my senior year in EE I decided to specialize in signal processing instead of controls.
For the software developers among you, the still-new matlab.unittest framework is growing rapidly, with support for new custom plugins, the Test Anything Protocol (TAP), and automatic test parameterization. (I recently used the automatic test parameterization feature to update the test for the suite of examples in the book Digital Image Processing Using MATLAB.)
There are many other capabilities introduced in the new MATLAB version. I've mentioned here only a few of them that particularly caught my eye. I encourage you to take a look at the MATLAB release notes or watch the the R2014a Highlights video.
I don't think I've ever mentioned product licensing on this blog, but I wanted to give a shout to MATLAB Home, a new license option intended to support MATLAB for personal and hobby use.
Finally, Cleve was in my office yesterday and said I should write a blog about these MATLAB bottlecap refrigerator magnets that my son made for me. They are very popular. Just about everyone picks them up to play with them. I think I like the one with the white background better.
What's your vote?
Get
the MATLAB code
Published with MATLAB® R2014a
We got completely hooked on this challenge and before we knew it, we had a fully detailed simulator allowing you to play complete curling games!
We thought you might be interested, so we decided to share our Curling Simulator on MATLAB Central File Exchange. In this post we describe how we modeled the dynamics of the moving stones in Simulink.
The Big Picture
Our submission primarily consists of a MATLAB App and a Simulink model.
In the App, you decide where you want to aim, how fast you want to throw the stone, and what spin you want to put on it. When you click the Go button, the MATLAB App starts the Simulink model using the sim command. While it runs, the Simulink model uses a Level-2 MATLAB S-function to redraw the position and orientation of every stone in play.
In this post we are not going into the details of the interactions between the App and the model… maybe in a future post.
The Simulink implementation of the game dynamics is centered around one Second-Order Integrator block. The size of the signal is proportional to the number of stones on the ice.
Inside the Curling Model subsystem, we have two main components working in parallel: a contact model and a friction model.
The Friction Model
When he started developing a model for the friction, Corey quickly discovered that there is a lot of interest in the scientific community in curling. Understanding the movement of the stones on the ice has important implications in the field of tribology.
Based on what he learned, Corey parameterized lookup tables to calculate the force exerted by friction on the stone as a function of the stone's speed, angular velocity and the friction coefficient (which is function of how much you sweep…. Be sure you try the sweep button in the App!).
There are a lot of things about the behavior of spinning objects on ice that are counter-intuitive at first. For example, the faster you spin a curling stone, the less it actually curls. But if you don’t spin the stone at all, it erratically zig-zags down the ice. Researchers have yet to come to a consensus on the exact mechanisms for all of this behavior (you can read more about it in this article on National Geographics News or this article on phys.org).
Let's see the behavior in our simulator of two consecutive stones, one with low spin, and one with high spin
Watching a simulated stone curl realistically across our virtual ice was pretty exciting, but we couldn’t stop there. We expanded the model to simulate all 16 regulation stones at once. Doing that was actually pretty easy… all we had to do was drop the friction model for one stone into a For Each Subsystem.
The Contact Model
But what about simulating collisions? The contact model, we first computes the distance between each pair of stones. If we detect a contact, we compute new initial velocities and reset the Second-Order Integrator (This logic is similar to the bouncing ball demo).
To manage all the possible combinations of stones hitting each other, I used two nested For Iterator Subsystems (not to be confused with the For Each Subsystem) and an Enabled Subsystem. That way, the model computes new velocities only for the pair of stones in contact.
To compute the transfer of momentum, I found the perfect example in my Meriam and Kraige Engineering Mechanics Dynamics book I had bought during my first undergrad semester.
In the App, I was pretty impressed when I first saw the stones colliding together:
Now it’s your turn
Head on over to MATLAB Central, download the Curling Simulator, and use Simulink to hone your curling strategy! The 2018 Olympic qualifiers start sooner than you think!
]]>Sean's pick this week is Straightforward Copy & Paste by Yvan Lengwiler.
When presenting seminars on MATLAB, I often ask customers how to get data into and out of Excel. The answers are typically commands such as xlsread and xlswrite. These are great, and are often the best way to do this for an application. However, they can be overkill if you just need to move a single array to your active Excel sheet every now and then.
Let's see what Yvan has done:
% A simple magic square:
A = magic(5)
A = 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9
Because I left the semicolon off, it was echoed to the command line.
I'll copy it from there and paste it into Excel:
It pasted it in as a string into cells A1:A5. That's not what I wanted!
Instead, use Ivan's copy:
copy(A)
Now in Excel, when I paste, it'll paste the way I intended:
We could get that same functionality with just MATLAB by interactively opening the variable in the Variable Editor (or by using openvar) and then copying from there.
So where does this really shine? What if we want to extract the middle part of some text or a table from a website and paste it into MATLAB?
Since Yvan's copy and paste allow linefeeds and separators, it's easy to bring in complicated data like this.
Let's get a list of all of MathWorks' products from the website into MATLAB. First, highlight and copy the table for all products:
% Pause so I can go copy it, then press a key when ready pause % Use Yvan's Paste and show first few results products = paste(); disp(products(1:10))
[1x0 char] [1x0 char] 'MATLAB® Product Family' [1x0 char] 'MATLAB ' [1x0 char] 'Parallel Computing' 'Parallel Computing Toolbox ' 'MATLAB Distributed Computing Server ' [1x0 char]
This gives us the whole table including headers and empty elements. We want only the products. The product strings have a space at the end of them that we can use to determine if the non-empty cells are products.
% Define an anonymous function to test if each cell is a product % I.e. non-empty and the last element a space isproduct = @(x)~isempty(x) && isspace(x(end)); % Keep only the cells which are products products = products(cellfun(isproduct,products)); % Deblank it to remove those precious spaces at the end. Keep only % unique products since some show up a few times under different categories. products = unique(deblank(products)); disp(products)
'Aerospace Blockset' 'Aerospace Toolbox' 'Bioinformatics Toolbox' 'Communications System Toolbox' 'Computer Vision System Toolbox' 'Control System Toolbox' 'Curve Fitting Toolbox' 'DO Qualification Kit (for DO-178)' 'DSP System Toolbox' 'Data Acquisition Toolbox' 'Database Toolbox' 'Datafeed Toolbox' 'Econometrics Toolbox' 'Embedded Coder' 'Filter Design HDL Coder' 'Financial Instruments Toolbox' 'Financial Toolbox' 'Fixed-Point Designer' 'Fuzzy Logic Toolbox' 'Gauges Blockset' 'Global Optimization Toolbox' 'HDL Coder' 'HDL Verifier' 'IEC Certification Kit (for ISO 26262 and IEC 61508)' 'Image Acquisition Toolbox' 'Image Processing Toolbox' 'Instrument Control Toolbox' 'LTE System Toolbox' 'MATLAB' 'MATLAB Builder EX (for Microsoft Excel)' 'MATLAB Builder JA (for Java language)' 'MATLAB Builder NE (for Microsoft .NET Framework)' 'MATLAB Coder' 'MATLAB Compiler' 'MATLAB Distributed Computing Server' 'MATLAB Production Server' 'MATLAB Report Generator' 'Mapping Toolbox' 'Model Predictive Control Toolbox' 'Model-Based Calibration Toolbox' 'Neural Network Toolbox' 'OPC Toolbox' 'Optimization Toolbox' 'Parallel Computing Toolbox' 'Partial Differential Equation Toolbox' 'Phased Array System Toolbox' 'Polyspace Bug Finder' 'Polyspace Code Prover' 'RF Toolbox' 'Real-Time Windows Target' 'Robust Control Toolbox' 'Signal Processing Toolbox' 'SimBiology' 'SimDriveline' 'SimElectronics' 'SimEvents' 'SimHydraulics' 'SimMechanics' 'SimPowerSystems' 'SimRF' 'Simscape' 'Simulink' 'Simulink 3D Animation' 'Simulink Code Inspector' 'Simulink Coder' 'Simulink Control Design' 'Simulink Design Optimization' 'Simulink Design Verifier' 'Simulink PLC Coder' 'Simulink Real-Time' 'Simulink Report Generator' 'Simulink Verification and Validation' 'Spreadsheet Link EX (for Microsoft Excel)' 'Stateflow' 'Statistics Toolbox' 'Symbolic Math Toolbox' 'System Identification Toolbox' 'SystemTest' 'Trading Toolbox' 'Vehicle Network Toolbox' 'Wavelet Toolbox'
And there we have it, MathWorks' 82 products!
How many of those products have you heard of (or not heard of)? I'm proud to say that I've only not heard of two of them; I'll send some MathWorks swag to the first person to guess each correctly :) Limit: three guesses per day
Do you have complicated data that you need to get into MATLAB? If so, what are your favorite tools or workflows to do this?
Give this a try and let us know what you think here or leave a comment for Yvan.
Get
the MATLAB code
Published with MATLAB® R2013b
At the recent SIAM conference on Parallel Processing for Scientific Computing, I met my friend Tim Davis. Tim is a professor at the University of Florida, a consultant to MathWorks, and author of many of the functions involved in MATLAB sparse matrix calculations. One of his hobbies involves the use of sparse matrices to translate music into visual art.
Tim maintains the University of Florida Sparse Matrix Collection, which he and others use to test their sparse matrix algorithms. The collection includes web pages of beautiful visualizations of the sparse matrices created by a package called Graphviz.
A London design agency saw the images, and contacted Tim to see if music could be visualized. Tim hadn't done music before, but agreed to give it a try. The agency liked what they saw, and commissioned Tim to do the artwork for the London Electronic Arts Festival held last November. His artwork appeared as the theme artwork for LEAF 2013, appearing on posters and billboards around London last fall, and still appears as the background on the agency's web page today.
To produce this art, Tim created an algorithm in MATLAB which processes the music via the FFT and then translates the frequencies into a MATLAB sparse matrix. The matrix is then drawn via Graphviz, a force-directed graph visualization package. To generate the visualization, Graphviz places electrical charges on the nodes and springs on the edges, and then does a physics simulation to find a low-energy state. The color reflects the notes, where low frequencies are blue and high ones are red, using the MATLAB jet colormap.
Link to gallery image. Link to music.
Metastaseis is a complex orchestral piece by Iannis Xenakis, a composer and architect who translated his architectural thoughts of space and structure into music, relying heavily on mathematics and algorithms. In Xenakis' preliminary sketches for Metastaseis, lines and their intersection in space form hyperbolic paraboloids, which he translated into musical notes and chords. Xenakis described his process as "vision -> rules -> works of art." He describes an inverse path as "rules -> vision." Tim's visualization completes Xenakis' inverse path by using algorithms to create visual artwork from Xenakis' music.
Link to gallery image. Link to music.
Blue Monday is by New Order, a British rock band that combines new wave and electronic dance music. Tim's rules do not specify "build a mesh here." Instead, mesh-like structures arise because of the repetitive bass notes that appear in electronic music.
Link to gallery image and music.
Electrodoodle is a piece of electronic music by Kevin Macleod. Like Blue Monday, this piece of music has lots of deep repetitive bass notes.
Link to gallery image. Link to music.
Morning Has Broken is a hymn written in 1931 by Eleanor Farjeon who set it to an old Gaelic tune. This version is arranged and performed by Cat Stevens. Tim created this piece for Iain and Di Duff. Iain's work appears in MATLAB as the algorithms underlying LDL and backslash when A is sparse and symmetric indefinite.
Link to gallery image and music.
Tim's visualization of Bach's Toccata and Fugue in D minor (BWV565) reflects the tensions in this complex orchestral piece.
All these images are reproduced here with Tim Davis's permission and may not be copied. For many more images and information about Tim's work, see NotesArtStudio.com.
Get
the MATLAB code
Published with MATLAB® R2014a
Jiro's pick this week is seconds2human by Rody Oldenhuis.
There was a time when I obsessed over waitbars. Why is that? I'm not sure, but Ned has a few guesses. But one of the features that I look for in these waitbar entries is the ability to show the remaining time, like some of these. I have created my own waitbar as well to include the time remaining, and I remember thinking quite a bit about how I display the time. I can easily calculate the time using tic and toc. But what's the best way to display the seconds? It's probably better to say "3 hours 44 minutes 13 seconds remaining" than to say "13453 seconds remaining". So I had to do some manipulations to have an appropriate display, and many of the entries I highlighted above do that as well.
Rody's seconds2human does this nicely. You simply pass in the time in seconds, and it will convert it to something human understandable. There's an approximate display:
seconds2human(13453)
ans = About 3 hours and 44 minutes.
and an accurate display:
seconds2human(13453,'full')
ans = 3 hours, 44 minutes, 13 seconds.
Now, I can just make use of this and not have to worry about decomposing the seconds.
hWait = waitbar(0, 'Calculating...'); total = 1000; % Start timer tstart = tic; for iter = 1:total SomeComputation(iter); % Calculate remaining time remaining = toc(tstart)*(total-iter)/iter; % Update waitbar waitbar(iter/total,hWait, ... ['Time remaining: ',seconds2human(remaining)]) end delete(hWait);
Comments
I also want to point out that Rody's code is very well-written. It is well-documented, has error checking, and allows for a vectorized call, i.e. you can pass an array of numbers. Give this a try and let us know what you think here or leave a comment for Rody.
Get
the MATLAB code
Published with MATLAB® R2013b
Let's see how to debug Simulink Coder executables in R2013b.
The Problem
If you want to debug an executable, you need debugging information. For example, when using Microsoft Visual Studio, this information is stored in a Program Database, a file with a .pdb extension.
For this .pdb file to be generated, you need to pass specific flags to the compiler.
For many releases, when I needed to debug a Simulink Coder executable, my workflow was to generate code , open the makefile, search for this line of code:
change the value to 1, and execute the .bat file to build the executable.
I have to admit, this process was not very efficient and had to be repeated every time your regenerated code, because the modified makefile would get overwritten. Of course, there were ways to automate some part of this workflow, but that required some pretty advanced knowledge of Simulink.
The R2013b workflow
In R2013b, enabling debugging can now be made directly from the Simulink interface. No more need to be a master hacker!
All you need to do is set the Build configuration option to Debug in the code generation section of the model configuration.
Build your model and a .pdb file will appear in your current directory along with the executable.
In Visual Studio, go to File -> Open -> Project/Solution and select the executable. Open the source code you want to debug, set breakpoints and click the Start Debugging button.
Now it's your turn
Let us know what you think of this new enhancement in R2013b by leaving a comment here.
]]>Arithmetic is associative, right? Well, in the world of paper and pencil, where you can often do calculations exactly, that can be true. However, in the computing world, where real numbers can't always be represented exactly because of working with finite precision datatypes, it turns out that you can't depend on the arithmetic to behave the way you were taught in grade school.
Suppose I want to check the following:
$$ \sqrt {2} = 2/\sqrt {2} $$
I can do this analytically using the Symbolic Math Toolbox.
symsqrt2 = sqrt(sym(2)); shouldBeZero = symsqrt2 - 2/symsqrt2
shouldBeZero = 0
Now let's perform the same calculation numerically.
mightBeZero = sqrt(2) - 2/sqrt(2)
mightBeZero = 2.2204e-16
What is happening here is that we are seeing the influence of the accuracy of floating point numbers and calculations with them. I discussed this in an earlier post as well.
Now let's try something a little different. First, let's find out what the value of eps is for %sqrt {2}%. This should be the smallest (in magnitude) floating point number which, when added to %sqrt {2}%, produces a number different than %sqrt {2}%.
sqrt2eps = eps(sqrt(2))
sqrt2eps = 2.2204e-16
Next, we want a number smaller in magnitude than this to play with. I'll use half its value.
halfsqrt2eps = sqrt2eps/2
halfsqrt2eps = 1.1102e-16
And now let's calculate the following expressions, symbolically and numerically.
$$expr1 = \sqrt{2} - \sqrt{2} + halfsqrt2eps$$
$$expr2 = (\sqrt{2} - \sqrt{2}) + halfsqrt2eps$$
$$expr3 = \sqrt{2} + (-\sqrt{2} + halfsqrt2eps)$$
First we do them all symbolically.
expr1 = symsqrt2 - symsqrt2 + sym(sqrt2eps)/2 expr2 = (symsqrt2 - symsqrt2) + sym(sqrt2eps)/2 expr3 = symsqrt2 + (-symsqrt2 + sym(sqrt2eps)/2) double(expr1)
expr1 = 1/9007199254740992 expr2 = 1/9007199254740992 expr3 = 1/9007199254740992 ans = 1.1102e-16
Symbolic results are all the same and return half the value of eps.
Now we'll calculate the same expressions numerically.
expr1 = sqrt(2) - sqrt(2) + halfsqrt2eps expr2 = (sqrt(2) - sqrt(2)) + halfsqrt2eps expr3 = sqrt(2) + (-sqrt(2) + halfsqrt2eps)
expr1 = 1.1102e-16 expr2 = 1.1102e-16 expr3 = 2.2204e-16
So what's going on here? As I stated earlier, this example illustrates that floating point arithmetic is not associative the way symbolic arithmetic is. There's on reason to get upset about this. But it is worth understanding. And it might well be worth rewriting a computation occasionally, especially if you are trying to compute a very small difference between two large numbers.
Have you found yourself in a situation where you needed to rewrite how to calculate a numeric result (like here, by different groupings) to ensure you got a more accurate solution. Let me know about it here.
Get
the MATLAB code
Published with MATLAB® R2013b
I was impressed to see how deep this assessment goes to identify the strengths and weaknesses of the methods, tools, practices, organizational structure, and environment used in your current MBD workflow... and of course at the end coming up with both immediate and long term improvement plans.
If your organization is using Model-Based Design for Production Real-Time Embedded Systems, you might be wondering how you benchmark against leaders in your industry. Maybe you need to know the best way to eliminate redundant processes or tools, and reduce duplication of efforts – and how to standardize these things across your entire company. Perhaps you now want to change your development process to meet a certification standard (like DO-178, ISO 26262 or IEC 61508). If you’re seeing yourself in any of those scenarios, or if you want to improve overall efficiency or check how you’re doing 3 or 5 years into your Model-Based Design efforts, a Process Assessment might be for you.
Model-Based Design Process Assessment and Maturity Framework
After working with multiple organizations in various domains, our team of experienced engineers identified six core competencies crucial to deploying a mature Model-Based Design development environment. These six pillars were then decomposed into Process Groups, and further broken down into over 200 Process Attributes:
After spending time in your organization, evaluating your processes and interviewing employees, our Consulting Services group will be able to tell you how your adoption of Model-Based Design compares with other companies in your industry and identify key areas that could be improved to add value to your organization. An implementation plan suggests a structured set of activities to achieve immediate and long term results.
This radar chart, with the Pillars as the main spokes, compares the maturity assessments between leaders (top 20%), average, and laggards (bottom 20%) in the automotive industry based on Process Assessments performed by MathWorks. Where does your organization fit in?
Now it's your turn
Look at the Model-Based Design Process Assessment and Maturity Framework to get more details on how it could help your organization.
Have you already used this service? Are you interested? Let us know by leaving a comment here.
]]>The Tektronix 4081 minicomputer at Argonne National Laboratory in the summers of 1977 and '78 played an important role in the early development of MATLAB.
Tektronix, Inc., is a company founded in Portland, Oregon in 1946 that was best known for its oscilloscopes. In 1976 Tektronix temporarily entered the computer business. They marketed a sophisticated display together with one of the first 32-bit minicomputers, the model 7/32 from Interdata, and their own operating system. The result was their model 4081.
The marketing material emphasized the display. It incorporated traditional storage tube graphics, where the tube phosphor provides the memory and which Tektronix had been using for years in its oscilloscopes, and the new refresh graphics, where digital memory drives the display. In this view of the system, the computer is simply the graphics controller. The Tektronix software was known as GOS, for Graphics Operating System.
Here is a screen shot of a lunar lander game on the machine. The background was generated by the storage tube graphics and remained fixed. The lander and the text at the bottom were raster graphics. The lander's rockets were controlled by arrow keys on the keyboard. Remember, this was 1976. Any computer with an interactive display being marketed at the time had to have a lunar lander. This one was fancier than most of the competition. Here is a link to a demo reel made at the time that shows the lander in action, and several other applications on the 4081.
I visited the Mathematics and Computer Science Division at Argonne National Laboratory, just west of Chicago, every summer in the late 1970s to work on Linpack with Jack Dongarra, Pete Stewart and Jim Bunch. We did most of our work on Argonne's central main frame computer, the IBM 370/195. This was a batch oriented system. We submitted jobs and got the output back some time later. The time lapse could be a few minutes or a few hours, depending upon how many other jobs there were. It was quicker at night and on weekends.
The jobs would be punched on IBM cards and the output printed on fan folded paper returned in bins at the computer center. Or, somewhat more conveniently, we would work at terminals, preparing files containing images of card decks, and have output returned in files viewable at the terminal.
One important goal of my first MATLAB was to make matrix computation interactive. I wanted to enter an expression involving matrices in a natural mathematical notation and have the result evaluated immediately. The first MATLAB was a Fortran program that used about a dozen subroutines from Linpack and Eispack to create a simple matrix calculator. The TSO time sharing facility on another IBM main frame, the 360/75, was reasonably interactive. As I remember it, if there were not too many other users on TSO, the response time was pretty good.
The MCS Division at Argonne acquired a Tektronix 4081 sometime before I arrived for my 1977 summer visit. I was excited. This was my first personal computer. There was a Fortran compiler and I could run the MATLAB that I was developing. There was just one difficulty. The 7/32 had only a 16-bit address space and was limited to 64K of memory -- that's 64 kilobytes. The designers of GOS knew that 64K was a problem and so had provided a swap feature. Fortran designates a section of global memory as COMMON. It was possible to retain the COMMON memory, swap all the executable code in the remainder of the memory, and keep the program running.
So, my MATLAB consisted of five or six blocks of code. The parser was one or two blocks, LU and QR were a block, SVD was a block, and so on. The matrices were kept in common and the various blocks swapped in and out as needed. The memory limitation turned out to be a blessing after all. This blocked structure of the software forced a valuable clarity of organization.
The 4081 was in a small office. I loved the feel of having the machine all to myself. It was the size of a desk, but there is not much difference between using a computer that is a desk and using a computer that is on a desk. The response time was good, and it was predictable and consistent. It would be four or five years before we would have the IBM PC and the Sun workstation, but this experience gave me a taste for what personal computing would be like.
Argonne's 4081 was connected to the ARPA-net, the predecessor of today's Internet. Jim Pool had been head of the MCS Division at Argonne and had been involved in the acquisition of the 4081, but had moved to the Department of Energy in Washington, DC before the machine was delivered. Jim was also on the ARPA-net. I sent my first ever email to Jim from the 4081. I told him that I was using the machine, I was working on this new program that I called MATLAB, and I included the entire text supporting the help command. The entire email was less than 700 lines.
Ned Thanhouser was Software Manager in charge of the Graphics Operating System at Tektronix. He came to visit Argonne one day to see what we were doing with the machine. Jack Dongarra got hold of me and the three of us ended up ordering pizza and staying far into the night. Ned had a briefcase full of cassette tapes and floppy discs with demos, games, and other goodies. They were 8-inch floppies in those days. The cassette tape reader could be troublesome because the tape alignment was delicate, but there was a thick Chicago telephone book on top of the unit and if you hit the reader with the phone book, it would relax its tolerances.
Ned later moved from Tektronix to Intel, where he spent most of his career. Ned's grandfather was a famous motion picture industry pioneer, and so after retiring from Intel, Ned has been involved with the preservation of early silent motion pictures, at Thanhouser Company Film Preservation.
Ned wrote the lunar lander and made the demo reel that I mentioned above. While preparing this blog, Ned and I had a couple of very pleasant phone conversations and email exchanges, remembering old times. He sent me this recollection by his colleague, Sam Malicot.
And, after all these years, I can reveal that Ned's visit to Argonne that day led to the anonymous quote on the Table of Contents page of the Linpack Users' Guide. Thanks, Ned.
Tektronix 4081 System Handles Both Refresh, Storage Graphics, Computer World, May 17,1976, page 29.
Ned Thanhouser, Intermixing Refresh and Direct View Storage Graphics, SigGraph, 1976.
Get
the MATLAB code
Published with MATLAB® R2014a
by Kent Millard
Let’s say you’re wondering how to plot some data contained in a 3D matrix. You’ve searched high and low to find an answer, but without satisfying results. Finally, you come across the MATLAB Answers page and there it is… the big blue “Ask a Question” button!
With gathering excitement, you press the button and begin to type.
This is where it gets interesting.
Before you even finish your query, similar questions that have already been answered start to appear. One of these might just have the information you need. If so, you get an instantaneous answer pulled from our database of more than 75,000 answered questions!
Here’s how it works.
Start typing your question—keywords or natural language work equally well. The most relevant, answered questions will appear below. The boxes with numbers tell you how many answers there are to a particular question. Green boxes mean the person asking the question confirmed that one of the answers solved their problem.
Questions with tiny blue icons (it’s a stylized MathWorks logo) have been posted by the MathWorks Support Team, and they always include an answer that solves the problem.
We think auto-suggest will save you time. But rather than read more about it, give it a try and tell us what you think!
]]>Sean's pick this week is WorkerObjWrapper by MathWorks' Parallel Computing Team.
MATLAB's Parallel Computing Toolbox provides you with the the ability to open a pool of MATLAB workers that you can distribute work to with high level commands like parfor.
When communicating with the workers in this pool, there will always be an overhead in data communication. The less data we can transmit to the workers the better speed improvements we'll see. This can be difficult when working with large arrays and can actually cause parallel computations to be slower than serial ones. WorkerObjWrapper has provided a convenient way to make data persist on a worker; this could be large arrays, connections to databases or other things that we need on each iteration of a parfor loop.
We're going to pull some financial data from Yahoo! using the connection from the Datafeed Toolbox.
I have a list of securities and the corresponding fields I want from them:
% Securities and desired fields securities = {'MAR','PG','MSFT','SAM',... 'TSLA','YHOO','CMG','AAL'}; fields = {'High',{'low','High'},'High','High',... {'Low','high'},'Low',{'low','Volume'},'Low'};
I first want to make sure there is an open parallel pool ( parpool ) to distribute computations to. I have a two core laptop, so I'll open two local workers by selecting the icon at the bottom left hand side of the desktop.
I've written three equivalent functions to pull the prices from Yahoo!
First, a sanity check to make sure they all do the same thing:
ff = fetchFOR(securities,fields);
fp = fetchPARFOR(securities,fields);
fw = fetchWOWPARFOR(securities,fields);
assert(isequal(ff,fp,fw)); % Errors if they're not equal
Since the assertion passed, meaning the functions return the same result, we can now do the timings. I'll use timeit.
t = zeros(3,1); % Measure timings t(1) = timeit(@()fetchFOR(securities,fields),1); t(2) = timeit(@()fetchPARFOR(securities,fields),1); t(3) = timeit(@()fetchWOWPARFOR(securities,fields),1); % Show results fprintf('%.3fs %s\n',t(1),'for',t(2),'parfor',t(3),'parfor with WorkerObjWrapper')
8.631s for 5.991s parfor 4.255s parfor with WorkerObjWrapper
So we can see that creating the connection once on each worker in the parallel pool and then using parfor gives us the best computation time.
Do you have to work with large data or repeat a process multiple times where parallel computing might help? I'm curious to hear your experiences and the challenges that you've faced.
Give it a try and let us know what you think here or leave a comment for our Parallel Computing Team.
Get
the MATLAB code
Published with MATLAB® R2013b
In our recent post, Mike Hosea and I talked about adjusting both the absolute and relative tolerances for getting more accurate results when calculating a double integral. Today we'd like to talk about choosing the method of integration as well as which order to choose for the first dimension (direction) of integration.
integral2 has, at present, two different integration methods, 'tiled' and 'iterated', not counting the 'auto' method that chooses between them.
Each integration method employs a type of "divide-and-conquer" approach to double integration but in very different ways. The 'tiled' method is based on quad2d's approach of dividing the region into quadrants and approximating the integral over each quadrant by a 2-D quadrature rule. If the error condition on a rectangle is not met, the rectangle is divided into quadrants and so forth. The 'iterated' method, as the name suggests, performs the integration as iterated one-dimensional integrals, so its "divide-and-conquer" strategy is to tackle one dimension first and then move on to the other.
So, why have two methods? Each has different strengths. Usually the tiled approach is the faster of the two, but it has its weaknesses. As an example of this, consider the way integrating over non-rectangular regions was formerly done in MATLAB. Because dblquad would only integrate over rectangles, one would "mask" the input function to make it zero inside of a rectangular region. For example, to integrate F over a triangle with the vertices (0,0),(1,0), and (1,1), you would write
F = @(x,y)exp(10*x).*tan(y) Q = dblquad(@(x,y)F(x,y).*(y <= x),0,1,0,1)
F = @(x,y)exp(10*x).*tan(y) Q = 1072.6
Now you can instead write
Q = integral2(F,0,1,0,@(x)x)
Q = 1072.6
This is more than a convenience. Telling the integrator where the boundary is gives it an enormous advantage. Conversely, masking the integrand has devastating effects on performance and accuracy. The underlying quadrature rule gives the best results when the integrand is smooth, and a discontinuity where the integrand drops precipitously to zero will cause the adaptive quadrature algorithm to subdivide and further subdivide in that area until the error estimates across the discontinuity are sufficiently small.
This is especially problematic for quad2d and integral2 with the 'tiled' method. With the 'iterated' method, any curvilinear discontinuity is covered by a single interval in any given integration, but in two dimensions, subdividing each tile into 4 pieces, the integrator ends up covering a curvilinear path with a large number of small rectangles. Usually the limit on the number of function evaluations will be reached before it can polish them all off.
We can visualize this using quad2d's FailurePlot option and use the MaxFunEvals option to control how much work is done before the code gives up. Suppose
F = @(x,y)sin(x./(y.*y + 1)).*(y <=x)
F = @(x,y)sin(x./(y.*y+1)).*(y<=x)
Then
Q = quad2d(F,0,1,0,1,'MaxFunEvals',60,'FailurePlot',true)
Warning: Reached the maximum number of function evaluations (60). The result fails the global error test. Q = 0.26562
There's nothing special about the larger rectangle on the bottom. quad2d preferentially attacks the rectangles with the largest error, so when a larger rectangle just barely fails the error test, it hangs around if the integrator has bigger problems elsewhere. Those bigger problems are signified by smaller rectangles. Now, if we increase the number of allowed function evaluations, we see that quad2d is still struggling with the boundary.
Q = quad2d(F,0,1,0,1,'MaxFunEvals',600,'FailurePlot',true)
Warning: Reached the maximum number of function evaluations (600). The result fails the global error test. Q = 0.26518
Since integral2's 'tiled' method is based on quad2d's algorithm, it has the same issue.
Q = integral2(F,0,1,0,1)
Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test. Q = 0.26514
Although it takes some time to do it, the 'iterated' method can deal with the problem.
tic Q = integral2(F,0,1,0,1,'method','iterated') tq = toc
Q = 0.26513 tq = 5.0793
If you find that the 'iterated' method seems slow, consider loosening the tolerances. integral2 uses the same tolerances for all subintegrations, and it turns out that this can be a very conservative thing to do. With the default tolerances, the previous example took over 4 seconds on the machine we publish the blog from.
tic Q = integral2(F,0,1,0,1,'method','iterated','AbsTol',1e-6,'RelTol',1e-3) tq(2) = toc
Q = 0.26514 tq = 5.0793 0.14705
And that took less than a tenth of second on the same machine.
Discontinuities may snake through the region of integration, but occasionally they might be linear and parallel to one axis or the other. In that case the order of integration matters when using the 'iterated' method. Since every outer integration performs many inner integrations, it tends to be more useful to make the outer integrations the more efficient ones if possible, and this means letting the inner integrations take care of jump discontinuities. In this example the integrand function alternates between z = 0 and z = y depending on whether floor(x) is even or odd.
f = @(x,y)mod(floor(x),2).*y tic integral2(f,0,20,0,1,'method','iterated') toc
f = @(x,y)mod(floor(x),2).*y ans = 5 Elapsed time is 5.625905 seconds.
That worked, but the "inner" integration is in the y-direction, and for a given x that means that the function is either identically zero or a straight line. Either is very easy to integrate. However, this means that the integrand for the outer integral has jump discontinuities from 0 to 0.5 depending on whether floor(x) is even or odd. We'd rather have the inner integrals deal with any (or at least most of) the discontinuities. Changing the order of integration is easy.
tic integral2(@(y,x)f(x,y),0,1,0,20,'method','iterated') toc
ans = 5 Elapsed time is 0.264443 seconds.
Unlike dblquad, integral2 can integrate singularities on the boundary if they are not too severe. If your problem has a singularity in the interior of the region of integration, you should divide the integral into pieces so that all singularities reside on the boundaries of regions of integration. For example,
f = @(x,y)1./sqrt(abs(x.*y))
f = @(x,y)1./sqrt(abs(x.*y))
If we integrate this over the square -1 <= x <= 1, -1 <= y <= 1, we have a problem because f is singular along both x and y axes. Obviously f is symmetric with respect to both axes, so we need only integrate over 0 <= x <= 1, 0 <= y <= 1 and multiply by 4, but let's set this aside and try to do it directly.
integral2(f,-1,1,-1,1,'RelTol',1e-12,'AbsTol',1e-14)
Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test. ans = 15.869
The 'iterated' method doesn't bail us out this time.
integral2(f,-1,1,-1,1,'RelTol',1e-12,'AbsTol',1e-14,'Method','iterated')
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 9.0e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy. Warning: The integration was unsuccessful. ans = NaN
However if we split the integration into four pieces, putting the singularity on the boundary, the problem is absolutely routine.
Q = 0; Q = Q + integral2(f,-1,0,-1,0,'RelTol',1e-12,'AbsTol',1e-14); Q = Q + integral2(f,-1,0,0,1,'RelTol',1e-12,'AbsTol',1e-14); Q = Q + integral2(f,0,1,-1,0,'RelTol',1e-12,'AbsTol',1e-14); Q = Q + integral2(f,0,1,0,1,'RelTol',1e-12,'AbsTol',1e-14)
Q = 16
What kinds of functions do you need to integrate? Will the flexibility from integral2 and it's other companions help you? Let us know here.
Get
the MATLAB code
Published with MATLAB® R2013b
Today I was talking to Spandan (a regular guest blogger here) about the history of the Image Processing Toolbox function regionprops. This function is a widely-used work-horse of the toolbox. Many of my blog posts here have shown it in practice.
The function regionprops computes properties of regions of pixels. Most often, these regions correspond to individual "objects" (connected components of foreground pixels) of a binary image. Here's an example.
bw = imread('text.png');
imshow(bw)
s = regionprops(bw)
s = 88x1 struct array with fields: Area Centroid BoundingBox
The output is an 88x1 struct array, elements of which contain measurements for objects in the binary image. (We didn't tell regionprops what measurements to compute, so it computed Area, Centroid, and BoundingBox by default.) Here are the measurements for the 10th object.
s(10)
ans = Area: 61 Centroid: [54 39] BoundingBox: [49.5000 33.5000 9 11]
This form of output is convenient in some ways but inconvenient or awkward in others. For example, it requires the obscure "comma-separated-list" syntax to get all the area measurements into a single vector to call a function such as mean, min, or max.
a = [s.Area]; mean(a)
ans = 54.2386
min(a)
ans = 6
max(a)
ans = 106
In the comma-separated list syntax for struct arrays, the expression [s.Area] automatically expands into the expression [s(1).Area, s(2).Area, ..., s(end).Area]. But relatively few people recognize or understand this syntax.
Another inconvenient aspect of the struct array aspect is that it's hard to get a tabular view of all the data. As shown above, you can see the data for the 10th object by displaying s(10), but it's hard to see all the data for multiple objects. For example:
s(1:5)
ans = 5x1 struct array with fields: Area Centroid BoundingBox
If I were designing regionprops from scratch today, based on the most recent version of MATLAB, I would have regionprops return a table. Tables were introduced as a new fundamental MATLAB type in R2013b. Here's a simple example to illustrate.
age = [38 43 38 40 49]'; height = [71 69 64 67 64]'; weight = [176 163 131 133 119]'; blood_pressure = [124 93; 109 77; 125 83; 117 75; 122 80]; t = table(age,height,weight,blood_pressure)
t = age height weight blood_pressure ___ ______ ______ _______________ 38 71 176 124 93 43 69 163 109 77 38 64 131 125 83 40 67 133 117 75 49 64 119 122 80
t.height is an ordinary 5-by-1 vector of heights.
t.height
ans = 71 69 64 67 64
t.blood_pressure is an ordinary 5-by-2 matrix. (Table "variables" are often column vectors, but they can have any size as long as the number of rows is the same as the height of the table.)
t.blood_pressure
ans = 124 93 109 77 125 83 117 75 122 80
We are excited about the introduction of tables in MATLAB, and we think they will gradually be used in more and more situations where previously people used struct arrays.
Such as regionprops.
So what to do?
Well, the function struct2table comes in very handy here! As you might have guessed from the name, it automatically converts a struct array into a table.
Try it:
t = struct2table(s)
t = Area Centroid BoundingBox ____ __________________ ____________ 66 11 13.5 [1x4 double] 41 7.6829 38.171 [1x4 double] 63 17.889 38.857 [1x4 double] 80 22.95 15.963 [1x4 double] 53 26.472 36.604 [1x4 double] 63 34.889 16.857 [1x4 double] 63 34.889 38.857 [1x4 double] 41 43.683 38.171 [1x4 double] 48 48.396 15.875 [1x4 double] 61 54 39 [1x4 double] 63 56.889 16.857 [1x4 double] 41 65.683 16.171 [1x4 double] 48 67.396 37.875 [1x4 double] 105 78.79 16.79 [1x4 double] 66 76.5 39 [1x4 double] 68 93.529 39.265 [1x4 double] 92 100 16.696 [1x4 double] 41 106.68 38.171 [1x4 double] 68 113.53 17.265 [1x4 double] 6 114 31.5 [1x4 double] 33 114 39 [1x4 double] 85 123.24 38.259 [1x4 double] 48 121.4 15.875 [1x4 double] 63 129.89 16.857 [1x4 double] 100 134.87 40.86 [1x4 double] 41 138.68 16.171 [1x4 double] 63 145.89 38.857 [1x4 double] 61 149 17 [1x4 double] 80 159.95 15.963 [1x4 double] 85 163.26 185.76 [1x4 double] 6 156.5 195 [1x4 double] 6 156.5 211 [1x4 double] 85 163.26 218.76 [1x4 double] 48 159.4 37.875 [1x4 double] 61 164 107 [1x4 double] 68 164.26 117.47 [1x4 double] 63 163.86 129.11 [1x4 double] 41 163.17 139.32 [1x4 double] 68 164.26 147.47 [1x4 double] 61 164 164 [1x4 double] 63 163.86 175.11 [1x4 double] 33 164 195 [1x4 double] 54 163.44 203 [1x4 double] 33 164 211 [1x4 double] 80 167.95 37.962 [1x4 double] 9 168 233 [1x4 double] 9 168 238 [1x4 double] 9 168 243 [1x4 double] 63 171.89 16.857 [1x4 double] 68 180.53 39.265 [1x4 double] 48 184.88 67.604 [1x4 double] 106 183.6 118.53 [1x4 double] 6 178.5 127 [1x4 double] 85 185.26 134.76 [1x4 double] 85 185.26 164.24 [1x4 double] 85 185.26 179.76 [1x4 double] 6 178.5 212 [1x4 double] 85 185.26 238.76 [1x4 double] 85 184.24 16.259 [1x4 double] 72 185.68 76.653 [1x4 double] 63 185.86 88.111 [1x4 double] 41 185.17 98.317 [1x4 double] 63 185.86 107.11 [1x4 double] 33 186 127 [1x4 double] 74 187.2 152.31 [1x4 double] 63 185.86 192.11 [1x4 double] 72 185.68 203.65 [1x4 double] 33 186 212 [1x4 double] 68 186.26 219.47 [1x4 double] 41 185.17 230.32 [1x4 double] 48 188.4 37.875 [1x4 double] 9 199 43 [1x4 double] 48 206.88 161.6 [1x4 double] 6 200.5 235 [1x4 double] 9 204 43 [1x4 double] 61 208 125 [1x4 double] 105 207.79 139.21 [1x4 double] 63 207.86 153.11 [1x4 double] 61 208 170 [1x4 double] 74 209.2 181.31 [1x4 double] 61 208 192 [1x4 double] 41 207.17 207.32 [1x4 double] 63 207.86 216.11 [1x4 double] 54 207.44 227 [1x4 double] 33 208 235 [1x4 double] 41 207.17 242.32 [1x4 double] 9 209 43 [1x4 double] 9 212 117 [1x4 double]
Immediately, you can see that the data is much more easily visible with a table. You get scrolling and editing if you open t in the Variable Editor.
You can plot directly with the data in a table variable:
imshow(bw) hold on plot(t.Centroid(:,1),t.Centroid(:,2),'r*') axis([.5 100.5 .5 100.5]) hold off
Logical indexing comes in handy for filtering objects based on various criteria. For example, the following code removes objects with area less than 10 from the table.
t(t.Area < 10,:) = [];
And you can easily export tables to various formats. Here is our table exported to a .CSV file and then read into Excel.
writetable(t,'regionprops.csv');
I encourage everyone who has used regionprops to play around with tables (assuming you have MATLAB R2013b) to see how they might be useful in your own work.
Get
the MATLAB code
Published with MATLAB® R2013b
I am surprised when many of the singular values of a nonsymmetric variant of the Hilbert matrix turn out to be nearly equal to $\pi$. The explanation involves the Fourier series for a square wave.
My very first project in matrix computation involved the Hilbert matrix. It had been introduced in 1894 by the famous mathematician David Hilbert as what we would today call the coefficient matrix in the normal equations for least squares approximations by polynomials with their monomial basis. The elements of an $n$ -by- $n$ Hilbert matrix are
$$ a_{i,j} = \frac{1}{i+j-1}, \ \ i,j = 1,...,n $$
Some years ago I wanted a nonsymmetric test matrix, so I changed the plus sign in the denominator to a minus sign, but those elements blew up along the subdiagonal, so I also changed the -1 to +1/2.
$$ a_{i,j} = \frac{1}{i-j+1/2}, \ \ i,j = 1,...,n $$
I refer to this as the nonsymmetric Hilbert Matrix. Here is the 5-by-5 with format rat.
format rat
n = 5;
[I,J] = ndgrid(1:n);
A = 1./(I-J+1/2)
A = 2 -2 -2/3 -2/5 -2/7 2/3 2 -2 -2/3 -2/5 2/5 2/3 2 -2 -2/3 2/7 2/5 2/3 2 -2 2/9 2/7 2/5 2/3 2
I was amazed when I happened to compute the singular values. Here are all of the singular values of the 30-by-30 nonsymmetric Hilbert matrix.
format long
n = 30;
[I,J] = ndgrid(1:n);
A = 1./(I-J+1/2);
sigma = svd(A)
sigma = 3.141592653589795 3.141592653589794 3.141592653589794 3.141592653589794 3.141592653589794 3.141592653589794 3.141592653589793 3.141592653589793 3.141592653589793 3.141592653589793 3.141592653589793 3.141592653589792 3.141592653589792 3.141592653589792 3.141592653589791 3.141592653589779 3.141592653589737 3.141592653588278 3.141592653559363 3.141592653039148 3.141592644613932 3.141592521935767 3.141590919627153 3.141572215205717 3.141378147613475 3.139604197477250 3.125546289482884 3.032461817019588 2.564498174433655 1.105241141304968
Incredibly, half of them are equal to pi to the full length of our long format, that's 16 digits. And all but the last few are equal to pi to at least a couple of digits. How did that happen?
Let's do the 30-by-30 SVD again, this time using the Symbolic Toolbox with 100-digit arithmetic. The computation takes about a minute and a half on my machine.
tic digits(100) n = 30; [I,J] = ndgrid(sym(1:n)); A = 1./(I-J+1/2); sigma = svd(A); toc
Elapsed time is 83.939032 seconds.
Here is the first singular value, to 100 digits.
sigma1 = sigma(1)
sigma1 = 3.141592653589793238462643383279502884197169352107009745551549597527462440105298598090594765999823398
Here is a logarithmic plot of the difference between the singular values and $\pi$.
semilogy(abs(pi - sigma),'.') xlabel('n') ylabel('|\pi - \sigma|');
The first singular value equals $\pi$ to 44 digits and, as our double precision computation predicted, fifteen of the thirty singular values equal $\pi$ to at least 15 digits.
Seymour Parter was a visiting professor at Stanford when I was a grad student there in the early 1960s. I took an individual reading course under his supervision. He went from Stanford to the University of Wisconsin, where he spent the rest of his career as a Professor of Mathematics and Computer Sciences. He was President of SIAM in 1981 and '82.
In 1985, the Second SIAM Conference on Linear Algebra was held at North Carolina State University in Raleigh. I demonstrated MATLAB during the coffee breaks. We didn't really have a booth, just a desk or table. I used this example in my demonstration and asked conference attendees if they could explain the appearance of $\pi$. Seymour said that he could. He had written two papers twenty years earlier about eigenvalues of Toeplitz matrices that were relevant. He soon produced another paper, referenced below, that explains the behavior of the singular values of my nonsymmetric Hilbert matrix. It is related to a result from 1920 by the famous Hungarian mathematician Gabor Szego.
Our matrix $A$ is constant along its diagonals, so it is a Toeplitz matrix. So is the $2n$ -by- $2n$ matrix $B$ constructed from $A$ by the MATLAB array operation
Z = zeros(size(A)); B = -i*[Z A; -A' 0];
The matrix $B$ is both Hermitian and Toeplitz. Its eigenvalues are real and occur in plus/minus pairs. The absolute values of the eigenvalues of $B$ are the singular values of $A$, each repeated twice.
Szego's 1920 theorem says that the eigenvalues of the Toeplitz matrix formed from the coefficients in the Fourier expansion of a function converge to the max and min of that function.
Parter realized two important facts. First, that he could analyze the singular values of my nonsymmetric matrix $A$ by looking at the eigenvalues of the Hermitian matrix $B$. And second, that the elements of $A$ and $B$ are the reciprocals of odd integers, which are the coefficients of the Fourier expansion of a square wave. So Szego's 1920 theorem implies that the largest eigenvalues of $B$ and hence the largest singular values of $A$ converge to the maximum of that square wave.
There is already a hint of what is going on in an old MATLAB demo named xfourier. The first line of xfourier says
This example shows (graphically) how the Fourier series expansion for a square wave is made up of a sum of odd harmonics.
The core of xfourier is a code fragment similar to this.
t = -pi:pi/256:pi; y = zeros(size(t)); for k = 1:2:99 y = y + sin(k*t)/k; end plot(t,y)
The plot shows two things. First you notice the overshoot at the ends of the intervals known as the Gibbs phenomenon. But more important for our purposes here is the convergence in the interior of the intervals to values near plus and minus 0.8.
Let's use more terms and take a closer look.
t = -pi:pi/256:pi; y = zeros(size(t)); for k = 1:2:599 y = y + sin(k*t)/k; end plot(t,y) axis([0 pi .76 .81])
We find the sum is converging pointwise to a value near 0.7854. The fragment code is approximating the Fourier sine expansion
$$ s(t) = \sum_{odd\ k}\frac{\sin{kt}}{k} $$
The coefficients are the reciprocals of odd integers, which are the same as the elements of the nonsymmetric Hilbert matrix.
The Fourier series converges pointwise to $\pm\pi/4$. If we trace this fact back through Parter's analysis and Szego's theorem, we have an explanation for the fact that the largest singular values of the nonsymmetric Hilbert matrix are converging to $\pi$.
Evgenij E. Tyrtyshnikov is a Professor and Deputy Director of the Institute of Numerical Mathematics of the Russian Academy of Sciences in Moscow. In 1991 he studied the behavior of the small singular values of this matrix. And in 1992 he considered the matrices where the 1/2 is replaced by any noninteger value.
Seymour V. Parter, On the distribution of the singular values of Toeplitz matrices, Linear Algebra and its Applications 80, 1986, 115-130, <http://www.sciencedirect.com/science/article/pii/0024379586902806>
Evgenij E. Tyrtyshnikov, Cauchy-Toeplitz matrices and some applications, Linear Algebra and its Applications 149, 1991, 1-18, <http://www.sciencedirect.com/science/article/pii/002437959190321M>
Evgenij E. Tyrtyshnikov, Singular values of cauchy-toeplitz matrices, Linear Algebra and its Applications 161, 1992, 99-116, <http://www.sciencedirect.com/science/article/pii/002437959290007W>
Get
the MATLAB code
Published with MATLAB® R2014a
Alfonso is the guy you want to have in your community. Brilliant, generous, and energetic, he moves through Cody like an artist, insightful answers flying from his fingertips. Through his Cody interactions alone, I suspect he has taught us more MATLAB than many a professor.
He was kind enough to let us do a virtual interview with him. Spoiler alert! In one of his answers he says he solves easy Cody problems while he’s on the phone. I like that. I’m lucky if I can draw a stick man while I’m talking on the phone.
A: I am a computational neuroscience researcher, and my main focus is on neuroimaging methods. I use Matlab constantly in my research for all sorts of projects. It may be simply using well-established neuroimaging tools (e.g. SPM), developing my own (e.g. Conn), building neural models in Simulink, performing statistical analyses, or even developing real-time speech synthesizers, Matlab is always my first choice. I learned Matlab during my undergraduate studies as a telecommunication engineer, and was quickly won over by the incredible pace at which I could develop a project in Matlab compared to C or other less specialized languages. By the time I was in my graduate studies I was using Matlab almost daily, and I do not believe I have stopped since.
A: I was just browsing Matlab Central when I first saw a banner for Cody, and I was quickly hooked by the simple mini-problems, doing “just one more” over and over again. But it was really the wonderful community of users that started growing around Cody what kept me coming back, learning from and being often surprised by other people solutions, and enjoying the continuous stream of imaginative problems that have since then been posted.
A: Both Cody and the Matlab Contest were fantastic ways of learning from the solutions and the thought-process of others. The old/classic Matlab contest had this fantastic pacing, and a great sense of community (and the rush of last minute solutions!). Cody’s permanence, on the other hand, may miss the rush of a more bounded competition, but it brings in return a never-ending source of enjoyment and learning opportunities. I often find myself revisiting old problems or solutions just to see what other people may have come up with in the meantime, or following specific players just because I enjoy their way of solving problems or because I enjoy the sort of problems that they create. Also the possibility to create your own problems adds a whole different perspective to Cody, and I may have found myself spending way too much time learning the inner workings of Cody just to see “what could be done” with it.
A: Not terribly much. It is true that at times I may have gone and solved some perhaps not-terribly-interesting problems just for the thrill of catching that first position, but mostly I come back to Cody looking forward to being surprised by some new problems or solutions, or just staying in touch with the great community here.
A: My time on Cody fluctuates a lot, depending on my work or life schedule, but if I am caught up in a series of interesting problems I do have to make a real effort not to let Cody take over too much time from either.
A: It depends a lot on the problem. There are some interesting problems where I would send one solution after another until I cannot find other ways to “improve” it. Also, when other players are actively sending solutions to the same problems I may get caught in this fun collaboration where we would be trying to build on each other’s ideas. It also depends on the setting, if I am on my phone I may go for simpler or shorter problems, or if I am just relaxing having a coffee perhaps I may take the time to think about some really complex or interesting problem that I read on Cody some time back and never had the time to solve. Last, sometimes I may just feel like revisiting old problems or solutions to see if I can think of a better solution or just wondering what others may have come up with since I last visited them.
A: For me typically a good Cody problem would have a strong testsuite (avoiding simple look-up table solutions), would not have an immediately obvious “best” way of solving it, and it would not over-engineer the range of expected solutions (e.g. allow some tolerance in error checking, do not forbid or discourage the use of alternative solutions, etc.) Of course I have seen many fantastic problems that break about every of these rules, so new problem creators should probably just follow their instinct and know that their problems will get better with practice. A good place to start is probably to try to find inspiration in the way other people create their problems, I would recommend Richard Zapor, Ned Gulley, Bryant Tran (bmtran), or Tim Vaughan (the cyclist) as some of my favorite problem creators out there. Also I would recommend playing with the format to see what works best. Lately I am increasingly enjoying creating problems that attempt to somewhat explore the limits of what can be done in Cody (e.g. freepasses, Blockland) and, while this might not bring as much fun to others, I am certainly enjoying it! So again, just follow your instincts and experiment a bit.
A: The answers I enjoy the most typically would surprise me by using an unexpected approach, some Matlab functionality I never heard of or considered in the context of this problem, or just by virtue of their goofiness or ingenuity.
A: This would be a very long list. I would say that some of the players who manage to very consistently create solutions that I am sure to enjoy or learn something from are Bryant Tran, with whom I shared the thrill of “discovering” Cody, learning to write highly optimized solutions, and unearthing quite a few hacking tricks, and Tim, who not only creates beautifully efficient code but also he is the one player I can always count on to be able to solve pretty much any problem I throw at him.
Among the problem creators that I enjoy the most I would “special mention” the prolific Richard Zapor (with a ton of fantastic problems in his GJam competition series; his PACMAT series was also ridiculously fun), and Ned Gulley (his problems consistently hit the perfect balance between simplicity of definition, and complexity of possible solutions/approaches; e.g. Remove the air bubbles, Return its own character count, Stride of the longest skip sequence, are really beautiful problems that I often come back to). I also had great fun with Khaled Hamed’s functional programming problems (problems 1198, 1210, 1224), as well as the many hacking problems out there –– starting with Óscar’s “Now!” series (problems 901, 907), Frequency Domain’s “Get me!” series (problems 1745, 1746, 1747, 1752) and my own “Singularity 2.0” series (problems 1769, 1770, 1764).
A: Those sorts of interactions make some of the most memorable moments for me in Cody (and they are always more like conversations or shared experiences, only perhaps disguised as battles). In fact they are probably what have made many of the problems and solutions that I mentioned above particularly meaningful or enjoyable to me. As an example when the “Get me!” series of hacking problems was created it brought together many players trying to find better ways to build strong testsuites that would protect against hacking attempts (and/or find better hacks to break stronger testsuites), and we built from each other’s ideas in what ended up being a remarkable collective experience of discovery. And even though it was somewhat sad to me in the end that we found some sort of all-purpose hack that could not be currently protected against, which unavoidably ended our quest, the entire process was a genuinely fun/bonding/instructive experience.
A: There are a few of these “tricks”. Using command syntax (e.g. print file) instead of function syntax (e.g. print(‘file’)), using str2num to define a “constant” vector or matrix, and the somewhat ugly but necessary ans trick, would probably fall into the “recommendable” category, and they will help you shave a few points of almost any solution. The infamous regexp trick (as well as all sorts of hacking tricks) would probably fall into the “not recommendable” category, as these inevitably flatten the score function (they give any solution, independently of how clever or absurd it is, the exact same score). Despite these tricks, my favorite pleasures really come from finding new simpler ways to approach a problem, or applying some linear algebra tricks, or these sort of discoveries which lead to some small improvements in score in a way that was somewhat unexpected or surprising.
A: This would again be a very long list. Most importantly I have widened considerably the range of approaches that I consider now for almost any “real world” problem, but in practical terms I have also learned a lot about regexp functionality, Java usage in Matlab, tricks to speed-up your code, and a lot of new Matlab functionality that I was not aware of (e.g. uint64 operations, timeseries objects, map containers, etc.), just to mention a few examples.
A: I believe mostly I would like to better be able to search/browse through old problems and solutions. For example, being able to search problems/solutions that I (or specifically somebody else) ‘liked’ would make the task of coming back to revisit some of the most memorable solutions/problems or discovering new solutions that you may have missed much simpler (in this regard something like a “self-like” marking would be useful, as I would love to hear/follow which of their own problems/solutions some players are most “proud” of).
A: I love the Blogs section (Cleve’s and Loren’s are some of my favorites), very sporadically contribute to the Newsgroup or Answers sections, and have been lately doing my first steps on Trendy (and of course was a big fan of the Matlab Contest).
A: accumarray, bsxfun, and cellfun/arrayfun/structfun are incredibly flexible and useful functions that are somewhat underused in my opinion.
A: svd (I would spend my day doing linear algebra tricks)
A: To close I would just like to reiterate how much I appreciate the excellent work of the Mathworks team which makes all of these shared learning and fun experiences possible, as well as to thank the amazing Cody community for sharing so much of their time, brainpower, and humor with the rest of us.
]]>I have written blogs before on coding best practices in MATLAB and they have always generated a lot of great feedback. I recently came across a language-neutral list of best practices for scientific computing, written by Greg Wilson of Software Carpentry.
If you don't want to start by reading the whole article, I strongly encourage you to read Greg's summary. He breaks the practices into eight categories. I have to say, I think they are very well stated, starting with "Write programs for people, not computers.".
I have to say, I think Greg hit the nail on the head.
Do You Violate these Best Practices?
The most frequent violations I see include not using a debugger and, a lot of code is trying to optimize it before being sure it is truly correct. Which practice is your Achilles' heel? Let me know here.
Get
the MATLAB code
Published with MATLAB® R2013b
At about the same time, File Exchange contributor Wolfgang Schwanghart was developing a similar approach, which he made available as the TopoToolbox.
I just heard from Wolfgang that a new version of the TopoToolbox is available. It is described in the paper Schwangart and Scherler, TopoToolbox 2 – MATLAB-based software for topographic analysis and modeling in Earth surface sciences, Earth Surface Dynamics, 2014.
The new version features improved memory efficiency and speed, new analysis algorithms, and interactive apps for visual exploration and interaction with DEMs. Here's a figure from the paper:
You can find TopoToolbox here. If you are interested in analytical and interactive tools for DEM analysis, you should give it a try.
]]>I'd like to welcome back guest blogger Spandan Tiwari for today's post. Spandan, a developer on the Image Processing Toolbox team, posted here previously about circle finding and homomorphic filtering.
A few weeks ago Steve wrote about watershed transform and how to use it for image segmentation. I thought I would continue the topic and highlight another neat application for watershed transform - Maze Solving. Steve discussed this problem in his earlier post on Exploring Shortest Paths and showed a solution based on Geodesic distance transform (bwdistgeodesic). In this post I show another way for solving a maze using watershed transform.
Here's a maze. We would like to find the solution , i.e. a path between the two points (entry and exit points) shown in red.
I = imread('http://blogs.mathworks.com/images/steve/2011/maze-51x51.png'); I = logical(I); imshow(I) hold on; plot([6 526], [497 517], 'ro', 'MarkerSize', 15,'LineWidth', 2); hold off;
I should mention that in my approach I am considering only 'standard' or 'perfect' mazes that have one and only one path from the entry point to the exit point, and some other nice properties such as absence of loops. I will call these mazes well-behaved mazes.
Here is the key idea behind using watershed transform for solving a maze. The landscape (or surface) of the image of a well-behaved maze has two catchment basins. The watershed between those catchment basins is the solution path for the maze! (For more details on the definition of watershed and catchment basin see Steve's newsletter article The Watershed Transform: Strategies for Image Segmentation). Let's look at the catchment basins as produced by watershed.
L = watershed(I); imshow(L,[])
The two basins are shown in grey and white. If you follow the interface (boundary) between the two basins (grey and white regions), you can see notice that it is the solution (path) to the maze.
The good news is that at this point we have already solved the bulk of the problem. All we need to do now is to extract the interface between the two basins. There are different ways to do this. Here's the one that I like, mainly because I find it to be robust. Another minor reason I like it is that the final path stays parallel to the maze lines which gives it a nice clean look.
First we create a new image retaining only the original image region from one of the catchment basins.
L1 = L == 2; I1 = L1.*I; imshow(I1)
Next we use watershed again to get catchment basins for this new image.
L2 = watershed(I1); imshow(L2,[])
Notice that the catchment basin in this new image has shrunk roughly by half the width of the maze paths as compared to the same catchment basin in the original image. This can be seen clearly using imshowpair. It highlights the area where the images are different using color (green, in this case)
imshowpair(L,L2)
As the last step we just extract the path shown as the green region, which is the region where the two catchment basins differ.
img1 = L == 2; img2 = L2 == 2; path = img1 - img2;
We can visualize the path on the maze nicely using Steve's imoverlay function.
P = imoverlay(I, path, [1 0 0]); imshow(P)
Here are few other results on different mazes.
Get
the MATLAB code
Published with MATLAB® R2013b
Magic Stars are a little like Magic Squares. My grandson was introduced to them in his junior high school math class.
This post is a family enterprise. My grandson Ben heard about five-point magic stars in his junior high school math class, where they tackled them by hand. Ben told his mom Kam, my daughter, who wrote a MATLAB program, and then told me.
Here is the function that I use to view the stars.
type observatory
function observatory(v) % observatory View Magic Stars. % observatory(v), where v is a vector of 10 or 12 integers. clf n = length(v)/2; d = pi/n; t = 0:2*n-1; r0 = .42 + .20*(n == 5); r = 1 - r0*(mod(t,2)==1); x = r.*sin(t*d); y = r.*cos(t*d); p = plot(x,y,'o','markersize',28,'markeredgecolor',[0 0 .66]); axis(1.25*[-1 1 -1 1]) axis square set(gca,'xtick',[],'ytick',[]) if n == 5 h = line([x(1:4:2*n) x(3:4:2*n) x(1)], ... [y(1:4:2*n) y(3:4:2*n) y(1)]); elseif n == 6 h = [line([x(1:4:2*n) x(1)],[y(1:4:2*n) y(1)]) line([x(3:4:2*n) x(3)],[y(3:4:2*n) y(3)])]; else error('Can only observe 10 or 12 points') end set(h,'color',[0 .66 .66]) for k = 1:2*n text(x(k),y(k),int2str(v(k)),'fontsize',18,'horiz','center') end end
There are lots of kinds of magic stars. Here are the most basic. Draw a regular n-point star the way you learned to draw a 5-point star. Connect the n points evenly spaced around a circle by n lines between alternate points. This produces a regular n-gon in the interior. The vertices of the outer points, together with the vertices of the interior polygon, gives you 2n points. Here are the five and six point stars, with the points numbered in consecutive clockwise order.
observatory(1:10) snapnow observatory(1:12)
To turn one of these stars into a magic star, you need to rearrange the numbers so that the sums along each of the n lines are equal, just like the sums along the rows and columns in a magic square. It turns out that this is impossible for the 5-point star and we will have to modify the rules. But it is possible for the 6-point star.
The sum of the integers from 1 to 10 is S = 55. There are five lines in the 5-point star and each integer is involved in the intersection of two of these lines, so if the sums along each line are to be equal, that sum must be L = (2/5)S = 22. We are seeking a permutation p of the first ten integers that satisfies the five equations.
$$ p_1+p_2+p_4+p_5 = L $$
$$ p_3+p_4+p_6+p_7 = L $$
$$ p_5+p_6+p_8+p_9 = L $$
$$ p_7+p_8+p_{10}+p_1 = L $$
$$ p_9+p_{10}+p_2+p_3 = L $$
The subscripts in each of these equations can be read off by following one of the lines across the star.
An exhaustive search of all 10! permutations of the integers 1 through 10 reveals that there are no solutions to these equations, so no perfect five-point magic star exists.
But do not despair. Relax the rules by allowing more players. Ben's class was given the challenge of finding a magic 5-point star using nine of the ten integers between 1 and 10, and also using 12. They had to be careful. The sum S needs to be a multiple of 5 because the line sum is always going to be L = (2/5)S. So they learned that they had to exclude either 2 or 7. Let's first try leaving out the 7.
I'm going to give the integer 12 a special role by fixing it at the apex of the star. This eliminates all the stars that are obtained from the ones that I will find by simple rotations, and it also eliminates all the stars where 12 is at an internal point. It leaves only nine integers left to permute. Here are the parameters.
clear
n = 5;
v = [1:6 8:10]; % Omit 7.
S = sum(v) + 12
L = 2/n*S
S = 60 L = 24
Use the MATLAB function perms to generate all possible permutations of the nine free integers.
P = perms(v);
Insert a column of 12's in front to produce a matrix with 9! (that's 9 factorial) rows and 10 columns.
P = [12*ones(size(P,1),1) P]; sizeP = size(P)
sizeP = 362880 10
See how much storage we're using.
whos
Name Size Bytes Class Attributes L 1x1 8 double P 362880x10 29030400 double S 1x1 8 double n 1x1 8 double sizeP 1x2 16 double v 1x9 72 double
My laptop can handle 29 megabytes.
We are ready to search for the solutions.
f = find(P(:,1)+P(:,2)+P(:,4)+P(:,5) == L & ... P(:,3)+P(:,4)+P(:,6)+P(:,7) == L & ... P(:,5)+P(:,6)+P(:,8)+P(:,9) == L & ... P(:,7)+P(:,8)+P(:,10)+P(:,1) == L & ... P(:,9)+P(:,10)+P(:,2)+P(:,3) == L)
f = 85418 99312 133775 139937 205213 222969 257655 272091 284491 298403 322934 361090
Out of the over one-third of a million permutations, these are the indices that satisfy the five line equations. Here are those solutions.
P = P(f,:)
P = 12 8 9 1 3 10 4 6 5 2 12 8 5 3 1 10 6 4 9 2 12 6 10 4 2 9 1 8 5 3 12 6 5 2 4 9 8 1 10 3 12 4 9 2 6 5 8 3 10 1 12 4 10 6 2 5 3 8 9 1 12 3 5 8 1 9 2 4 10 6 12 3 10 1 8 9 4 2 5 6 12 2 9 4 6 10 1 3 5 8 12 2 5 6 4 10 3 1 9 8 12 1 9 8 3 5 2 6 10 4 12 1 10 3 8 5 6 2 9 4
Twelve solutions were found. The last solution has the integers from 1 to 5 in the interior.
p = P(12,:) observatory(p)
p = 12 1 10 3 8 5 6 2 9 4
Our search for stars formed from the integers from 1 to 12, leaving out 7 and 11, found 12 solutions. But the last six are actually reflections of the first six. Look at the solution array P, and ignore the first column with the 12s. If you read the last six rows backwards, from right to left, you will get the first six rows in a different order. You can also think of this as clockwise versus counter clockwise indexing.
So, we have found six distinctly different solutions with 12 at the apex. There are three kinds of transformations that can be made on these solutions. This will produce all of the possible magic stars that can be made from this set of integers.
Reflect, about the vertical axis.
p = p([1 10:-1:2]) observatory(p) snapnow
p = 12 4 9 2 6 5 8 3 10 1
Rotate, counter clockwise.
p = p([3:10 1:2]) observatory(p) snapnow
p = 9 2 6 5 8 3 10 1 12 4
Invert, interchange interior with exterior.
p = p([6:10 1:5]) observatory(p) snapnow
p = 3 10 1 12 4 9 2 6 5 8
Applying two reflections, five rotations, and two inversions to each of the original six solutions produces a grand total of 6*2*5*2 = 120 different magic stars.
You can modify my code to leave out the 2 instead of the 7. Use v = [1 3:10] instead of v = [1:6 8:10]. You will find that the search for solutions comes up empty. Magic stars formed from the integers 1 through 12 without the 2 and 11 do not exist.
You can also try including 11, but leaving out 2 and 6. Use v = [1 3:5 7:11] instead of v = [1:6 8:10]. Is there an isomorphism between these stars and we ones I described above?
We can make a 6-point star using the integers from 1 to 12, so we don't have to bend the rules. But I am going to restrict the search space and be more careful about memory usage.
Here is the line sum.
clear n = 6; S = sum(1:2*n); L = 2/n*S
L = 26
I have decided to restrict the search space to stars with 11 at the apex and 12 at the base. So I need to find all permutations of the remaining integers from 1 to 10. This would use ten times as much memory as I needed to do the 5-point star and my laptop can't handle it. But instead of storing the integers as the MATLAB default eight-byte double precision numbers, I can use one-byte uint8s and reduce the memory requirement by a factor of eight.
Find all permutations of ten integers, stored in single bytes.
v = uint8(1:10); P = perms(v);
Insert a column of 11s and a column of 12s in the appropriate places.
c = ones(size(P,1),1,'uint8'); P = [11*c P(:,1:5) 12*c P(:,6:10)]; clear c
The matrix has 10! rows, which is about 3.6 million.
sizeP = size(P)
sizeP = 3628800 12
But the storage required is not much more than above.
whos
Name Size Bytes Class Attributes L 1x1 8 double P 3628800x12 43545600 uint8 S 1x1 8 double n 1x1 8 double sizeP 1x2 16 double v 1x10 10 uint8
There are six line equations to be satisfied.
f = find(P(:,1)+P(:,2)+P(:,4)+P(:,5) == L & ... P(:,3)+P(:,4)+P(:,6)+P(:,7) == L & ... P(:,5)+P(:,6)+P(:,8)+P(:,9) == L & ... P(:,7)+P(:,8)+P(:,10)+P(:,11) == L & ... P(:,9)+P(:,10)+P(:,12)+P(:,1) == L & ... P(:,11)+P(:,12)+P(:,2)+P(:,3) == L)
f = 235278 520708 641908 1321733 1522350 1930425 2012865 2613900
Here are the eight solutions found in the 3.6 million rows. Again, the last four are reversals of the first four.
P = P(f,:)
P = 11 10 4 2 3 8 12 6 9 1 7 5 11 9 6 1 5 7 12 4 10 2 8 3 11 9 3 1 5 10 12 4 7 2 8 6 11 7 4 2 6 8 12 3 9 1 10 5 11 6 8 2 7 4 12 10 5 1 3 9 11 5 7 1 9 6 12 8 3 2 4 10 11 5 10 1 9 3 12 8 6 2 4 7 11 3 8 2 10 4 12 7 5 1 6 9
I think the first one is quite attractive. I have not investigated all the possible symmetries and transformations.
observatory(P(1,:))
Marian Trenkler, Magic Stars, PME Journal 11, 2004, 549-554, <http://math.ku.sk/~trenkler/MagicStars.doc>
Harvey Heinz, Magic Stars Index Page, <http://www.magic-squares.net/magic_stars_index.htm>
Get
the MATLAB code
Published with MATLAB® R2014a
Today, David Garrison and Sarah Wait Zaranek, today's guest bloggers, would like to ask for your help in tracking down interesting examples using MATLAB graphics and UI building capabilties.
We are always looking for interesting MATLAB code. Your examples help us understand all the things people do in MATLAB and they help us plan for future enhancements. With that in mind, we'd like to ask you to think about interesting graphics or UIs you've produced with MATLAB. We are looking for highly customized UIs or complicated graphics that show off the power of MATLAB.
Any code that you send us will be used for internal purposes only.
However, we may contact you to see if you are interested in letting us use your example as inspiration for an addition to the MATLAB Plot Gallery or as a part of a follow-up guest blog post. If you don't want to be contacted, please let us know.
To help inspire you, here are some complex visualizations from the MATLAB File Exchange
and here are some interesting applications examples.
And...if you are just getting started using graphics or making applications in MATLAB, here are some of our favorite resources to help get you up and running in no time:
So please look for graphics or GUI examples that you are really proud of and that show off the power of MATLAB. Please email them to us at david.garrison@mathworks.com or post a link in the comments here. If possible, try to include all the code and data required to run your application or to create your special plot or graphic.
Thanks to Loren for letting us make our request on her blog.
Get
the MATLAB code
Published with MATLAB® R2013b
I'd like to welcome back my fellow MATLAB Central blogger Brett Shoelson for the last in a three-part series on extracting curve values from a plot. You can find Brett over at the File Exchange Pick of the Week blog, or you can check out his many File Exchange contributions. -Steve
If you've followed Steve's blog for the past couple weeks, you'll know that Steve has graciously allowed me to demonstrate how one might extract real data from a graphical depiction of the data. In this final post in this three-part series, I will use the coordinates I extracted from the curve of interest to fit a predictive model to those data. Recall that I showed two approaches to determine the x- y- coordinates of the efficiency curve:
First, load the variables xs, ys, and bb, which were computed in the previous post.
tempfile = [tempname '.mat']; url = 'http://blogs.mathworks.com/images/steve/2013/curve.mat'; urlwrite(url,tempfile); s = load(tempfile); xs = s.xs; ys = s.ys; bb = s.bb; delete(tempfile);
Clearly, the results obtained using regionprops were better than those I got using bwtraceboundary. We may as well use those values for our curve fit. Before we do that, though, it will be useful to transform the data to account for the fact that the extracted coordinates are in units of pixels, rather than units of flow rate and efficiency. I do that manually here:
flowLims = [0 240]; efficiencyLims = [0 100]; % Scale data to specified limits for curvefit % Scale xs: xs = (xs-bb(1))/bb(3)*diff(flowLims)+flowLims(1); % Scale ys: ys = (bb(2)+bb(4)-ys)/bb(4); %Percentages = Efficiencies
(By the way, if you ever need to fit data to anything other than a very simple polynomial, and if you don't have the Curve Fitting Toolbox, I hope this motivates you to get it; it is is an excellent tool that will make your life much easier!)
You can fit to just about any type of equation that makes sense. (You can even provide a custom equation, if you have one in mind.) I can, for instance, get a reasonably good robust "Bisquare" fit using a two-term exponential. Or, since I just want to represent my data, and not force them to a specific form, I could fit them to a smoothing spline.
cftool(xs,ys)
When you've interactively selected your fit options, you can readily tell MATLAB to generate code, with which you can apply the same routine to subsequent images (data sets).
Here's the code that cftool told me to use:
[xData, yData] = prepareCurveData( xs, ys ); % Set up fittype and options. ft = fittype( 'smoothingspline' ); opts = fitoptions( 'Method', 'SmoothingSpline' ); opts.SmoothingParam = 0.004; % Fit model to data. [fitresult, gof] = fit( xData, yData, ft, opts ); % Plot fit with data. figure( 'Name', 'untitled fit 1' ); h = plot( fitresult, xData, yData ); legend( h, 'ys vs. xs', 'untitled fit 1', 'Location', 'NorthEast' ); % Label axes xlabel( 'xs' ); ylabel( 'ys' ); grid on
We could, for illustration, calculate the efficiency of the pump when the flow rate is 100 m^3/hr:
flow100 = fitresult(100)
flow100 = 0.5997
If you refer to the original image, you'll see that that value perfectly reflects the efficiency at the specified flow rate!
I don't really have a good sense for how well this exact approach will work for subsequent images. I tried to make the code fairly general, but the automatability of the problem depends markedly on the ability to exploit similarities in the data set (images). Therein lies the art of image processing!
Get
the MATLAB code
Published with MATLAB® R2013b
Sometimes I need to construct a matrix with certain symmetries. There are a bunch of tools in MATLAB that are well suited for such tasks. In today's post, I will mention the ones I use most often.
I want to start with perhaps the most ignored funtionality for ensuring symmetry of a certain kind is the transpose operators .' (ctranspose) and ' (transpose). The reason these can be so so helpful is that during the construction of a matrix A, that we know should be symmetric (such as a covariance matrix) or hermetian (such as a hessian), the corresponding elements A(i,j) and A(j,i) may be computed with different expressions, or with the terms being aggregated in different orders. If that happens, these elements may differ due to numerical considerations. However, if the next step in the algorithm requires A to be symmetric, then the code may fail. So, to ensure symmetry, assuming A is real, you can do something like this:
A = (A + A')/2
Of course, another thing you might do is only calculate the elements on the diagonal and either above or below, and then fill in the remaining elements. A third, strategy, if you have control of the algorithm following the matrix construction, is to only use the upper or lower triangular elements, reducing the need to symmetrize the input.
Also, if for some reason I was working in higher dimensions and needed symmetry in a "plane", I could use permute and its companion ipermute.
There are other kinds of symmetries that arise in numerical computations. One example is for windowing signals for certain signal processing applications. In this case, the window that is applied is generally symmetric, but in a different way. For a window W of length n, the symmetry can be expressed like this:
W(end:-1:(ceil(n/2)+1)) = W(1:floor(n/2))
NOTE: I was careful in the above expression to not copy the middle element of an odd-length array :-) !
If I had just computed the first half of the window (being careful about odd and even lengths), I could just use one of these functions to get the window symmetry instead of working out the expression above, using the : (colon) operator and some math.
Here's a list (likely not exhaustive) of other tools to help produce symmetric arrays.
And here's one I basically never use, because I find that I have to do a lot of special case coding if I do.
When do you need and use symmetry? In what application? And how do you achieve that in your code. Let us know here.
Get
the MATLAB code
Published with MATLAB® R2013b
I'd like to welcome back my fellow MATLAB Central blogger Brett Shoelson for the second in a three-part series on extracting curve values from a plot. You can find Brett over at the File Exchange Pick of the Week blog, or you can check out his many File Exchange contributions. -Steve
Continuing from last week's guest post, we're heading down the path of extracting real data from the "efficiency curve" in a graphical representation of those data. To recap, here's the original image:
Here are the steps we need to repeat from part 1
url = 'http://blogs.mathworks.com/images/steve/2013/example-pump-perf-curve.jpg'; chart = imread(url); bw = im2bw(imcomplement(rgb2gray(chart))); filled = imfill(bw,'holes'); cc = bwconncomp(filled); stats = regionprops(cc,'Area'); A = [stats.Area]; [~,biggest] = max(A); filled(labelmatrix(cc)~=biggest) = 0; bb = regionprops(filled,'BoundingBox'); %We'll use this later! bb = bb.BoundingBox; chart(~repmat(filled,[1 1 3])) = 0; curveOfInterest = rgb2ycbcr(chart); curveOfInterest = curveOfInterest(:,:,2);
By the end of the first post, I had managed to reduce the SNR so that I ended up with an image (curveOfInterest) that contains a single region...that partially isolates the two blue curves in the top axes:
Many (perhaps most) image processing problems require some approach to "segmenting" the image--creating a binary mask that shows white ("true," or 1) where we have signal, and black ("false," or 0) elsewhere. Now that I've masked and converted the original image to YCbCr ("luminance/chrominance/chrominance"), and selected the second (Cb) plane of the resulting image, segmentation is relatively easy. I'm going to use SegmentTool to explore options for doing that:
Simple thresholding works well here (try it!), but I'm going to use a different approach--one that sometimes succeeds where thresholding fails.
curveOfInterest = imextendedmax(im2double(curveOfInterest),0.30,8);
imshow(curveOfInterest)
title('Binary Curves')
For subsequent analysis, and extraction of data, it will be useful to morphologically alter the image. MorphTool will be your best friend in that regard! ;)
After playing around a bit, I decided to dilate the image, and then to skeletonize it and remove spurs. (bwmorph makes both of those latter processes easy.):
curveOfInterest = imdilate(curveOfInterest, strel('Disk',6)); curveOfInterest = bwmorph(curveOfInterest, 'skeleton', Inf); curveOfInterest = bwmorph(curveOfInterest, 'spur', 11); imshow(curveOfInterest)
If the curves in that image don't look contiguous, recognize that that's just a screen-resolution issue. There's now exactly one region in this image:
cc = bwconncomp(curveOfInterest); cc.NumObjects
ans = 1
And now for the difficult part
Okay, we're well on our way. We've eliminated the black and red curves that we don't need, and we've got the curve of interest reduced to a skeleton that contains the minimum information that we do need. But we still have to find a way to extract only the efficiency curve, and to discard the curve representing "334-mm" of head. We can't use color; the curves in the original image are identical. So we're left with differentiating the curves based on their morphological (shape) characteristics.
I'm going to do this two different ways
Sometimes it's useful to be able to trace the boundary of a region, or (as we do here), to trace along a particular path. bwtraceboundary will allow us to do just that.
If you've never used it before, take a moment to read the documentation for bwtraceboundary. You'll find that it allows you to operate on a binary image by specifying a starting point (P), an initial trajectory (fstep) along a path in the image, and a direction (dir) in which you want to trace the path. ("Direction" here allows you to specify the overall orientation of the boundary you want to trace.)
Automating the determination of the starting point
Recalling that the curve of interest always starts in the lower left of the upper axes, we can readily automate the determination of the starting point. (It will be the point closest to the bottom left of the masked axes.) I'm going to use regionprops to determine the bounding box of the masked axes (currently contained in variable "filled"). However, the doc for regionprops indicates that BoundingBox is measured from the upper left; I need the lower left so I'll have to manipulate the calculation of the bounding box a bit:
targetPoint = regionprops(filled,'BoundingBox');
targetPoint = targetPoint.BoundingBox(1:2) + [0 targetPoint.BoundingBox(4)];
[ys,xs] = find(curveOfInterest);
Now it remains to get the distance from all points to the lower left corner of the bounding box, and to find the point that minimizes that distance. (I'm going to use a Statistics Toolbox function here, but this is a relatively easy function to write if you don't have that Toolbox.)
[~,ind] = min(pdist2(targetPoint,[xs,ys]));
P = [ys(ind),xs(ind)]; %(Careful here: P is represented as [row,col], not [x,y]!)
Now we need to orient the trace
Consider the results using two different orientations:
N = 1000; %If I leave this unbounded (Inf), all points on the curve will be detected! % Counterclockwise Bccw = bwtraceboundary(curveOfInterest,P,'NE',8,N,'counterclockwise'); imshow(curveOfInterest); hold on plot(P(2),P(1),'go','MarkerSize',12) plot(Bccw(:,2),Bccw(:,1),'r','LineWidth',5); title('Counterclockwise') hold off
Now clockwise.
Bcw = bwtraceboundary(curveOfInterest,P,'NE',8,N,'clockwise'); imshow(curveOfInterest); hold on plot(P(2),P(1),'go','MarkerSize',12) plot(Bcw(:,2),Bcw(:,1),'r','LineWidth',5); title('Clockwise') hold off
That's interesting; I see that if I specify a clockwise orientation for my curve (and if I use a sufficiently large N), I can constrain the boundary-tracing algorith to roughly follow the efficiency curve I'm after. However, there's a significant spur on the data that may present difficulties when we try to fit the data to a curve that reasonably represents the efficiency-versus-flow rate. If I break this up into smaller pieces and examine the output at each stage, I find that I can significantly improve this result. Here, for instance, I "step" in an easterly direction, initially specifying a counterclockwise orientation. Then I iterate over small chunks, and evaluate the results after each iteration. When the curve stops increasing in the _x-_direction (i.e., when the largest column value isn't equal to the step size), I change the orientation to keep the trace following the desired curve. (Note that there are lots of combinations of parameters that will work here; this step took some experimentation!)
sz = 25; ind = sz; nIters = 100; allBs = nan(sz*nIters,2); Pnew = P; %reset h = text(200,385,'Direction:','FontSize',14,... 'FontWeight','bold','color','w','HorizontalAlignment','left') for ii = 1:nIters if ind == sz direction = 'counterclockwise'; else direction = 'clockwise'; end set(h, 'string',['Direction: ',direction]) %B yields [ROW,COL]: B = bwtraceboundary(curveOfInterest,Pnew,'E',8,sz,direction); hold on; plot(B(:,2),B(:,1),'r.');drawnow; allBs((ii-1)*sz+1:(ii-1)*sz+sz,:) = B; [~,ind] = max(B(:,2)); % Update the starting point, Pnew: Pnew = B(ind,:); end
That result looks great, and will facilitate a good curve fit!
Anyone who has ever heard me lecture about image processing has likely heard me mention that regionprops can make difficult problems easy. Let's first use bwmorph to find the intersection of the two curves; it will present as a branch point in the skeletonized curve.
branchPoints = bwmorph(curveOfInterest,'branchpoints');
Now dilate that point and use it to break the curves at their branches.
branchPoints = imdilate(branchPoints,strel('disk',3));
curveOfInterest = curveOfInterest & ~branchPoints;
imshow(curveOfInterest)
axis([100 420 180 400])
In this zoomed view, we can see that we've now broken the curve into branches. We can also verify that we now have multiple ROIs...
cc = bwconncomp(curveOfInterest)
cc = Connectivity: 8 ImageSize: [738 1067] NumObjects: 11 PixelIdxList: {1x11 cell}
...and we can use their shape properties to select the ones we want. (I will select the ones comprising more than 10 connected pixels, and that have a positive orientation (angle from the horizontal).
stats = regionprops(cc, 'Area','Orientation','Centroid'); idx = find([stats.Area] > 10 & [stats.Orientation] > 0); curveOfInterest = ismember(labelmatrix(cc),idx); clf imshow(curveOfInterest)
figure [ys,xs] = find(curveOfInterest); plot(xs,ys,'r.','MarkerSize',7) axis ij axis equal
So that essentially wraps up the "image processing aspect" of the problem, and from here, it becomes a problem of fitting a curve. For completeness, I'll hijack Steve's blog one more time to show how the Curve-Fitting Toolbox facilitates that process. See you next week!
Get
the MATLAB code
Published with MATLAB® R2013b
The Rosser matrix is a classic matrix eigenvalue test problem.
You may not have noticed the rosser function, although it has been on your path for as long as you have had MATLAB.
which rosser
C:\Program Files\MATLAB\R2014a\toolbox\matlab\elmat\rosser.m
Here is the Rosser matrix.
R = rosser
R = 611 196 -192 407 -8 -52 -49 29 196 899 113 -192 -71 -43 -8 -44 -192 113 899 196 61 49 8 52 407 -192 196 611 8 44 59 -23 -8 -71 61 8 411 -599 208 208 -52 -43 49 44 -599 411 208 208 -49 -8 8 59 208 208 99 -911 29 -44 52 -23 208 208 -911 99
Here is most of the help entry.
rosser Classic symmetric eigenvalue test problem. This matrix was a challenge for many matrix eigenvalue algorithms. But LAPACK's DSYEV routine used in MATLAB has no trouble with it. The matrix is 8-by-8 with integer elements. It has: * A double eigenvalue. * Three nearly equal eigenvalues. * Dominant eigenvalues of opposite sign. * A zero eigenvalue. * A small, nonzero eigenvalue.
Here are the eigenvalues.
format long
e = eig(R)
e = 1.0e+03 * -1.020049018429998 0.000000000000000 0.000098048640722 1.000000000000000 1.000000000000000 1.019901951359278 1.020000000000000 1.020049018429997
As the help entry says, there is a double eigenvalue at 1000; there are three nearly equal eigenvalues around 1020; the largest positive eigenvalue is matched in magnitude by the largest negative eigenvalue; there is an exact zero eigenvalue; and there is one relatively small eigenvalue near .098.
When I was a student fifty years ago just learning about matrix computation, there were several competing methods for matrix eigenvalue computation. Francis and Wilkinson were just publishing their work on the QR algorithm, and EISPACK and MATLAB were years in the future. So the Rosser matrix represented a serious test of the methods we had available at the time. I have always admired the fact that Rosser had constructed one elegant matrix whose eigenvalue problem contained all of these challenges.
Today, the symmetric QR algorithm with Wilkinson's shift can easily handle the Rosser matrix. But we still like to have it in our collection to use in testing other algorithms.
Using the Symbolic Toolbox, we can compute the exact eigenvalues. First, find the characteristic polynomial.
p = poly2sym(charpoly(R),'x')
p = x^8 - 4040*x^7 + 5080000*x^6 + 82518000*x^5 - 5327676250000*x^4 + 4287904631000000*x^3 - 1082852512000000000*x^2 + 106131000000000000*x
The matrix is singular, so the constant term is zero. The polynomial has degree eight, but it factors nicely into linear and quadratic factors.
f = factor(p)
f = x*(x - 1020)*(x^2 - 1040500)*(x^2 - 1020*x + 100)*(x - 1000)^2
These low degree factors make it possible to obtain exact symbolic expressions for the eigenvalues.
e = solve(f)
e = 0 1000 1020 1000 10*10405^(1/2) -10*10405^(1/2) 100*26^(1/2) + 510 510 - 100*26^(1/2)
However, these symbolic expressions do not immediately reveal the numerical properties that interest us.
J. Barkley Rosser, Sr., was a colleague of my mentors George Forsythe and John Todd at the Institute for Numerical Analysis at UCLA in the 1950s. He is the most prominent figure in the front row of this photo of the INA members that I use near the beginning of my talk about "The Evolution of MATLAB". Forsythe is in the center of the group and Todd is in the back row, looking over Forsythe's shoulder.
Rosser spent much of his career as director of the famous Army Mathematics Research Center at the University of Wisconsin in Madison. He is best known for his work in number theory and logic. In 1936, he developed a version of Godel's incompleteness theorem that is now referred to as Rosser's trick. Stated informally, it says "For every proof of a theorem, there is a shorter proof of its negation.".
Get
the MATLAB code
Published with MATLAB® R2014a
These short days always make me think about the geometry of earth and sun that gives rise to seasons. And this relates to MATLAB because it gave me the opportunity to answer a question and write a little code. The question is this: for any given day, let’s say today, where are the people who will see the sun rise at exactly the same moment as me? Answering this question gave me the chance to learn a little Mapping Toolbox code.
First some basic geometry. The earth is always half dark and half illuminated by the sun. So the sunset/sunrise line (sometimes called the gray line) that separates day from night is a great circle, which is to say a circle that divides the earth evenly into two halves. This gray line is centered on the point on the earth that is precisely under the sun right now. Since we are near the winter solstice (for the northern hemisphere), the sun can be found hovering close to the Tropic of Capricorn at around 22.8 degrees south latitude. So now we have a goal: draw a great circle with its center positioned at -22.8 deg. latitude and a longitude such that the sunrise line passes through Natick.
Now let’s learn some basic Mapping Toolbox skills. This was fun for me, because two hours ago I didn’t know anything about the Mapping Toolbox. After a quick search for great circles, I started with an example from this page on circles in the documentation. From there I quickly got a map like this.
clf axesm('MapProjection','ortho') gridm on framem on setm(gca, ... 'Origin', [30 -50 20], ... 'MLineLimit', [75 -75], ... 'MLineException',[0 90 180 270]) landareas = shaperead('landareas.shp','UseGeoCoords',true); geoshow(landareas,'FaceColor',[1 1 .5],'EdgeColor',[.6 .6 .6]);
Let’s put Natick on the map. That’s me under the red dot.
natick = [42.3 -71.4]; plotm(natick(1),natick(2),'ro','MarkerFaceColor','r')
Here’s something really cool about the Mapping Toolbox: I can do a completely different map projection by changing one line in the above code. Here is the Miller Cylindrical projection. Fun!
clf axesm('MapProjection','miller') gridm on framem on setm(gca, ... 'Origin', [30 -50 20], ... 'MLineLimit', [75 -75], ... 'MLineException',[0 90 180 270]) geoshow(landareas,'FaceColor',[1 1 .5],'EdgeColor',[.6 .6 .6]); plotm(natick(1),natick(2),'ro','MarkerFaceColor','r')
I could also do a Murdoch III Minimum Error Conic Projection, but that would just be showing off.
Now down to business. Sunrise in Natick is at 7:15 Eastern Standard Time. At that moment, the sun is above the south Atlantic at 22.8 deg. south latitude and 4 deg. west longitude. We’ll use the function SCIRCLE1 to generate the great circle.
clf axesm('MapProjection','ortho') gridm on framem on setm(gca, ... 'Origin', [30 -50 20], ... 'MLineLimit', [75 -75], ... 'MLineException',[0 90 180 270]) sunLocation = [-23 -4]; grayLine = scircle1(sunLocation(1), sunLocation(2), 90); geoshow(landareas,'FaceColor',[1 1 .5],'EdgeColor',[.6 .6 .6]); plotm(sunLocation(1), sunLocation(2),'ko','MarkerFaceColor','y') plotm(grayLine(:,1), grayLine(:,2),'red') plotm(natick(1),natick(2),'ro','MarkerFaceColor','r')
Now we’re getting somewhere! Let’s switch back to the Miller projection and zoom in to the east coast of the U.S.
clf axesm('MapProjection','Miller', ... 'MapLatLimit',[-11 60], ... 'MapLonLimit',[-50 40]) gridm on framem on setm(gca, ... 'Origin', [11 -59 8], ... 'MLineLimit', [75 -75], ... 'MLineException',[0 90 180 270]) sunLocation = [-23 -4]; grayLine = scircle1(sunLocation(1), sunLocation(2), 90); geoshow(landareas,'FaceColor',[1 1 .5],'EdgeColor',[.6 .6 .6]); plotm(sunLocation(1), sunLocation(2),'ko','MarkerFaceColor','y') plotm(grayLine(:,1), grayLine(:,2),'red') plotm(natick(1),natick(2),'ro','MarkerFaceColor','r')
And there’s the answer we’ve been after. Tomorrow, in the unlikely event I get up early enough to see the sun rise, I’ll be sharing that experience and that instant with people in Arroyos De Mantua, Cuba and Ísafjörður, Iceland. There’s something comforting about that.
Published with MATLAB® R2013b
]]>I'd like to welcome back my fellow MATLAB Central blogger Brett Shoelson for the first in a three-part series on extracting curve values from a plot. You can find Brett over at the File Exchange Pick of the Week blog, or you can check out his many File Exchange contributions. -Steve
In my role as an application engineer, I frequently have the opportunity to help customers with their image processing-related problems. In fact, this is one of the aspects of my job that I enjoy the most; each customer's image poses unique--and sometimes difficult--challenges. I get to help solve them, and move on to other things.
Recently, a customer shared with me an image of a chart showing, among other things, the efficiency of a pump as a function of flow rate. He asked if I could help him automate the extraction of real data from this chart; he has to do this over and over, and is currently following a very manual workflow--clicking along the desired curve and recording the positions of his mouseclicks. If you're ever faced with a similar challenge and you can find a way to get the original data from which the chart was created, that's not only significantly easier, but less noisy than extracting the data from an image. But in this case, the original efficiency-versus-flow data were not available, and the customer had no alternative.
url = 'http://blogs.mathworks.com/images/steve/2013/example-pump-perf-curve.jpg';
chart = imread(url);
imshow(chart)
There are three axes in this image. In the top-most of the axes (which is also the largest), there are curves representing both "Head" (in m) and "Efficiency" (in percent) versus flow rate (m^3/h). (The efficiency is labeled on the right side of the top axis.) The customer was interested only in the data reflected by the "efficiency" curve, as labeled in the image.
What information can we rely on from chart to chart? After asking a few questions, I learned that the efficiency curve always starts in the lower left corner of the top axis. Moreover, the axes and curves in this graph are similar--in form, shape, and color--to the corresponding curves in the many other graphs from which he needs to extract these data. (Otherwise, automating this data extraction would be impossible, and we'd have to try for a more "semi-automated" approach to facilitate manual analyses.)
I like to make the point that solutions to image processing problems are typically non-unique. Your own approach may be different--perhaps radically so--from my own. In fact, I often reflect on problems that I've already solved, and realize that I might solve it differently if I were to face it again.
Before you read my solution, take a moment to look at the customer's image and think about how you would approach the problem. When I'm done sharing my workflow, I will encourage you to comment on it, or even to improve upon it. That way, we can learn from each other!
I'm going to step through the steps I took to automate this data-extraction process. This was somewhat tricky, so to keep this blog post from being too long, I'm going to break it up into three parts:
The information we're after (the "signal") is represented graphically in the "Efficiency" curve in the top axes. The rest of the curves are just "noise" here, and our initial task is to ascertain how to discard all of those noise sources, while keeping the signal intact.
Recognize that each axis in the original image is delineated by a black border (though that's difficult to see here, with the screen capture). With little effort, we can convert those black borders to white; that will allow us to create masks of each axis, and to select (or discard) specific ones.
imshow(chart)
title('Original RGB')
Let's create a reversed, binary version of the original; this will cast the black borders to white:
bw = im2bw(imcomplement(rgb2gray(chart))); % I combined some steps here imshow(bw) title('Grayscaled/Reversed/Binarized')
Now that the axes are bordered in white, we can fill the image regions to create solid regions of interest:
filled = imfill(bw,'holes'); imshow(filled) title('All Regions Filled')
Since I know that the axis I want to keep is the largest one. I can use that information to create a mask that will be useful for discarding a lot of the noise in the original image:
cc = bwconncomp(filled); stats = regionprops(cc,'Area'); A = [stats.Area]; [~,biggest] = max(A); filled(labelmatrix(cc)~=biggest) = 0; imshow(filled); title('Axis mask') bb = regionprops(filled,'BoundingBox'); %We'll use this later! bb = bb.BoundingBox;
(Note that the selection of the region [axes] of interest [axes 4, bottom right] was achieved simply by applying a logical constraint, and that I could easily have elected to use any combination of logical conditions to select the desired axes. For instance, I could have calculated both areas and centroids of the "blobs" in the "All-Regions-Filled" black and white image. Then I could have picked the region that is larger than a threshold area value, and that has a y-centroid above those of the other axes.)
Now we can use this mask to pare the original color image, plane-by-plane:
chart(~repmat(filled,[1 1 3])) = 0;
imshow(chart)
title('Masked RGB Original')
We've now discarded a lot of the noise, but the difficult part remains. Unfortunately, several non-interest curves intersect the curve we need. Even more unfortunately, one of those curves is the same color as the curve of interest! But perhaps we can continue along the path of noise reduction by exploring and using the color information in the image; it's possible that something in a grayscale representation of the image (rgb2gray, or individual red, green, or blue colorplanes) could be useful here. To explore that possibility, I like to use ExploreRGB.
ExploreRGB(chart)
Nothing that's immediately obvious jumps out at me that will make the selection of the efficiency curve any easier. But let's click on the "Advanced-Mode" toolbar button (upper left) to explore different colorspace representations:
Now, by studying the "blue chrominance image" (in the 4th row, 2nd column), we can see that the blue curves can readily differentiated from non-blue curves simply by changing colorspaces.
(Recognize that I could simply click on the blue-chrominance image, and then right-click to export that image to the base workspace, but I'll show the equivalent code here:) Convert to YCbCr, extract blue chrominance image:
curveOfInterest = rgb2ycbcr(chart);
curveOfInterest = curveOfInterest(:,:,2);
imshow(curveOfInterest);
title('Cb/masked');
By the way, I didn't have to create handles to the subplot axes (ax(1)--ax(4), above). But doing so allows me to use expandAxes, which conveniently lets me click on any of those axes to expand it to full screen, and then to right-click to export the image it contains to my MATLAB workspace.
Let's stop there. Next week, we'll resume the problem and isolate only the efficiency curve. (That, it turns out, is the most difficult part of this problem!)
Get
the MATLAB code
Published with MATLAB® R2013b
In today's post, I am joined by Mike Hosea, a developer who occasionally works on integration routines for MATLAB. In recent releases, we added new integration routines, including integral2 for double integrals. Today we'll talk about some nuances in using this routine and adjusting the absolute and relative tolerances.
There is at least one requirement on the integrand when you call any integration function, and that is that it can be called with multiple inputs and it will return corresponding outputs, 1-for-1, with the inputs. So, to integrate $x * y$, you need to supply the vectorized form, i.e., x .* y.
Instead of creating a file with the integrand definition, we instead create an anonymous function. If you are unfamiliar with anonymous functions and want to learn more about them, please check the categories to the right of the blog for discussions and examples.
F0 = @(x,y) x .* y;
We can now integrate F between 0 and 1 for both directions, x and y.
Q = integral2(F0,0,1,0,1)
Q = 0.250000000000000
Here is another example. First let's integrate the function F1
F = @(x,y)exp(10*x).*tan(y)
F = @(x,y)exp(10*x).*tan(y)
between 0 and 1 for x and 0 and $\pi/4$ for y.
Q1 = integral2(F,0,1,0,pi/4)
Q1 = 7.633444730985692e+02
The analytic answer for this integral is
Q0 = log(sqrt(2))/10*(exp(10) - 1)
Q0 = 7.633444758094897e+02
and you can see that Q1 differs from Q0, perhaps more than we might want. We first adjust the numeric display to see enough decimal digits, before displaying the abolute error of the integration Q1.
format long
intabserr1 = abs(Q1-Q0)
intabserr1 = 2.710920512072335e-06
One of the advantages of the integral2 routine (and its companions for other dimensions) over most of the older routines is that they provide mixed relative and absolute error control, which is to say, you can supply two tolerances instead of just one. Not only can you, you should. You actually need both, or rather sometimes the one and sometimes the other, but it is good practice to think about what you want for both of them and to be aware of the default values in case you choose not to set them. What you would rarely want to do is change only the absolute error tolerance, AbsTol.
Now suppose we want more accuracy and try to crank down the absolute tolerance.
Q2 = integral2(F,0,1,0,pi/4,'AbsTol',1e-14)
intabserr2 = abs(Q2 - Q0)
Q2 = 7.633444730985692e+02 intabserr2 = 2.710920512072335e-06
We got the same answer as before! The reason is that integral2 uses a mixture of relative and absolute error control, and it is always the least restrictive of the two that matters. The default AbsTol is 1e-10, and the default RelTol is 1e-6. The absolute error is bigger than the absolute tolerance, true, but the relative error is
intrelerr2 = abs(Q2 - Q0)/abs(Q0)
intrelerr2 = 3.551372411777180e-09
and that's smaller than RelTol, so integral2 returned that answer. Since the relative error is between 1e-8 and 1e-9, there should be about 8 correct significant digits, and indeed, we have that.
disp([Q0;Q2])
1.0e+02 * 7.633444758094897 7.633444730985692
The default RelTol is 1e-6, so let's try making that smaller, instead of AbsTol, to get more accuracy.
Q3 = integral2(F,0,1,0,pi/4,'RelTol',1e-12)
intrelerr3 = abs(Q3 - Q0)/abs(Q0)
Q3 = 7.633444758094937e+02 intrelerr3 = 5.212639177138190e-15
Now since the relative error is between 1e-14 and 1e-15, we should have about 14 significant digits correct, and we do
disp([Q0;Q3])
1.0e+02 * 7.633444758094897 7.633444758094937
Usually RelTol is the more important tolerance to manipulate, and the more useful one, because it indirectly controls the number of correct significant digits. But there is always a catch. Suppose the correct answer is zero. Then relative error is undefined! That's why integral2 has to use the least restrictive of AbsTol and RelTol. The actual error condition it tries to meet is
ErrEst <= max(AbsTol,RelTol*abs(Q))
where ErrEst is an approximate upper bound on the absolute error.
It turns out that the ratio R = AbsTol/RelTol is a cutoff. When abs(Q) < R, AbsTol is the tolerance that matters. When abs(Q) > R, RelTol is the one that matters. So it is good practice to set both tolerances to something reasonable by asking yourself two questions:
So when you choose your own tolerances, the calls should look something like
Q4 = integral2(F,0,1,0,pi/4,'AbsTol',1e-12,'RelTol',1e-10)
Q4 = 7.633444758095101e+02
In this post, we've only covered the topic of setting tolerances to insure a reasonable estimate for your integrals. We plan to talk more about which method to choose, order of integration (e.g., x then y), and dealing with singularities and discontinuities.
In the meantime, do you have any experiences to share about performing your own numerical integrations with MATLAB? Let us know here.
Get
the MATLAB code
Published with MATLAB® R2013b
The Fiedler companion matrix distributes the coefficients of a polynomial along the diagonals of an elegant pentadiagonal matrix whose eigenvalues are equal to the zeros of the polynomial.
I spent a very pleasant day in Berkeley several weeks ago with my friend Beresford Parlett. During the visit he told me about an elegant version of the companion matrix that Miroslav Fiedler introduced in 2003. Fiedler is Czech mathematician known for his contributions to linear algebra, graph theory, and the connections between the two. He has been a frequent participant in the Householder symposia and, at various times, a member of its organizing committee.
MATLAB usually represents a polynomial of degree $n$ by a vector with $n+1$ components. With such a representation, the leading term has its own coefficient, which might be zero. In this article, I want to have monic polynomials, where the leading coefficient is one, and so there are just $n$ components in the corresponding vector.
$$ p(x) = x^n + a_1 x^{n-1} + a_2 ^{n-2} + ... + a_{n-1} x + a_n $$
Associated with every such polynomial is its companion matrix with the coefficients strung out along the first row and ones along the first subdiagonal. This traditional companion matrix is readily generated by a slick bit of MATLAB code.
type companion
function C = companion(a) n = length(a); C = [-a; eye(n-1,n)];
Here is the traditional companion matrix for the coefficients 1:10.
companion(1:10)
ans = -1 -2 -3 -4 -5 -6 -7 -8 -9 -10 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0
And here is the spy plot for 30 coefficients.
spy(companion(1:30))
Fiedler showed that we could generate a different companion matrix by arranging the coefficients, alternating with zeros, along the super- and subdiagonal, together with ones and zeros along the supersuper- and subsubdiagonal of a pentadiagonal matrix. I love the MATLAB code that generates the Fiedler companion matrix.
type fiedler
function F = fiedler(a) n = length(a); b = ones(n-2,1); b(1:2:n-2) = 0; c = -a(2:n); c(1:2:n-1) = 0; c(1) = 1; d = -a(2:n); d(2:2:n-1) = 0; e = ones(n-2,1); e(2:2:n-2) = 0; F = diag(b,-2) + diag(c,-1) + diag(d,1) + diag(e,2); F(1,1) = -a(1);
Here is the Fiedler companion matrix for the coefficients 1:10.
fiedler(1:10)
ans = -1 -2 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 -3 0 -4 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -5 0 -6 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -7 0 -8 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -9 0 -10 0 0 0 0 0 0 0 1 0 0
And here is the spy plot for 30 coefficients.
spy(fiedler(1:30))
When I wrote the first Fortran MATLAB in the late 1970s, I included the roots function. I will admit now that it was hardly more than a quick hack. roots finds the zeros of a polynomial by setting up the (traditional) companion matrix and calling the eig function to compute its eigenvalues. I included this function for two reasons. First, it gave me a polynomial zero finder at very little cost in the way of additional code. And second, it turned the tables on the linear algebra textbook approach of finding matrix eigenvalues by computing the zeros of the characteristic polynomial. But I wasn't sure that roots was always accurate.
I was lucky. roots turned out to be a good idea. Papers published 15 years later by Kim-Chuan Toh and Nick Trefethen, and by Alan Edelman and H. Murakami, established that, as long as balancing is used, roots is numerically reliable.
On the other hand, roots can be criticized from an efficiency point of view. A polynomial has $O(n)$ coefficients. I'm not sure what the computational complexity of the zero finding problem should be, but the complexity of the matrix eigenvalue problem is $O(n^2)$ storage and $O(n^3)$ time, which is certainly too much.
This is where Fiedler comes in. Fiedler can replace the traditional companion matrix in roots. Now, if we could find a matrix eigenvalue algorithm that works on nonsymmetric pentadiagonal matrices in $O(n)$ storage and $O(n^2)$ time, then we would have the ultimate roots.
For an example, let's use the Wilkinson polynomial with roots 1 through 10.
syms x
p = prod(x-(1:10))
p = expand(p)
a = coeffs(p);
F = fiedler(double(fliplr(a(1:end-1))))
p = (x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 8)*(x - 9)*(x - 10) p = x^10 - 55*x^9 + 1320*x^8 - 18150*x^7 + 157773*x^6 - 902055*x^5 + 3416930*x^4 - 8409500*x^3 + 12753576*x^2 - 10628640*x + 3628800 F = Columns 1 through 6 55 -1320 1 0 0 0 1 0 0 0 0 0 0 18150 0 -157773 1 0 0 1 0 0 0 0 0 0 0 902055 0 -3416930 0 0 0 1 0 0 0 0 0 0 0 8409500 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 Columns 7 through 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 -12753576 1 0 0 0 0 0 0 10628640 0 -3628800 0 1 0 0
Look at those coefficients. Can you imagine what contortions a vector $v$ must go through in order to satisfy $Fv = \lambda v$ for some integer $\lambda$ between 1 and 10? As a result, the computed eigenvalues of $F$ show a loss of accuracy caused by the sensitivity of the roots of this polynomial.
format long
lambda = eig(F)
lambda = 9.999999999992546 8.999999999933312 8.000000000317854 6.999999999479941 6.000000000423509 4.999999999814492 4.000000000042952 2.999999999995175 2.000000000000225 0.999999999999996
Symmetric matrices have perfectly well conditioned eigenvalues. Balancing, an idea introduced by Parlett and Reinsch, is the attempt to find a diagonal similarity transformation that makes a nonsymmetric matrix "more nearly symmetric", in the sense that the norms of its rows and the norms of its columns are closer to being equal. Balancing is not always possible. For example, an upper or lower triangular matric cannot be balanced. Triangular matrices are hopelessly nonsymmetric, but that's OK because their eigenvalues are on the diagonal.
By default, MATLAB, using LAPACK, balances nonsymmetric matrices before computing eigenvalues. As Toh and Trefethen, and Edelman and Murakami, emphasized, this is essential for the numerical stability of roots. The same is true for roots based on a Fiedler matrix. When we called eig above, MATLAB balanced the input matrix. It helps a lot.
Balancing is done by powers of two, so as not to introduce any roundoff errors. The large elements are scaled down significantly, while the ones on the outer diagonals are scaled up by only a few powers of two. Here is the balanced Fiedler matrix.
format short
B = balance(F)
B = Columns 1 through 7 55.0000 -20.6250 32.0000 0 0 0 0 64.0000 0 0 0 0 0 0 0 8.8623 0 -4.8148 8.0000 0 0 0 16.0000 0 0 0 0 0 0 0 0 3.4411 0 -3.2586 4.0000 0 0 0 4.0000 0 0 0 0 0 0 0 0 2.0050 0 0 0 0 0 0 2.0000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Columns 8 through 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1.5203 2.0000 0 0 0 0 0.6335 0 -0.4326 0.5000 0 0
Balancing decreases the norm of the matrix by five orders of magnitude. Here is the before and after.
format short g disp(' no balance balance') disp([norm(F) norm(B)])
no balance balance 1.8092e+07 88.433
Sometimes, but not often, it is desirable to bypass balancing. This is the situation when the matrix contains small elements that are the result of roundoff or other errors that should not be scaled up to become comparable with the other elements. There are not any such elements in this matrix.
Here are the computed eigenvalues, without and with balancing.
format long disp(' no balance balance') disp([eig(F,'nobalance') eig(F)])
no balance balance 9.999999330734642 9.999999999992546 8.999996381627081 8.999999999933312 8.000028620785038 8.000000000317854 6.999932559108236 6.999999999479941 6.000079411601993 6.000000000423509 4.999948356986658 4.999999999814492 4.000018256147182 4.000000000042952 2.999996926380218 2.999999999995175 2.000000155398907 2.000000000000225 1.000000001522463 0.999999999999996
You have to look at the trailing digits to see that the eigenvalues computed without balancing are less accurate. It becomes clearer when we compute the errors. Balancing improves the accuracy by the five orders of magnitude that is achieved by reducing the norm of the input.
format short e disp(' balance no balance') disp([eig(F)-(10:-1:1)' eig(F,'nobalance')-(10:-1:1)'])
balance no balance -7.4536e-12 -6.6927e-07 -6.6688e-11 -3.6184e-06 3.1785e-10 2.8621e-05 -5.2006e-10 -6.7441e-05 4.2351e-10 7.9412e-05 -1.8551e-10 -5.1643e-05 4.2952e-11 1.8256e-05 -4.8255e-12 -3.0736e-06 2.2515e-13 1.5540e-07 -3.7748e-15 1.5225e-09
The traditional companion matrix and the Fiedler companion matrix are similar to each other. Can you find the similarity transformation? It has a very interesting nonzero structure.
If we could figure out how to take advantage of the pentadiagonal structure, Miraslav Fiedler's elegant companion matrix would be an all-around winner, but we must not forget to balance.
Miroslav Fiedler, A note on companion matrices, Linear Algebra and its Applications 372 (2003), 325-331, <http://citeseerx.ist.psu.edu/viewdoc/download>
Alan Edelman and H. Murakami, Polynomial Roots from Companion Matrix Eigenvalues, Mathematics of Computation. 64 (1995), 763--776, <http://www-math.mit.edu/~edelman/homepage/papers/companion.pdf>
Kim-Chuan Toh and Lloyd N. Trefethen, Pseudozeros of polynomials and pseudospectra of companion matrices, Numer. Math. 68: (1994), 403-425, <http://people.maths.ox.ac.uk/~trefethen/publication/PDF/1994_61.pdf>
Cleve Moler, Cleve's Corner: ROOTS -- Of polynomials, that is, <http://www.mathworks.com/company/newsletters/articles/roots-of-polynomials-that-is.html>
Beresford Parlett and Christian Reinsch, Balancing a matrix for calculation of eigenvalues and eigenvectors, Numer. Math., 13, (1969), 293-304, <http://link.springer.com/article/10.1007/BF02165404>
Get
the MATLAB code
Published with MATLAB® R2014a
Even though I am one of the founders of the MathWorks, I only acted as an advisor to the company for its first five years. During that time, from 1985 to 1989, I was trying my luck with two Silicon Valley computer startup companies. Both enterprises failed as businesses, but the experience taught me a great deal about the computer industry, and influenced how I viewed the eventual development of MATLAB. The second startup was known as Ardent Computer for most of its brief existence.
As soon as we could, we ported MATLAB to the Titan, and MathWorks established a working business relationship with Ardent Computer. On a trip to the West Coast, Jack Little and Jim Tung from MathWorks visited Ardent. In a meeting at Ardent, I gave Jack and Jim my Ardent business card and I gave the Ardent folks my MathWorks business card. The Ardent Titan eventually became the only computer in the world to ever have MATLAB delivered with its bundled software.
MATLAB and the Titan were excellent companions. The Titan's computational power came from its vector processor and, back then, MATLAB was still "Matrix Laboratory". We could demonstrate the Linpack Benchmark with one MATLAB command, or with one mouse click. The Linpack Benchmark provided the starting point of a deeper discussion of the capabilities of MATLAB.
MATLAB ran on top of Ardent's math libraries. Recall that Linpack was more than just a benchmark. Linpack and Eispack were the Fortran matrix libraries that got us started in this whole business. Randy Allen, who was one of Ken Kennedy's students, had come to Ardent to write vectorizing Fortran and C compilers. Vectorizing the Linpack code was easy because we had written it with vectoriziation in mind. Vectorizing the Eispack code is harder, both because the algorithms are more complex and because the coding style did not anticipate vector operations.
We also needed to have vectorized versions of the transcendental math functions. I had worked on the trig and exponential functions before Ardent because I had a research project at the University of New Mexico sponsored by IBM and a consulting gig at Convex that both involved vectorized math libraries. That work carried over readily to the MIPS chip that we were using on the Titan.
I had avoided the power function, $x^y$. That's harder because it involves two arguments. Which is the vector? I had rationalized not thinking about it saying by "Nobody ever wants to call the power function on a vector." But then one day we had a visit from a potential customer -- a large, international bank who shall remain nameless except to say their headquarters is in San Francisco. They were impressed with how easy it was to do their mortgage portfolio simulations in MATLAB, but we were concerned about performance. It turns out that most of the time was being spent in the power function on a vector. Fortunately, it was interest calculations of the form $s^{(n/365)}$ where $s$ is a scalar and $n$ is an integer vector containing numbers of days. That was easy to vectorize and add to the library.
Gustave Dore was a nineteenth century French engraver who illustrated works by Milton, Lord Byron, and Edgar Allen Poe, among others. Dore Avenue is an exit off the Bayshore Freeway on the way from Silicon Valley to the SFO Airport. When Michael Kaplan drove past Dore, the exit, it reminded him of Dore, the artist, and he decided to use the acronym for his Dynamic Object Rendering Environment.
I had met Michael in Oregon when he was trying to put Dore on the Intel Hypercube. That seemed like a good idea because ray tracing is computationally intensive and embarrassingly parallel. But the iPSC was pretty hopeless because, among other things, it didn't yet have a graphics output device. But Dore and the Titan were made for each other. The Titan had hardware support for all the tasks that Dore required, as well as a large frame buffer and a fast, high resolution display. Dore became the cornerstone of Ardent's graphics software environment.
In 1988 the graphics in MATLAB itself was still pretty basic. We did not yet have three dimensional color in our shipping product, and would not have it until 1992, with MATLAB version 4.0. So by linking MATLAB to Dore on the Titan we were able to anticipate capabilities that would not be available on other machines for another three or four years.
The September, 1988, issue of Computer, published by the IEEE Computer Society, included an extensive technical article by Ardent engineers about the Titan. The graphic on the cover shown above was chosen by the journal's editors. It is a generalization of the Moebius strip and was computed in MATLAB and rendered by Dore. Today we can create this graphic in MATLAB itself by applying shading and lighting to a surf plot
For me the even more exciting graphic is the vibrating membrane, an animation of the solution to the wave equation that underlies the MathWorks logo. This is my recreation in an animated .gif of something I first saw in Michael Kaplan's darkened office at Ardent sometime in 1988. It was the first time I had ever seen the membrane moving on a computer screen and it nearly brought tears to my eyes.
Much of what I know about computer graphics, I learned by accessing Dore from MATLAB while I was at Ardent. And many of the demos that Ardent used to show off the Titan at trade shows and customers visits involved MATLAB driving Dore. It was a terrific marriage -- while it lasted.
It took longer to bring the Titan to market than the original business plan called for. The machine was ultimately more expensive than the original business plan called for. It was bigger and heavier. It was noisier. And it required 220 volts.
It turned out to be harder to vectorize and much harder to parallelize most programs than everybody had hoped. The big-time applications on the big-time supercomputers like the Crays had been doing it for years, but most scientific and engineering applications that had been running on Vaxes and Sun workstations had not. We were beginning to see that Linpack performance was not so strongly correlated with the performance of more complicated programs as the machine architecture grew more complex. This is certainly true today.
Forest Baskett is a friend of mine. He has been at Stanford, at Los Alamos, and for a while he ran a Palo Alto research lab for DEC. He is now a venture capitalist in The Valley. But for a long time, including the time I was at Ardent, he was a VP at Silicon Graphics. In early 1989, he gave a seminar talk at Stanford about RISC with a deliberately provocative title like "Vector Computing is Dead". In the question period, I marched to the front of the seminar room and triumphantly produced my own overhead transparency with the latest performance numbers from the Titan. But Forest was right. We were the last of a dying breed.
But the biggest problem was competition. Direct competition from Stellar. And increasing competition from the workstations as Silicon Graphics and Sun increased their computing power. It took too long for customers to decide what to buy. Our sales people had to work too hard for each sale. Competitive benchmarks were often involved. These machines usually cost over \$100K and were not impulsive decisions.
On the East Coast, Stellar Computer was also having difficulties selling their machines. They did not have a vector machine, and I know less about the specifics of their frustrations, but I do know they were also suffering.
One morning in the fall of 1989, everybody at Ardent was summoned to an auditorium near Mountain View that could hold us all. We were connected via satellite to a similar auditorium near Waltham, Mass, that could hold everybody from Stellar. In that auditorium Al Michaels, CEO of Ardent, stood next to Bill Poduska, CEO of Stellar. I forget whether either was smiling. They announced that Stellar Computer and Ardent Computer had been merged to form Stardent Computer.
It was a shotgun marriage with the financial backers calling the shots. If there was any chance of saving their investments, it would have to be through combining the two companies.
Production and marketing of the first generation machines of both companies were terminated. Ardent had a second generation machine well along in development, and that was going to be Stardent's next product. It was an actual desk side workstation, using Intel's ill-fated excursion into RISC processing, the i860. There was no vector processor. There was considerable reuse of the successful graphics components from the Titan. This machine was to be known as the "Stiletto". I still have a hand-out from a preview -- an actual stiletto, a serious weapon with a long slender blade. I know we wanted to be viewed as cutting edge, but that was going too far. Fortunately, that machine never reached the marketplace.
From its beginning, Stardent had no chance. Key people began to leave almost immediately. Whatever management was left became concentrated in Massachusetts, but development was centered in California. The subsequent implosion and eventual failure was the largest in the history of the computer industry, at least before the dot com era. The combined companies had lost hundreds of millions of dollars of investment capital.
There is an interesting Wikipedia article about Stardent that is primarily about the period after I left. It is all about law suits and counter law suits. I know some of the players, but I have no personal knowledge of the events in that article.
The NA Digest is an electronic newsletter for the community of numerical analysis and mathematical software that has been going since the early 1980's. I was its editor for a long time, including the period around the breakup of Stardent. In late October, 1989, I ran this announcement. (Note the MathWorks phone number at the time.)
From: Cleve Moler <na.moler@na-net.stanford.edu> Date: Sun Oct 29 10:39:38 PST 1989 Subject: Positions at The MathWorks
The MathWorks is the company which develops and markets MATLAB. The company currently employs about 30 people and expects to add three or four more soon. The company headquarters is in East Natick, Massachusetts, which is about a half hour drive west of Boston.
The background and interests expected for the various positions available range from numerical linear algebra and matrix computation to systems programming and graphics. Educational level and experience expected range from inexperienced beginner willing to learn to seasoned Ph.D. with a personal library of M-files.
For more information, send email to na.moler@na-net.stanford.edu or phone me at 408-732-0400. Or, contact the MathWorks' president, John Little, with email to ******@mathworks.com, phone 508-653-1415, or write to: The MathWorks 21 Eliot Street South Natick, MA 01760
A couple of weeks after posting this announcement in the NA Digest, I decided to resign my position at Stardent and take one of the jobs at MathWorks that I had advertised. I've been a regular full time employee of MathWorks ever since.
Michael R. Kaplan, The Design of the Dore Graphics System, in E. H. Blake and Peter Wisskirchen, editors, Advances in Object-Oriented Graphics I, Springer, pp 177-198, 1991, <http://link.springer.com/chapter/10.1007%2F978-3-642-76303-8_10?>
Get
the MATLAB code
Published with MATLAB® R2013b
I'd like to introduce this week's guest blogger Edric Ellis. Edric works for the Parallel Computing development team here at The MathWorks. In this post he will talk about using the parfeval command in Parallel Computing Toolbox.
In release R2013b, we introduced a new function parfeval to Parallel Computing Toolbox which allows you to run MATLAB® functions on workers in a parallel pool without blocking your desktop MATLAB. This allows you to create parallel programs that are more loosely structured than those using either parfor or spmd. You can also have your desktop MATLAB process the results as they become available, allowing you to create more interactive parallel programs.
We're going to take a look at using parfeval to speed up accessing a web service by making multiple requests in parallel, and then displaying each result as soon as it is available. We're going to query the Flickr image sharing site for astronomical images that have structured data associated with them. We'll approach the problem in three phases:
Flickr provides a web-based API that allows you to query its large database of images. To use the API, you submit a structured query to their servers, and an XML document is returned to you. Retrieving the resulting image data from Flickr takes some time, and this can run more quickly by performing the retrievals in parallel. This is an example of an "I/O-bound" operation.
The Flickr web API is fully documented here. To use the API, you need a Flickr account, and you need to request an 'API Key' which is included in each request to the API. Once you've got your own key, you can place it in a function called flickrApiKey which simply returns the key as a string.
appKey = flickrApiKey;
We're going to search for images which have a specific tag: astrometrydotnet:status=solved. These images are submitted to a Flickr group, and then astrometry.net analyses and annotates the images. We build a search url, and then let xmlread read the results from that URL.
url = ['http://api.flickr.com/services/rest/', ... % API Base URL '?method=flickr.photos.search', ... % API Method '&api_key=', appKey, ... % Our API key '&tags=astrometrydotnet%3Astatus%3Dsolved', ... % Tag to search '&license=4', ... % Creative Commons '&sort=interestingness-desc', ... % Sort order '&per_page=12', ... % Limit results to 12 '&media=photos']; % Only return photos response = xmlread(url); % Get all the 'photo' elements from the document. photos = response.getDocumentElement.getElementsByTagName('photo'); numPhotos = getLength(photos); % Ensure our search query return some results. if numPhotos == 0 error('Failed to retrieve photos from Flickr. XML response was: %s', ... xmlwrite(response)); end % Otherwise, convert the search results into a cell array of structures % containing the search result information. photoInfo = cell(1, numPhotos); for idx = 1:numPhotos node = photos.item(idx-1); photoInfo{idx} = struct('farm', char(node.getAttribute('farm')), ... 'server', char(node.getAttribute('server')), ... 'id', char(node.getAttribute('id')), ... 'secret', char(node.getAttribute('secret')), ... 'owner', char(node.getAttribute('owner'))); end
For each search result, we need to run a function to retrieve the image data. The function also searches through the Flickr comments for the structured information from astrometry.net, and returns the direction the telescope was pointing. In the next stage, we'll submit multiple invocations of this function using parfeval.
function [photo, location, username] = getAstrometryPhotoFromFlickr(appKey, info) % First, build up the photo URL as per the Flickr documentation % here: <http://www.flickr.com/services/api/misc.urls.html> url = sprintf('http://farm%s.staticflickr.com/%s/%s_%s.jpg', ... info.farm, info.server, info.id, info.secret); try photo = imread(url); catch E % Sometimes the read fails. Try one more time. photo = imread(url); end % % Get the photo info to extract the username url = ['http://api.flickr.com/services/rest/', ... '?method=flickr.photos.getInfo', ... '&api_key=', appKey, '&photo_id=', info.id]; response = xmlread(url); owner = response.getDocumentElement.getElementsByTagName('owner'); if getLength(owner) == 1 username = char(owner.item(0).getAttribute('username')); else username = 'Unknown'; end % % Next, look through the comments for the photo to extract % the annotation from astrometry.net. url = ['http://api.flickr.com/services/rest/', ... '?method=flickr.photos.comments.getList', ... '&api_key=', appKey, '&photo_id=', info.id]; response = xmlread(url); % % Get all the actual comment elements from the XML response. comments = response.getDocumentElement.getElementsByTagName('comment'); % % Loop over all the comments, looking for the first by astrometry.net % which contains information about the photo. for idx = 0:(getLength(comments)-1) comment = comments.item(idx); if strcmp(char(comment.getAttribute('authorname')), 'astrometry.net') % We've found the comment, extract the comment text. commentText = char(comment.getTextContent()); % Pick out the location information into a row vector by % using a regular expression. locationText = regexprep(... commentText, '.*center: *\(([-0-9,. ]+)\) *degrees.*', '$1'); location = sscanf(locationText, '%g,', [1, 2]); return end end % % We didn't find the astrometry.net comment, so location is unknown. location = [NaN, NaN]; end
For each search result, we will make a parfeval request for that particular result. This will cause one of the workers in our parallel pool to invoke that function. If a parallel pool is not already open, then one will be created automatically. Each parfeval call returns a parallel.Future instance. The parfeval requests are executed in order by the parallel pool workers. The syntax of parfeval is rather similar to that of the MATLAB function feval, except that because the evaluation doesn't take place immediately, you must specify how many output arguments you want to request. In this case, we want 3 outputs from getAstrometryPhotoFromFlickr.
for idx = 1:numPhotos futures(idx) = parfeval(@getAstrometryPhotoFromFlickr, 3, ... appKey, photoInfo{idx}); end
Starting parallel pool (parpool) using the 'local' profile ... connected to 6 workers.
While the workers are busy retrieving images from Flickr, we can set up a figure to display those images.
% Create a new figure, initially not visible. figHandle = figure('Position', [200, 200, 600, 800], ... 'Units', 'pixels', 'Visible', 'off'); % Use FEX submission 'tight_subplot' to set up the axes % http://www.mathworks.com/matlabcentral/fileexchange/27991 axesHandles = tight_subplot(4, 3, [0.06, 0.01], [0.01, 0.06], 0.01); axis(axesHandles, 'off');
We now enter a loop where we wait for results to become available on the workers. We use fetchNext to achieve this. We pass fetchNext our array futures, and it returns when another element has completed and a new result is available. We can also specify a maximum time to wait for each new result, so that if the web service takes a really long time to run, we can abort execution gracefully.
Note that the results from fetchNext can arrive in any order, and we use the first output argument originalIdx to work out which element of futures the result corresponds to.
We're going to specify an upper bound on how long we're going to wait for all the futures to complete. If we exceed that time limit and there are still futures running, we can cancel them.
overallTimeLimit = 10; % seconds t = tic; set(figHandle, 'Visible', 'on'); numCompleted = 0; while numCompleted < numPhotos try % This blocks for up to 1 second waiting for a result. [originalIdx, photo, location, user] = fetchNext(futures, 1); catch E % Sometimes, the URL simply cannot be read, so an % error is thrown by the worker. Let's display that % and carry on. warning('Failed to read an image: %s', getReport(E)); originalIdx = []; end % If 'fetchNext' completed successfully, originalIdx % will be non-empty, and will contain the index into 'futures' % corresponding to the work that has just finished. if ~isempty(originalIdx) % Display attribution for photo and link to the original fprintf('Photo %2d by Flickr user: %s\n', originalIdx, user); info = photoInfo{originalIdx}; fprintf('Original at: http://www.flickr.com/%s/%s/\n', info.owner, info.id); % Pick our axes: axesToUse = axesHandles(originalIdx); % Display the image image(photo, 'Parent', axesToUse); axis(axesToUse, 'off'); axis(axesToUse, 'equal'); % Set the title of the axis to be the location title(axesToUse, num2str(location)); numCompleted = numCompleted + 1; elseif toc(t) > overallTimeLimit % We have exceeded our time budget! disp('Time limit expired!'); % Cancelling the futures stops execution of any running % futures, but has no effect on already-completed futures. cancel(futures); break; end end toc(t);
Photo 6 by Flickr user: s58y Original at: http://www.flickr.com/45032885@N04/4145493112/ Photo 2 by Flickr user: davedehetre Original at: http://www.flickr.com/22433418@N04/4950431138/ Photo 5 by Flickr user: s58y Original at: http://www.flickr.com/45032885@N04/8541349810/ Photo 1 by Flickr user: s58y Original at: http://www.flickr.com/45032885@N04/4553112538/ Photo 4 by Flickr user: Torben Bjørn Hansen Original at: http://www.flickr.com/21067003@N07/6105409913/ Photo 3 by Flickr user: Harry Thomas Photography Original at: http://www.flickr.com/84598277@N08/8882356507/ Photo 7 by Flickr user: s58y Original at: http://www.flickr.com/45032885@N04/4144380775/ Photo 8 by Flickr user: Jody Roberts Original at: http://www.flickr.com/8374715@N06/8719966332/ Photo 9 by Flickr user: s58y Original at: http://www.flickr.com/45032885@N04/4145140448/ Photo 11 by Flickr user: davedehetre Original at: http://www.flickr.com/22433418@N04/4954464378/ Photo 12 by Flickr user: s58y Original at: http://www.flickr.com/45032885@N04/4535865643/ Photo 10 by Flickr user: cfaobam Original at: http://www.flickr.com/33763963@N05/7060054229/ Elapsed time is 2.897368 seconds.
We have seen how using the parfeval command can be used to retrieve data from a web service, and display the results interactively as they become available. parfeval allows you more flexibility than parfor to:
If you've ever wanted to show a waitbar or plot partial results during a parfor loop, then you might be able to adapt your code to use parfeval. If you have any further ideas, let us know here.
Get
the MATLAB code
Published with MATLAB® R2013b
Last month, MathWorks held MATLAB EXPO UK at the Silverstone Circuit, the home of British motor racing. It was a great opportunity for the MATLAB and Simulink community to come together and the event provided something for everyone across multiple industries, applications and level of experience. The background roar of cars racing round the track certainly added to the atmosphere of the event.
A photo montage from MATLAB EXPO
The breadth of use of MATLAB in industry and academia never fails to impress me. This event featured top customer speakers from the aerospace, automotive, power, finance, and education industries:
Rolls-Royce takes the stage
The range of applications and technology showcased at the event was simply fantastic. Some of my personal favourites were:
A NAO Humanoid Robot takes a break
This year we added some introductory sessions for delegates wanting a primer on MathWorks core technology. These sessions proved very popular and we plan to expand them for next year. There was plenty for the advanced user too with a great ‘What’s new in MATLAB?’ talk from Chief Scientist Joe Hicklin, and masterclass tracks for both MATLAB and Simulink.
An animated Chief Scientist makes his point
There was also plenty of geeky fun on offer at MATLAB EXPO through games and competitions. In particular,
Concentrating hard in the F1 Simulator
The presentations from MATLAB EXPO are now available for download from the website:
http://www.matlabexpo.com/uk/
I hope to see some of you at next year’s community event!
]]>To see this in practice, type two or more characters in the command window:
Suggestions for MATLAB Functions (iPhone)
Tap on an option to complete it to the window.
You can leverage autocomplete for a variety of tasks, like navigating folder hierarchy:
Navigating Folder Hierarchy (Android)
Or displaying the methods that an object exposes:
Suggestions for an Object’s Methods (iPhone)
Autocomplete is available on the latest releases of MATLAB Mobile on the App Store and Google Play.
Has autocomplete improved your experience? Let us know by leaving a comment below.
]]>In the past MATLAB users had to go to two different places to find questions and answers. MathWorks-authored technical support solutions were located in the support area on the MathWorks website, while community questions and answers were located in Answers on MATLAB Central. Now both are located in MATLAB Answers. This creates a one-stop shop for MATLAB users to find answers to their questions, making the process faster and easier.
At first glance, these appear to be questions and answers with the same contributor. In reality these Support Answers are questions that were previously asked to MathWorks Technical Support by customers. The “MathWorks Support Team” contributor name will appear as both the questioner and answerer, even though the question was originally asked by someone else.
A lot of Support Answers were recently added to MATLAB Answers, and we want to make sure it is still easy for Community members to find what they are looking for. There is a toggle filter located in the left navigation panel titled “Refine by Source.” Clicking on one of the two options will make search results display either Community or MathWorks Support contributions.
Another way to differentiate between Community and MathWorks Support content is this little blue MathWorks icon located in the top right corner of all MathWorks Support contributions:
Below we show an example of a simple product download question that MATLAB users have asked. The three most relevant solutions found in MATLAB Answers are shown for “How to download older version of MATLAB”
As you can see, the first and third questions with accepted answers is from the MathWorks Support Team, and the second question is from someone within the MATLAB Central Community. If you drill into the top question the next screen will show the question previously asked by a MATLAB user, and the answer MathWorks Support provided the user. Any additional Community comments and answers will also be visible on this page.
We think this will make MATLAB Answers even more useful as a resource for MathWorks-related inquiries.
]]>For the longest time, and dating back to the day Cody launched, Bryant Tran (@bmtran) was the undisputed king of the hill. We even featured him in an interview on this blog. Eventually Richard Zapor was to surpass him. Richard is a tireless competitor and has made something of a specialty of creating high-quality challenges for the community (173 and counting!).
If you look at the blue line that represents Richard, you can see that he has steadily moved upward. But now look at Alfonso Nieto-Castanon. His green line progresses in enormous leaps punctuated by periods of relative inactivity. Alfonso finally passed Bryant Tran earlier this year, and today I was surprised to see that he was also converging rapidly with Richard at the number one spot.
When I saw this graph, I dashed straight over to the Cody leaderboard and what should I see but an exact tie for first place between Richard and Alfonso!
I was so impressed I took this snapshot. It’s probably changed by now and who knows who will be in the lead at the end of the day? But still, this snapshot depicts a great moment in Cody and sporting history.
Well done, gentlemen! It’s an honor to compete in the same game, though not in the same league.
]]>