MATLAB Central BlogsPipes Output
http://pipes.yahoo.com/pipes/pipe.info?_id=e72c9f53f10ba5cad505dbbed8d501cb
Wed, 01 Apr 2015 09:23:08 +0000http://pipes.yahoo.com/pipes/Face Coder Product Preview
http://feedproxy.google.com/~r/mathworks/desktop/~3/5CEsQwAPQUs/
<p>We don’t ordinarily talk about features under development, but today’s post comes from one of our internal product teams providing a very special sneak peak at an upcoming MATLAB feature.
There are 42 muscles in the human face. When different combinations of these muscles are flexed to varying degrees, a... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/community/2015/04/01/face-coder-product-preview/">read more >></a></p>http://blogs.mathworks.com/community/?p=3182Wed, 01 Apr 2015 07:00:59 +0000We don’t ordinarily talk about features under development, but today’s post comes from one of our internal product teams providing a very special sneak peak at an upcoming MATLAB feature.

There are 42 muscles in the human face. When different combinations of these muscles are flexed to varying degrees, a single face can take on millions of different positions. While even the most skilled experts can’t discern what someone is thinking based on their facial expressions, computers can build models based on the aggregation of billions of data points.

MATLAB already has a host of tools to build machine learning models and process big data. It also has the ability to identify and track human faces and facial features. Bringing all these technologies together, we’ve built a revolutionary tool that will change the way you use MATLAB. It’s called Face Coder, and early trials have shown promising results.

To use Face Coder, all you have to do is set up a webcam and look at your computer screen. Instead of having to type all your code out by hand, a slow and error-prone process, MATLAB will interpret your facial expressions and automatically generate the code you intend to write.

Most of our Beta feedback has been very positive:

“Face Coder has changed the way I write code. I can stare at my screen for hours thinking about MATLAB and never have to wait for my fingers to catch up to my thoughts.”

“With Face Coder, I feel a much closer connection to my computer. It’s almost as if MATLAB can read my mind.”

“I hope Face Coder will come out on MATLAB Mobile. I love to code, but I also love to maintain an active lifestyle. My dream is to someday be able to write code while I’m out cycling or kayaking.”

But other feedback has shown that we still have some work left to do:

“I often get distracted while I’m working. With Face Coder, I sometimes refocus on my work only to find that MATLAB has typed out my thoughts about bills I need to pay or what I want to eat for lunch.”

“Much of my work is classified, and it would be dangerous if unauthorized persons gained access to my MATLAB code. Therefore, I’m concerned that this feature will be a security vulnerability. I’ll need to wear sunglasses or a mask while programming to prevent anyone from recording and then interpreting my face.”

Though Face Coder is not yet available, all the technologies it’s built on are available in current MATLAB tools. These include:

]]>UncategorizedExperiencing the MATLAB Watch
http://feedproxy.google.com/~r/mathworks/moler/~3/vg5kYoIQ2ew/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/cleve/watch_blog_01.png"/></div><p>My experience with the forthcoming MATLAB watch exceeds my most optimistic expectations. The watch has not yet been announced officially, but a few of us have been testing prototypes for the past several weeks.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2015/04/01/experiencing-the-matlab-watch/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1183Wed, 01 Apr 2015 05:01:50 +0000

My experience with the forthcoming MATLAB watch exceeds my most optimistic expectations. The watch has not yet been announced officially, but a few of us have been testing prototypes for the past several weeks.

The watch features the latest mobile technology. The microprocessor is a low power dual core 64 bit Intel chip, running at 1.1 GHz, with IEEE 754 floating point arithmetic. There are 2 gigabytes of RAM, 250 gigabytes of SSD, and Wi-Fi and LTE wireless connections to the Internet. The primary battery lasts about 6 to 8 hours.

The unique new technology is in the watch band. The large piece on the underside of the band is a secondary battery that provides up to 24 additional hours of operation. (It requires about an hour to charge.) The links in the band are memory modules, similar to photo camera memory cards. Each of the 18 links holds 32 gigabytes, so a full watch band provides over half a terabyte of secondary storage.

The watch runs release R2015a of MathWorks software, including MATLAB version 8.5. It is possible to install Simulink and all of the Toolboxes, except those that rely upon connections to external hardware.

Commander

At first, as a dedicated user of the MATLAB command window, I was concerned about the lack of a full sized keyboard. I can't even effectively use the on-screen keyboard on my cell phone, so the tiny one on the watch is very frustrating. But the voice recognition command window, which is known as "Commander", is a joy to use. Since Commander only has to deal with a limited vocabulary, she is much more reliable than general purpose recognizers such as Siri or "OK Google".

Commander also has a promising social networking aspect. Users of other MATLAB watches who are also your friends on Facebook or in your Google+ circle can easily highlight snippets of code and share them with you on your watch. Since we had only a few users in our test group, we were not able to fully exploit this feature, but it should become very interesting as the size of the user community grows.

Crown

The crown on the side of the watch replaces a mouse in many ways. Scroll down rows or across columns of a matrix in the variable editor. Scroll through lines of code in the history window. Brush and select data plotted in the figure window.

Magic rank

Here is a simple demonstration executed on the watch.

for n = 3:20
r(n) = rank(magic(n));
end
watch_bar(r)

Presentations

For years I have used MATLAB instead of PowerPoint to give lectures and talks. Now I can do it with my watch, broadcasting to a wireless receiver connected to the projection equipment in the lecture hall. I am scheduled to kick off a High Performance Computing Day at Virginia Tech on Monday, April 6. I hope to debut the MATLAB Watch there.

]]>FunSimulink Easter Eggs
http://feedproxy.google.com/~r/SethOnSimulink/~3/3mRxxvmrQz8/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/seth/2015Q1/AprilFoolsDate.gif"/></div><p>With Easter just a few days away, we thought it was a good time to share some little known Easter eggs that exist in Simulink.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/04/01/simulink-easter-eggs/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4531Wed, 01 Apr 2015 05:01:10 +0000With Easter just a few days away, we thought it was a good time to share some little known Easter eggs that exist in Simulink.

Philosophical Annotations

Many MATLAB users are aware of the interesting response you can get by typing WHY in the command prompt. But, did you know that Simulink Annotations offer similar behavior? Fans of a certain science fiction series might appreciate this one.

Appreciation for Irrational Numbers that are Homophones for Baked Goods

We have a deep appreciation for pi (and for pie). I was recently told that if you string enough of the digits of pi in a series of gain blocks, you could achieve some interesting results. That sounded like a lot of work, so I wrote a script using add_block to construct the model. In the end, the result was definitely interesting.

% Initialize based on current selected block
lastBlock = get_param(gcb,'Name');
blockPos = get_param(gcb,'Position');
for ii = 1:100
% Calculate next digit of pi
lastDigit = rem(floor(pi*power(10,ii-1)),10);
% Generate a new block
blockPos = blockPos+[50 0 50 0];
newBlock = sprintf('Gain%d',ii);
newGain = add_block('built-in/Gain',[gcs,'/',newBlock],'Position',blockPos);
% Populate value and connect
set_param(newGain,'Gain',num2str(lastDigit))
add_line(gcs,[lastBlock,'/1'],[newBlock,'/1'])
lastBlock = newBlock;
end

Little-known Support Package

There are many Support Packages available through the Add-Ons button in the MATLAB toolstrip. Some users have stumbled upon the Advanced Support Package. We can’t reveal the installation steps, but you can see it in action below.

Wait, What Day is it Today?

When you save your Simulink model, the time and date are stored in the Model Properties. Did you also know that on certain dates, an additional comment gets added to the Model History? For example, here’s what you might see if you saved your model today, April 1st.

Now It's Your Turn

What are the best April Fools jokes you've played using Simulink?

]]>FunBlinking text in MATLAB (Experimental Feature!)
http://feedproxy.google.com/~r/mathworks/pick/~3/hlJ_YtEFQ78/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/pick/files/blinky.gif"/></div><p>I'm excited to be back for a special guest visit to the Pick of the Week blog, more than eleven years after I founded it with Doug. I'm extra excited because I have the opportunity to introduce to you a feature I've been playing with that we are hoping to... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/04/01/blinking-text-in-matlab-experimental-feature/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5899Wed, 01 Apr 2015 04:01:18 +0000Doug. I'm extra excited because I have the opportunity to introduce to you a feature I've been playing with that we are hoping to get into MATLAB really soon.
We hear requests for adding features to MATLAB graphics all of the time. Readers of this blog are hopefully well aware now that we released a massive update to MATLAB graphics last year in R2014b. This offered a lot of the enhancements you've been asking for, but missed the one that I hear about the most.
MATLAB has unfortunately lagged a bit behind HTML in providing the ability to build gorgeous grabbing user interfaces (like this one). HTML has long offered the ability make webpages really come to life, with embedded music, animated GIFs, dynamic backgrounds, and best of all - blinking text. I'm slowly chipping away at these, trying to make the case to the MATLAB developers that we need to bring all of these features into MATLAB.
As part of this effort, I've put a prototype of a function on the File Exchange that makes it easy to make any text on a UI blink:

Blink multiple elements on one figure at the same time, all at different speeds

Despite the complete obviousness (to me) that this is the next great graphics feature for MATLAB, I haven't been able to convince too many developers that I am right. Here's where I need your help. If you agree, please post a comment below explaining how important blinking text is to your use of MATLAB. Or, download blink. Lots of times. Every day. When I show that we've got 1,000 downloads a month, I'm sure people will start to listen.
If for some unfathomable reason you think I've got it wrong and that we should be working on something else for graphics and apps, I guess you could say what you think in a comment instead. Just be sure to check out the beautiful App Designer tech preview first to get a sense of what we are already working on. (Note: The App Designer preview only works in MATLAB R2014b.)
Here's to your good health this April, and to many more years of beautiful MATLAB Apps!
One more thing - I would be remiss if I were to forget to wish you a very happy April Fool's Day!
]]>Cache Your Nuts to Eat Later
http://feedproxy.google.com/~r/mathworks/pick/~3/4fZQvn1V0d8/
<p>
Greg's pick this week is cachedcall by Aslak Grinsted.
Contents
Call That Expensive... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/03/27/cache-your-nuts-to-eat-later/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5910Fri, 27 Mar 2015 13:00:09 +0000

Aslak has created an ingenious utility that stores the results to a function call, so that if the same function is called
again with the same inputs, it will return the stored results instead of recomputing the answer.

The definition is to store a resource for future use. In computing, it generally refers to storing data in memory in a fashion that it
can be rapidly accessed at a later time.

In fact developing algorithms for implicitly creating caches is always one of those areas under research in development.

However, in this particular example, the creation of the cache is explicit.

How Is It Done?

The simplified version is, for unique instance of a function call, the results of that function call are stored in a file.
A data table is developed to link the unique function calls to specific data files. When a function call is applied, CACHECALL
determines if that particular function call situation has been performed before. If so, return the stored results, otherwise
call the desired function and store the results for later use.

For me, the part I find most fascinating is Aslak is able to develop unique identifiers for each function call situation based
on the function, inputs, classes, etc. He actually ends up leveraging another File Exchange entry DataHash by Jan Simon.

The DataHash uses MATLAB's inherent ability to call Java code to develop his hash tables and unique identifiers.

Doesn't MATLAB Cache Functions Already?

Yes, but generally not based on specific inputs. When a function is called the first time, if it hasn't been used yet in the
current MATLAB session, then that function must be loaded into memory. Subsequent calls to the same function require less
time because the function is already in memory, but it doesn't store the results of specific computations for future use.

The same thing applies if the function has been edited. It must be reloaded into memory the first time it is called.

Should Every Function Call Be Cached?

Probably not. There is a penalty for the caching process. If the resulting cached results are never used, then that penalty
is wasted. Also caching involves writing to disc, which causes the slowdown in the following results.

It probably only makes sense in cases where there are a limited number of input scenarios that will occur multiple times,
and takes longer than writing to disc.

That's not always predictable, so it may require some profiling to determine if it's worth the cost of caching.

Why Did I Choose This Entry?

This is a clever method to speed up repetitive computations when there are a relatively finite number of possible inputs.

I'm impressed that Aslak was able to figure out how to identify unique function calls. When I saw this entry, it was one
of those "I wish I had thought of that!" moments.

In addition, Aslak does a nice job of architecting the software and documenting his public interface to the caching features.

Could You Use Function Caching?

Have you tried this out? What do you think about it? Should this type of caching be available in MATLAB?

]]>PicksNew online documentation system for R2015a
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/MyubIawWed8/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/steve/2015/repelem-doc-screenshot-phone.png"/></div><p>Earlier this month we shipped R2015a, the first of our two regularly scheduled annual releases. Typically, when there's a new release, I spend some time talking about features that interest me. The feature I want to mention today, though, is a little unusual because it benefits users who haven't even... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2015/03/25/new-online-documentation-system-for-r2015a/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1307Wed, 25 Mar 2015 20:45:33 +0000

Earlier this month we shipped R2015a, the first of our two regularly scheduled annual releases. Typically, when there's a new release, I spend some time talking about features that interest me. The feature I want to mention today, though, is a little unusual because it benefits users who haven't even upgraded to the new release yet.

With R2015a, our Documentation and Documentation Tools teams have overhauled the online documentation on mathworks.com. We've learned a lot from your feedback since the last major documentation overhaul a few years ago, and we are excited about the changes.

The left-side navigation also helps you easily see and get around to all the information that's available on the page you're looking at. Below, the left-side navigation is shown directing you to information about one of the input arguments for repelem.

Now if you haven't upgraded to R2015a yet, and if you've tried to use repelem in your MATLAB, then you've already noticed that repelem isn't there. If you scroll all the way down on repelem reference page, it shows you why:

You can see that repelem is new in R2015a! Customers have long been asking us to provide this information on function reference pages. (Please note that the "introduced in" information is not available for all functions yet. It will take us some time to update every page.)

Finally, the new document system displays reference pages and other content very nicely on mobile devices. Here is a screenshot from my phone:

We'd like to know what you think about the new online doc. That would help our teams as they bring these updates to the documentation system that's included with the product. If you have any feedback you'd like to share, please leave a comment below.

]]>UncategorizedBayesian Brain Teaser
http://feedproxy.google.com/~r/mathworks/loren/~3/2lQyRxA-M5w/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/doors.jpg"/></div><p>Winter in Boston can get quite cold. When we get a lot of snow, we need to take a break after shoveling, and solving puzzles is nice way to spend time indoors. Today's guest blogger, Toshi Takeuchi, gives you an interesting brain teaser, written during one of the many 2015 snowstorms in Boston.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/03/25/bayesian-brain-teaser/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1117Wed, 25 Mar 2015 18:41:31 +0000

Winter in Boston can get quite cold. When we get a lot of snow, we need to take a break after shoveling, and solving puzzles is nice way to spend time indoors. Today's guest blogger, Toshi Takeuchi, gives you an interesting brain teaser, written during one of the many 2015 snowstorms in Boston.

In The Signal and the Noise, the famed data scientist Nate Silver discusses the precarious business of making predictions. The key message of the book is that we should adopt the Bayesian way of thinking - deal with uncertainty probablistically and be flexible to change our prior beliefs as we gather more evidence.

How does it actually work? Let us use a fun brain teaser to work through an example.

The Monty Hall Problem - Simulation

The Monty Hall Problem is loosely based on a game show, "Let's Make a Deal," where the host Monty Hall gives you a choice of three doors, and he tells you that there is a car behind one of them, and you get to keep it if you pick the right door. If you pick a wrong door, you will see a goat, and you get nothing. After you made your choice, Monty opens one of the two doors you didn't pick that reveals a goat. Then he asks, if you would like to stick with your original choice or switch your choice to the other unopened door. At that point the options have been reduced to two doors, and only one of them has the car, so it appears that it is a 50/50 proposition.

Let's test this intuition by running simulations. We will play this game over many trials.

trials = 1000;

First, Monty randomly picks a door from 1, 2, or 3 to hide the car.

car_door = randi(3,trials,1);

You make a guess and a pick a door you think hides the car.

At this point, you can stick with your original choice or switch. My intuition was that it would be a 50/50 proposition. Here is the simulated win rate if you choose to stay.

simulation = zeros(3,1);
simulation(1) = sum((car_door - choice) == 0)/trials*100;
fprintf('Win Rate if stick to your original choice: %.2f%%\n',simulation(1))

Win Rate if stick to your original choice: 33.70%

As you can see, staying with your original choice gives you a much lower win rate than I expected. Let's see what happens if you switch.

Counterintuitively, the simulation shows that it is actually better to switch your choice. Now let's turn to Bayesian inference to figure out why, using Bayes' Rule.

$$P(H|E) = \frac{P(E|H)\cdot P(H)}{P(E)}$$

In this analysis let's assume that:

You initially picked Door 1. This is your hypothesis, H

Monty opened Door 2 to reveal a goat. This is your evidence, E

$$P(H)$$

The car is hidden randomly behind any of the three doors and you picked Door 1. The probability of Door 1 being the car door is 1/3. This is the prior, or P(H).

prior = ones(3,1)*1/3

prior =
0.33333
0.33333
0.33333

$$P(E|H)$$

Now Monty has to choose Door 2 or 3 as the goat door, and he picked Door 2. It's time to think about the conditional probability of the evidence given hypothesis, the likelihood, or P(E|H).

You observed that Monty opened Door 2. His choice is not random - he actually knows the car door, and he needs to avoid it. What does his choice imply?

If the car is behind Door 1, then Monty can choose either Door 2 or 3. So the probability that he chose Door 2 was 1/2.

If the car were behond Door 2, then Monty couldn't have chosen it. So the probability was 0. Getting confused? I will revisit this case below.

If the car is behind Door 3, then Monty had to choose Door 2. So the probability that he chose Door 2 was 1.

Let's go back to case #2. We already know Monty opened Door 2, so why are we talking about the probability of him opening it? Yeah, I know, the probability should be 1 because he actually opened it, right? Well, in this case, the probability is conditioned on the hypothesis under consideration to be true. So in each case you need to figure out what Monty could do based on this assumption. So we know that Door 2 is not the car door, but we still would like to know what Monty might have done if that was the car door.

likelihood = [1/2;0;1]

likelihood =
0.5
0
1

$$P(E|H)\cdot P(H)$$

Now we are ready to update the probabilities of hypotheses based on the new evidence using the equation. Let's calculate the joint probabilities of the prior and likelihood, the numerator of the equation.

joint_prob = likelihood .* prior

joint_prob =
0.16667
0
0.33333

$$P(E)$$

The joint probabilities are not true probabilities because they don't sum up to 1. To get the updated probabilities of hypotheses after the evidence, the posterior, we just need to normalize them by first adding up all the joint probabilities to get the total, and then use that to divide each joint probability. This total is called the normalizing constant, the probability of the evidence under all hypotheses, P(E).

$$P(H|E)$$

The posterior, P(H|E), should be pretty close to the simulation result.

In his book, Nate Silver uses a different formulation to calculate the posterior, but it is actually the same thing. Let's use the example from the book...

"Suppose you are living with a partner and come home from a business trip to discover a strange pair of underwear in your dresser drawer. You will probably ask yourself: what is the probability that your partner is cheating on you?"

x is the prior, the probability of your spouse cheating before you found the underwear, which he puts to 4%

y is the conditional probability or likelihood that you see the underwear because your spouse is cheating, 50%

z is the the conditional probability or likelihood that you see the underwear even though your spouse is not cheating, 5%.

The resulting posterior of your spouse cheating is only 29%, so don't jump to the premature conclusion!

x = 0.04; y = 0.5; z = 0.05;
posterior = x*y/(x*y + z*(1-x))

posterior =
0.29412

The method we use, however, computes the posteriors for both hypotheses - cheating and not cheating.

Bayes' Theorem is used in powerful machine learning algorithms. One fairely straight-forward and useful example is Naive Bayes Classification. It is often used for text classification tasks such as spam filter or sentiment analysis.

Let's do a quick walk-through using a toy example of sentiment analysis. We want to infer the sentiment, positive or negative, of a statement, based on the words contained. Naive Bayes classifier is available as fitcnb function from Statistics Toolbox in MATLAB R2014b.

Turn the training data into a numeric matrix. This is a bag-of-words model that uses frequency of words as predictors.

tokens = unique([training_exmples{:}]);
X = zeros(length(training_exmples),length(tokens));
for i = 1:length(training_exmples)
X(i,:) = ismember(tokens,training_exmples{i});
end

Train a Naive Bayes model with a multinomial distribution, which is the best choice for the bag-of-words example we have. Even though we know that words have grammatical relationships and they are not completely independent, we make a simplifying assumption that our features are independent so we can apply Bayes' Theorem - thus naive. It works just fine none the less.

Test example: "she was a very sweet girl"
Sentiment: positive
Posterior: 0.62

In this toy example, the trained model was able to classify it correctly.

Let's examine the details in the trained model Mdl to see how Bayes' Theorem is used internally.

Here is the prior, and it is 1/2 for each class since we have 3 positive and 3 negative examples. Note that the numbers are listed in order of negative and positive.

What's the actual frequency of positive and negative statements? I have no idea. Often, we don't know the prior, and in such cases we just apply equal probabilities for each hypothesis. This is called a uniform prior.

Here is the likelihood for words 'good' and 'bad'. You see that 'bad' is twice as negative as 'good' and 'good's twice as positive as 'bad'.

good bad
________ ________
negative 0.025 0.05
positive 0.044444 0.022222

Since we are simply interested in relative scores between the negative and positive classification, we can ignore the normalizing constant, which is the denominator of the equation and the same for either class. Therefore all you need is the prior and conditional probabilities of the tokens to calculate relative posterior scores between two hypotheses - postive or negative.

"So what happens if you have 'not' in the statement?" you may ask. A very good point. This simplistic approach doesn't work so well with negation, or irony, for that matter. However, if you train your model on a large number of training examples, it actually works surprisingly well on a majority of real world cases.

To deal with negation like 'not good' or 'not bad' in a more sophisticated way, you need to use word pairs as features instead of single words. For the work you need to put in, however, the performance gain would be incremental.

Likewise, it is easy to deal with irony - just come up with an appropriate way to capture the ironical features! (This is irony).

It's Your Turn

I hope you now see how powerful and useful Bayesian reasoning is and ready to apply to your problems. Do you have an interesting example of Bayesian reasoning? Please share here.

3
% 1 + 3 = 4 -> 2
% 2 + 3 = 5 -> 1
%
% Otherwise he can randomly choose either of the unopened doors to open.
%
% 1 + 1 = 2 REPLACE_WITH_DASH_DASH> 2 or 3
% 2 + 2 = 4 REPLACE_WITH_DASH_DASH> 1 or 3
% 3 + 3 = 6 REPLACE_WITH_DASH_DASH> 1 or 2
%
% This solution is not so elegant, but here it goes.
goat_door = zeros(trials,1);
goat_door(car_door + choice == 3) = 3;
goat_door(car_door + choice == 4 & car_door ~= choice) = 2;
goat_door(car_door + choice == 5) = 1;
goat_door(goat_door == 0 & choice == 1) = randsample([2,3],1);
goat_door(goat_door == 0 & choice == 2) = randsample([1,3],1);
goat_door(goat_door == 0 & choice == 3) = randsample([1,2],1);
%%
% Let's update the preview with the goat door.
T = table(car_door,choice,goat_door);
disp(T(1:5,:))
%%
% At this point, you can stick with your original choice or switch. My
% intuition was that it would be a 50/50 proposition. Here is the simulated
% win rate if you choose to stay.
simulation = zeros(3,1);
simulation(1) = sum((car_door - choice) == 0)/trials*100;
fprintf('Win Rate if stick to your original choice: %.2f%%\n',simulation(1))
%%
% As you can see, staying with your original choice gives you a much lower
% win rate than I expected. Let's see what happens if you switch.
switch_choice = zeros(trials,1);
switch_choice(goat_door + choice == 3) = 3;
switch_choice(goat_door + choice == 4) = 2;
switch_choice(goat_door + choice == 5) = 1;
%%
% Let's update the preview with the switch and compute the win rate.
T = table(car_door,choice,goat_door,switch_choice);
disp(T(1:5,:))
simulation(3) = sum((car_door - switch_choice) == 0)/trials*100;
fprintf('Win Rate if switch your choice: %.2f%%\n',simulation(3))
%% Why Switching Is Better - Bayesian Analysis
% Counterintuitively, the simulation shows that it is actually better to
% switch your choice. Now let's turn to
% to
% figure out why, using Bayes' Rule.
%
% $$P(H|E) = \frac{P(E|H)\cdot P(H)}{P(E)}$$
%
% In this analysis let's assume that:
%
% * You initially picked Door 1. This is your hypothesis, |*H*|
% * Monty opened Door 2 to reveal a goat. This is your evidence, |*E*|
%
%%
% $$P(H)$$
%
% The car is hidden randomly behind any of the three doors and you picked
% Door 1. The probability of Door 1 being the car door is 1/3. This is the
% *prior*, or |P(H)|.
prior = ones(3,1)*1/3
%%
% $$P(E|H)$$
%
% Now Monty has to choose Door 2 or 3 as the goat door, and he picked Door
% 2. It's time to think about the conditional probability of the evidence
% given hypothesis, the *likelihood*, or |P(E|H)|.
%
% You observed that Monty opened Door 2. His choice is not
% random - he actually knows the car door, and he needs to avoid it. What
% does his choice imply?
%
% # If the car is behind Door 1, then Monty can choose either Door 2 or 3.
% So the probability that he chose Door 2 was 1/2.
% # If the car were behond Door 2, then Monty couldn't have chosen it. So
% the probability was 0. Getting confused? I will revisit this case below.
% # If the car is behind Door 3, then Monty had to choose Door 2. So the
% probability that he chose Door 2 was 1.
%
% Let's go back to case #2. We already know Monty opened Door 2, so why are
% we talking about the probability of him opening it? Yeah, I know, the
% probability should be 1 because he actually opened it, right? Well, in
% this case, the probability is conditioned on the hypothesis under
% consideration to be true. So in each case you need to figure out what
% Monty could do based on this assumption. So we know that Door 2 is not
% the car door, but we still would like to know what Monty might have done
% if that was the car door.
likelihood = [1/2;0;1]
%%
% $$P(E|H)\cdot P(H)$$
%
% Now we are ready to update the probabilities of hypotheses based on the
% new evidence using the equation. Let's calculate the joint probabilities
% of the prior and likelihood, the numerator of the equation.
%
joint_prob = likelihood .* prior
%%
% $$P(E)$$
%
% The joint probabilities are not true probabilities because they don't sum
% up to 1. To get the updated probabilities of hypotheses after the
% evidence, the *posterior*, we just need to normalize them by first adding
% up all the joint probabilities to get the total, and then use that to
% divide each joint probability. This total is called the
% ,
% the probability of the evidence under all hypotheses, |P(E)|.
%
% $$P(H|E)$$
%
% The posterior, |P(H|E)|, should be pretty close to the simulation result.
posterior = joint_prob/sum(joint_prob);
compare = table(posterior,simulation,'RowNames',{'Stay','N/A','Switch'});
disp(compare([1 3],:))
%% Nate Silver's Method
% In his book, Nate Silver uses a different formulation to calculate the
% posterior, but it is actually the same thing. Let's use the example from
% the book...
%
% "Suppose you are living with a partner and come home from a
% business trip to discover a strange pair of underwear in your dresser
% drawer. You will probably ask yourself: what is the probability that your
% partner is cheating on you?"
%
% * |x| is the prior, the probability of your spouse cheating before you
% found the underwear, which he puts to 4%
% * |y| is the conditional probability or likelihood that you see the
% underwear because your spouse is cheating, 50%
% * |z| is the the conditional probability or likelihood that you see the
% underwear even though your spouse is not cheating, 5%.
% * The resulting posterior of your spouse cheating is only 29%, so don't
% jump to the premature conclusion!
x = 0.04; y = 0.5; z = 0.05;
posterior = x*y/(x*y + z*(1-x))
%%
% The method we use, however, computes the posteriors for both hypotheses -
% cheating and not cheating.
prior = [x; 1-x]; likelihood = [y; z];
joint_prob = likelihood .* prior;
posterior = joint_prob ./ sum(joint_prob)
%% Bayesian in Machine Learning
% Bayes' Theorem is used in powerful machine learning algorithms. One
% fairely straight-forward and useful example is
%
. It is often used for text classification
% tasks such as spam filter or sentiment analysis.
%
% Let's do a quick walk-through using a toy example of sentiment analysis.
% We want to infer the sentiment, positive or negative, of a statement,
% based on the words contained. Naive Bayes classifier is available as
% function from
% Statistics Toolbox in MATLAB R2014b.
%
% Let's load the training examples.
training_exmples = {{'i','love','this','food','very','much'}; ...
{'that','left','a','bad','taste'}; ...
{'the','movie','was','boring'}; ...
{'i','had','a','very','good','time'};
{'she','has','such','a','nice','sweet','dog'}; ...
{'the','ending','was','very','sad'}};
%%
% Those examples are pre-labeled as positive or negative.
sentiment_class = {'positive';'negative';'negative'; ...
'positive';'positive';'negative'};
%%
% Here is the unlabeled example we will use to test the classifier.
test_example = {'she','was','a','very','sweet','girl'};
%%
% Turn the training data into a numeric matrix. This is a bag-of-words
% model that uses frequency of words as predictors.
tokens = unique([training_exmples{:}]);
X = zeros(length(training_exmples),length(tokens));
for i = 1:length(training_exmples)
X(i,:) = ismember(tokens,training_exmples{i});
end
%%
% Train a Naive Bayes model with a multinomial distribution, which is the
% best choice for the bag-of-words example we have. Even though we know
% that words have grammatical relationships and they are not completely
% independent, we make a simplifying assumption that our features are
% independent so we can apply Bayes' Theorem - thus _naive_. It works
% just fine none the less.
Mdl = fitcnb(X,sentiment_class,'Distribution','mn','PredictorNames',tokens);
%%
% Now check the trained model with the test example.
[label,post,~]= predict(Mdl,double(ismember(Mdl.PredictorNames,test_example)));
fprintf('Test example: "%s"\nSentiment: %s\nPosterior: %.2f\n',...
strjoin(test_example,' '),label{1},post(strcmp(Mdl.ClassNames,label)))
%%
% In this toy example, the trained model was able to classify it correctly.
%
% Let's examine the details in the trained model |Mdl| to see how Bayes'
% Theorem is used internally.
%
% Here is the
% prior, and it is 1/2 for each class since we have 3 positive and 3
% negative examples. Note that the numbers are listed in order of
% negative and positive.
disp(table(Mdl.Prior(1),Mdl.Prior(2),'VariableNames',Mdl.ClassNames))
%%
% What's the actual frequency of positive and negative statements? I have
% no idea. Often, we don't know the prior, and in such cases we just apply
% equal probabilities for each hypothesis. This is called a uniform prior.
%
% Here is the likelihood for words 'good' and 'bad'. You see that 'bad' is
% twice as negative as 'good' and 'good's twice as positive as 'bad'.
disp(table(cell2mat(Mdl.DistributionParameters(:,ismember(Mdl.PredictorNames,'good'))),...
cell2mat(Mdl.DistributionParameters(:,ismember(Mdl.PredictorNames,'bad'))),...
'VariableNames',{'good','bad'},'RowNames',Mdl.ClassNames))
%%
% Since we are simply interested in relative scores between the negative
% and positive classification, we can ignore the normalizing constant,
% which is the denominator of the equation and the same for either class.
% Therefore all you need is the prior and conditional probabilities of the
% tokens to calculate relative posterior scores between two hypotheses -
% postive or negative.
%
% "So what happens if you have 'not' in the statement?" you may ask. A very
% good point. This simplistic approach doesn't work so well with negation,
% or irony, for that matter. However, if you train your model on a large
% number of training examples, it actually works surprisingly well on a
% majority of real world cases.
%
% To deal with negation like 'not good' or 'not bad' in a more
% sophisticated way, you need to use word pairs as features instead of
% single words. For the work you need to put in, however, the performance
% gain would be incremental.
%
% Likewise, it is easy to deal with irony - just come up with an
% appropriate way to capture the ironical features! (This is irony).
%% It's Your Turn
% I hope you now see how powerful and useful Bayesian reasoning is and
% ready to apply to your problems. Do you have an interesting example of
% Bayesian reasoning? Please share
% .
##### SOURCE END ##### 6db174c34e5a4e19b4971377057aebfd
-->]]>MATLAB Virtual Conference 2015
http://feedproxy.google.com/~r/SethOnSimulink/~3/Yx-cKNqwMmI/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/seth/2015Q1/VirtualConference.png"/></div><p>In case your were not aware, this March 25th is the MATLAB Virtual Conference.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/03/24/matlab-virtual-conference-2015/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4513Wed, 25 Mar 2015 00:56:39 +0000In case your were not aware, this March 25th is the MATLAB Virtual Conference.

Whether you are a beginner or an expert with MATLAB and Simulink, I am sure you can find interesting sessions for you. Here are the sessions I think will interest most Simulink users.

Introduction to Simulink: Quadcopter Simulation and Control

In this talk, Ryan Gordon will go through the steps he went through to design a controller for a Parrot® AR.Drone® quadcopter:

Modeling the quadrotor using Simulink and SimMechanics

Importing CAD files for visualization

Designing and tuning a controller using Simulink Control Design

Deploying the controller on the ARM Cortex processor of the Parrot® AR.Drone®

I like Arduino boards, they offer many possibilities... especially when programmed with MATLAB and Simulink.

If you are interested in getting started with Arduino, don't miss this session. Dan Seal will go through a complete example that will give you all the basic knowledge you need to get started programming your Arduino boards with Simulink.

What’s New in MATLAB and Simulink R2015a

If you are using Simulink every day, you will benefit from learning about new features you can integrate into your work.

I plan to blog about this new toolbox soon, but if you can't wait, this talk by Yanliang Zhang is for you. You will learn how the Robot Operating System (ROS) interface allows you to deploy your Simulink models on any ROS-enabled device.

Presentation of the Platform for the Development of Aircraft Engine Monitoring Algorithms: SAMANTA

In this presentation, Jérôme Lacaille and Aurélie Gouby, from Snecma, will talk about an environment the created to design and develop health monitoring algorithms for aircraft engines.

I don't know much about this topic, but I will attend this session to learn more about this very interesting project.

]]>Expand
http://feedproxy.google.com/~r/mathworks/pick/~3/S9yc23g1gk4/
<p>
Sean's pick this week is Expand by Matt Fig.
Contents
Background
... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/03/20/expand/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5895Fri, 20 Mar 2015 13:00:38 +0000

Matt Fig is one of MATLAB Central's top contributors. Even though he's been on hiatus for over two years, he's still the
10th highest contributor on Answers and the 14th most downloaded author on the File Exchange.

Expand

My pick today is a contribution of Matt's that made an appearance somewhere in the code base for my Master's thesis. I was
working with large 3-dimensional CT images and computation speed was important. Here's what it does:

As you can see x has been expanded thrice along rows and twice along columns ([3 2]). Expand also works along N-dimensional arrays but does
have a limitation that the expansion must match the number of dimensions. This can be circumvented with a second call to
repmat to expand along higher dimensions.

Kron

The 2-dimensional expand operation has always been possible with kron, but it requires creating a ones matrix and, at least back in my grad school era, was slower than expand. Since I was working with 3-dimensional images, this was not an option.

Starting in R2015a, which shipped the first week of March, there is a new function that does the same thing as expand and
more built in to MATLAB: meet repelem.

repelem on a vector can also vary the number of repeats, very nice for building up a vector or run length encoding.

v123 = repelem([2 4 7],1:3)

v123 =
2 4 4 7 7 7

R2015a Other

Browsing through the release notes for R2015a and having had the opportunity to beat on the R2015a Prerelease, here are a few of my other favorite new features:

discretize - similar to histcounts or histc but skips the histogram counting part.

Array Size Limitation - Helps accidentally creating large arrays that would cause you to run out of RAM.

Tab Complete for Object Authoring - Finally!

ismembertol and uniquetol - ismember and unique with tolerance; these will save me a few calls to discretize.

Comments

Thanks Matt for providing the expand functionality; it has been very useful for nearly 6 years!

Give expand, repelem, and R2015a a try and let us know what you think about any of them here or leave a comment for Matt.

]]>Getting the most out of Rapid Accelerator mode
http://feedproxy.google.com/~r/SethOnSimulink/~3/xcY7HsuH2Rg/
<div class="overview-image"><img src="http://blogs.mathworks.com/seth/files/feature_image/rapidAccelerator.png" class="img-responsive attachment-post-thumbnail wp-post-image" alt="rapidAccelerator"/></div><p>When working with Simulink, it's important to configure your model and workflows to be as efficient as possible for the task you are working on.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/03/19/getting-the-most-out-of-rapid-accelerator-mode/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4486Thu, 19 Mar 2015 17:32:25 +0000When working with Simulink, it's important to configure your model and workflows to be as efficient as possible for the task you are working on.

Today I will show how to setup a model to get the maximum performance when you need to run many simulations and vary parameters. This is often the case if you do monte-carlo simulation, or system optimization.

Rapid Accelerator

As explained here, the rapid accelerator mode can speed up simulation by generating an executable for your model. The speedup can vary between models for various reasons. Personally, I have seen some models run more than 10 times faster in rapid accelerator mode.

RapidAcceleratorUpToDateCheck off

When you click the play button in a model or use the sim command, Simulink verifies if your model has changed or not. If the model has not changed structurally, the rapid accelerator executable is not re-generated.

One nice thing with rapid accelerator is that, once the executable has been generated, it is possible to skip this part where Simulink verifies if the model has changed or not. For large models, this can save a lot of time, bringing the initialization time to zero.

Let's see how this can be done using a simple demo like sldemo_bounce where I want to vary the coefficient of restitution, which I specified as k in the block dialog.

If I want to be able to tune the value of k, I need to define it as a Simulink.Parameter object:

k = Simulink.Parameter
k.CoderInfo.StorageClass = 'SimulinkGlobal';
k.Value = -0.9;

and I enable the Inline Parameters option in the model configuration.

Once this is done, I explicitly build the rapid accelerator executable:

As you can see, this returns a structure containing all the information relevant to tunable parameters in the model. Using this structure, we can create an array of run-time parameter structures for all the values we want to run:

k_values = [-0.9:0.1:-0.1];
for i = 1:length(k_values)
paramSet(i) = Simulink.BlockDiagram.modifyTunableParameters(rtp, 'k', k_values(i));
end

And we are ready to simulate, by passing to the sim command, we pass the parameters RapidAcceleratorUpToDateCheck and RapidAcceleratorParameterSets.

for i = 1:length(k_values)
simout(i) = sim(mdl,'SimulationMode','rapid',...'RapidAcceleratorUpToDateCheck','off', ...'RapidAcceleratorParameterSets',paramSet(i));
end

Now it's your turn

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

]]>MATLAB, Landsat 8, and AWS Public Data Sets
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/FWdGHItNqQQ/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/steve/files/AWS-image-300x168.jpg"/></div><p>A few weeks ago, a fellow developer pulled me into his office to show me something he had been working on. It was very cool! Even better, it's now available for you to try. I'd like to introduce guest blogger Bruce Tannenbaum, product marketing manager for image processing and test... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2015/03/19/matlab-landsat-8-aws/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1294Thu, 19 Mar 2015 15:00:37 +0000A few weeks ago, a fellow developer pulled me into his office to show me something he had been working on. It was very cool! Even better, it's now available for you to try. I'd like to introduce guest blogger Bruce Tannenbaum, product marketing manager for image processing and test & measurement applications, to tell you all about it. Thanks, Bruce!

MathWorks is excited to announce a freely-downloadable MATLAB based tool for accessing, processing, and visualizing Landsat 8 data hosted by AWS as part of its Landsat on AWS Public Data Sets. With this tool, you can create a map display of scene locations with markers that show each scene’s metadata. You can download and visualize thumbnails and any of the 11 spectral bands. The tool includes clickable links for automatically combining and processing spectral bands in a variety of approaches, such as NDVI and color infrared. You can then visualize processed results in MATLAB and create map displays that help convey context about where the data is located on the Earth. This interactive tool is available on MATLAB Central. Watch the video below to see it in action.

How MATLAB and Landsat 8 Work Together

MATLAB provides a great environment for working with Earth observation data. It supports a wide range of file formats, including the GeoTIFF format used for Landsat 8 imagery. It supports accessing Earth science data over the Internet, such as the http requests used by AWS Public Data Sets. The Image Processing Toolbox extends MATLAB with a comprehensive set of image processing algorithms and tools for image processing, analysis, visualization, and algorithm development. The interactive tool for Landsat 8 imagery takes advantage of processing techniques in the toolbox, such as adaptive histogram equalization.

One challenge of working with Landsat 8 imagery is that you need to know where each scene is located. A scene represents a specific time and a rectangular surface region of about 185 by 180 kilometers. To make it easier to locate scenes, our tool creates a map display with markers showing the centroid of each scene within the field of view. Additionally, the tool color-codes the markers to show cloud coverage, which is encapsulated in the metadata file for each scene. The most useful images have low cloud coverage; the tool makes it easy to find images with 0-10% or 10-20% cloud coverage, which are the areas that are clearest.

In summary, Landsat 8 is an incredible resource for global change research and has been used in a diverse array of scientific endeavors including the monitoring of deforestation, population growth, and glacier recession. The tool offers a great way for MATLAB users to build on the foundation of AWS support for Landsat 8 imagery. It can also be run on an EC2 instance, which avoids the download time and allows you to process many images in a shorter amount of time. MathWorks would like to hear your feedback, including any additional features you would like to have.

]]>UncategorizedThe Ultimate Pi Day is Over… Now What?
http://feedproxy.google.com/~r/mathworks/loren/~3/Z6LxojvlyMQ/
<p>Today's post is from guest blogger Dan Seal who works in the MATLAB Product Marketing team here at MathWorks. Dan recently celebrated Pi Day on 3/14/15 and wanted to know what other holidays he might be able to celebrate in the future.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/03/16/the-ultimate-pi-day-is-over-now-what/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1150Mon, 16 Mar 2015 18:28:09 +0000

Today's post is from guest blogger Dan Seal who works in the MATLAB Product Marketing team here at MathWorks. Dan recently celebrated Pi Day on 3/14/15 and wanted to know what other holidays he might be able to celebrate in the future.

Every year, mathematicians and their geeky cohort celebrate Pi Day on the fourteenth of March, a date whose written representation (3/14) looks like the first three digits of the mathematical constant pi (3.14). But this year’s Pi Day was extra special and generated additional buzz because when you include the year in the written date of 3/14/15, it gives us two more digits of everyone’s favorite irrational number. And if you were like me, you may have even marked the minute when it was 3/14/15 9:26, both in the AM and again in the PM.

But the opportunity to observe these additional digits of Pi on the calendar is something that comes around only once per century, so I don’t expect to see it again. This got me wondering if there are any other special dates I might be able to look forward to. M/DD/YY is not the only representation of a date, and Pi is not the only special constant out there with multiple digits. So I collected all my favorite mathematical and physical constants along with a slew of possible date representations, and threw them all into MATLAB to see what it would give me.

Starting small

I decided to start out with just two constants which I could use to test my work as I went along. I selected Pi and the square root of 2. Fortunately, MATLAB can calculate a very accurate approximation of these values, so I didn’t have to look up the digits.

names = {
'Archimedes'' Constant, or Pi''Pythagoras'' Constant, or Root 2'
};
numbers = [
pi
sqrt(2)
];

Defining the possible date formats

Next, I wanted to construct all the allowable date formats. I wanted to consider only dates that included the month, day, and year. However, I figure it could be a 2- or 4-digit year. Similarly, I allowed the month and day to each be 1 or 2 digits. And finally, I allowed the month, day, and year to be ordered in any of four different arrangements:

the American version of month/day/year,

the common standard outside the US, day/month/year,

the common year-first version, useful for sorting, year/month/day,

and for good measure, the alternate year-first arrangement, year/day/month.

MATLAB uses the datetime data type, introduced in R2014b to easily manipulate and format dates and times. Based on my formatting requirements above, I was able to build all the possible datetime format strings that would meet my criteria.

yearFormats = {'yy''yyyy'}; % 2 or 4 digit year
monthFormats = {'M''MM' }; % 1 or 2 digit month
dayFormats = {'d''dd' }; % 1 or 2 digit day
availableFormats{length(yearFormats)*length(monthFormats)*length(dayFormats)*4} = '';
c = 1;
for y = yearFormats
for m = monthFormats
for d = dayFormats
availableFormats{c } = [d{:} '/' m{:} '/' y{:}];
availableFormats{c+1} = [m{:} '/' d{:} '/' y{:}];
availableFormats{c+2} = [y{:} '/' m{:} '/' d{:}];
availableFormats{c+3} = [y{:} '/' d{:} '/' m{:}];
c = c + 4;
endendend
disp(availableFormats')

The next task was to extract the digits from my numbers. I figured I’d be using some quantities with large positive or negative exponents, so my technique had to be robust to different types of numbers. I decided the best approach was to convert the numbers to strings using a consistent format so I could easily remove any decimal points and exponential terms. I found that sprintf with an exponential format could achieve this. Keeping nine digits after the decimal point ensured that I had enough precision to build even my longest date (which required 8 digits). The char and arrayfun commands allowed me to run sprintf on each number individually and then build it into a character array.

Next, I removed the decimal point and keep only the first eight digits, which is the most I will need.

digitsStr = digitsStr(:,[1 3:9]);
disp(digitsStr)

31415926
14142135

Building the list of special dates

Finally, I built a loop to test each of my mathematical constants against all the possible date formats and see which ones could be constructed as real dates. I tried out each possible date format using a try/catch statement. If the date I tried to construct wasn’t a real date, I would catch that particular error. But if it was, it would be added to the list of special holidays to celebrate and the name of that constant would be recorded as well. I also specified a pivot year to define which century two-digit years belong to. Since I’m looking or future holidays, I decided that all two-digit years should represent dates sometime after 2015.

% Initialize arrays
maxPossibleDates = length(numbers)*length(availableFormats);
holidayDates(maxPossibleDates) = datetime;
holidayDates.Format = 'MMM dd, yyyy';
holidayFormats{maxPossibleDates} = '';
holidayNames{maxPossibleDates} = '';
% Initialize counter
c = 1;
for ii = 1:length(numbers) %Consider each numberfor jj = 1:length(availableFormats) %For each number, consider each format
fmt = availableFormats{jj};
testFmt = strrep(fmt,'/','');
strLength = length(testFmt);
testStr = digitsStr(ii,1:strLength);
try%If format is valid, add it to the list
holidayDates(c) = datetime(testStr,'InputFormat',testFmt,'PivotYear',2015);
holidayFormats{c} = fmt;
holidayNames{c} = names{ii};
c = c + 1;
catch err %If format errors, make sure it is the expected errorif ~strcmp(err.identifier,'MATLAB:datetime:ParseErr')
rethrow(err)
endendendend% Trim empty values from end of lists
holidayDates = holidayDates(1:c-1);
holidayFormats = holidayFormats(1:c-1);
holidayNames = holidayNames(1:c-1);

Because my dates are represented as datetimes, I can easily put them in chronological order with the sort function.

Let’s see how many opportunities there are to celebrate Pi Day or Root 2 Day.

for ii = 1:c-1
specialDate = holidayDates(ii);
specialDate.Format = holidayFormats{ii};
fprintf('On %s (%s), we can celebrate %s Day\n',...
char(holidayDates(ii)),char(specialDate),holidayNames{ii})
end

It looks like there are a few Pi Days we can celebrate over the next 30 years. And then there are some Root 2 Days to celebrate after that.

On Jan 02, 1414 (1414/2/1), we can celebrate Pythagoras' Constant, or Root 2 Day
On Feb 01, 1414 (1414/2/1), we can celebrate Pythagoras' Constant, or Root 2 Day
On Feb 13, 1414 (1414/2/13), we can celebrate Pythagoras' Constant, or Root 2 Day
On Mar 21, 1414 (1414/21/3), we can celebrate Pythagoras' Constant, or Root 2 Day
On Jan 04, 1421 (1/4/1421), we can celebrate Pythagoras' Constant, or Root 2 Day
On Apr 01, 1421 (1/4/1421), we can celebrate Pythagoras' Constant, or Root 2 Day
On Mar 14, 1592 (3/14/1592), we can celebrate Archimedes' Constant, or Pi Day
On Mar 14, 2015 (3/14/15), we can celebrate Archimedes' Constant, or Pi Day
On Jan 04, 2031 (31/4/1), we can celebrate Archimedes' Constant, or Pi Day
On Apr 01, 2031 (31/4/1), we can celebrate Archimedes' Constant, or Pi Day
On Apr 15, 2031 (31/4/15), we can celebrate Archimedes' Constant, or Pi Day
On Jan 03, 2041 (3/1/41), we can celebrate Archimedes' Constant, or Pi Day
On Mar 01, 2041 (3/1/41), we can celebrate Archimedes' Constant, or Pi Day
On Jan 14, 2042 (14/1/42), we can celebrate Pythagoras' Constant, or Root 2 Day
On Jan 04, 2114 (1/4/14), we can celebrate Pythagoras' Constant, or Root 2 Day
On Jan 04, 2114 (14/1/4), we can celebrate Pythagoras' Constant, or Root 2 Day
On Feb 14, 2114 (14/14/2), we can celebrate Pythagoras' Constant, or Root 2 Day
On Apr 01, 2114 (1/4/14), we can celebrate Pythagoras' Constant, or Root 2 Day
On Apr 01, 2114 (14/1/4), we can celebrate Pythagoras' Constant, or Root 2 Day
On May 09, 3141 (3141/5/9), we can celebrate Archimedes' Constant, or Pi Day
On Sep 05, 3141 (3141/5/9), we can celebrate Archimedes' Constant, or Pi Day
On Jan 03, 4159 (3/1/4159), we can celebrate Archimedes' Constant, or Pi Day
On Mar 01, 4159 (3/1/4159), we can celebrate Archimedes' Constant, or Pi Day
On Jan 14, 4213 (14/1/4213), we can celebrate Pythagoras' Constant, or Root 2 Day

Generalizing to a larger set of inputs

Since this gave me good results, I could put the entire analysis into a function and then call it with a much longer list of my favorite mathematical and physical constants. What are the next holidays coming up? Be sure to mark your calendar.

names = {
'Archimedes'' Constant, or Pi''Euler''s Number, or e''Pythagoras'' Constant, or Root 2''Golden Ratio, or Phi''Root 3''Root 5''Speed of Light, or c''Avogadro''s Number, or N_A''Standard Gravity, or g''Gravitational Constant, or G''Planck Constant, or h''Ideal Gas Constant, or R''Coulomb''s Constant, or k_e''Elementary Charge, or e''Bohr Radius, or a_0''Electron Mass, or m_e''Proton Mass, or m_p''Atomic Mass Constant, or m_u''Boltzmann Constant, or k_B''Faraday Constant, or F''Stefan-Boltzmann Constant, or Theta''Standard Atmosphere, or atm''Planck Length, or l_P''Planck Mass, or m_P''Planck Time, or t_P''Planck Charge, or q_P''Planck Temperature, or T_P'
};
numbers = [
pi
exp(1)
sqrt(2)
(1+sqrt(5))/2
sqrt(3)
sqrt(5)
299792458
6.02214129e23
9.80665
6.67384e-11
6.62606957e-34
8.3144621
8.9875517873681764e9
1.602176565e-19
5.2917721092e-11
9.10938291e-31
1.672621777e-27
1.660538921e-27
1.3806488e-23
96485.3365
5.670373e-8
101325
1.616199e-35
2.17651e-8
5.39106e-44
1.875545956e-18
1.416833e32
];
holidayForFamousConstants(numbers,names)

On Jan 06, 0217 (1/6/0217), we can celebrate Elementary Charge, or e Day
On Jun 01, 0217 (1/6/0217), we can celebrate Elementary Charge, or e Day
On Jun 16, 0538 (16/6/0538), we can celebrate Atomic Mass Constant, or m_u Day
On Aug 13, 0648 (13/8/0648), we can celebrate Boltzmann Constant, or k_B Day
On Aug 09, 0665 (9/8/0665), we can celebrate Standard Gravity, or g Day
On Sep 08, 0665 (9/8/0665), we can celebrate Standard Gravity, or g Day
On Jan 09, 0938 (9/1/0938), we can celebrate Electron Mass, or m_e Day
On Sep 01, 0938 (9/1/0938), we can celebrate Electron Mass, or m_e Day
On Feb 05, 1013 (1013/2/5), we can celebrate Standard Atmosphere, or atm Day
On May 02, 1013 (1013/2/5), we can celebrate Standard Atmosphere, or atm Day
On Apr 06, 1380 (1380/6/4), we can celebrate Boltzmann Constant, or k_B Day
On Jun 04, 1380 (1380/6/4), we can celebrate Boltzmann Constant, or k_B Day
On Jan 02, 1414 (1414/2/1), we can celebrate Pythagoras' Constant, or Root 2 Day
On Feb 01, 1414 (1414/2/1), we can celebrate Pythagoras' Constant, or Root 2 Day
On Feb 13, 1414 (1414/2/13), we can celebrate Pythagoras' Constant, or Root 2 Day
On Mar 21, 1414 (1414/21/3), we can celebrate Pythagoras' Constant, or Root 2 Day
On Mar 08, 1416 (1416/8/3), we can celebrate Planck Temperature, or T_P Day
On Aug 03, 1416 (1416/8/3), we can celebrate Planck Temperature, or T_P Day
On Jan 04, 1421 (1/4/1421), we can celebrate Pythagoras' Constant, or Root 2 Day
On Apr 01, 1421 (1/4/1421), we can celebrate Pythagoras' Constant, or Root 2 Day
On Mar 08, 1446 (8/3/1446), we can celebrate Ideal Gas Constant, or R Day
On Aug 03, 1446 (8/3/1446), we can celebrate Ideal Gas Constant, or R Day
On Mar 14, 1592 (3/14/1592), we can celebrate Archimedes' Constant, or Pi Day
On Jan 07, 1602 (1602/1/7), we can celebrate Elementary Charge, or e Day
On Jun 17, 1602 (1602/17/6), we can celebrate Elementary Charge, or e Day
On Jul 01, 1602 (1602/1/7), we can celebrate Elementary Charge, or e Day
On Jan 09, 1616 (1616/1/9), we can celebrate Planck Length, or l_P Day
On Sep 01, 1616 (1616/1/9), we can celebrate Planck Length, or l_P Day
On Sep 19, 1616 (1616/19/9), we can celebrate Planck Length, or l_P Day
On Mar 03, 1618 (1618/03/3), we can celebrate Golden Ratio, or Phi Day
On Mar 03, 1618 (1618/03/3), we can celebrate Golden Ratio, or Phi Day
On Jan 06, 1619 (1/6/1619), we can celebrate Planck Length, or l_P Day
On Jun 01, 1619 (1/6/1619), we can celebrate Planck Length, or l_P Day
On Mar 05, 1660 (1660/5/3), we can celebrate Atomic Mass Constant, or m_u Day
On May 03, 1660 (1660/5/3), we can celebrate Atomic Mass Constant, or m_u Day
On Feb 06, 1672 (1672/6/2), we can celebrate Proton Mass, or m_p Day
On Jun 02, 1672 (1672/6/2), we can celebrate Proton Mass, or m_p Day
On Jun 21, 1672 (1672/6/21), we can celebrate Proton Mass, or m_p Day
On Jan 04, 1683 (1/4/1683), we can celebrate Planck Temperature, or T_P Day
On Apr 01, 1683 (1/4/1683), we can celebrate Planck Temperature, or T_P Day
On May 08, 1732 (1732/05/08), we can celebrate Root 3 Day
On Aug 05, 1732 (1732/05/08), we can celebrate Root 3 Day
On Feb 16, 1765 (16/02/1765), we can celebrate Elementary Charge, or e Day
On May 29, 1772 (5/29/1772), we can celebrate Bohr Radius, or a_0 Day
On Jan 06, 1803 (1/6/1803), we can celebrate Golden Ratio, or Phi Day
On Jun 01, 1803 (1/6/1803), we can celebrate Golden Ratio, or Phi Day
On Feb 07, 1828 (2/7/1828), we can celebrate Euler's Number, or e Day
On Jul 02, 1828 (2/7/1828), we can celebrate Euler's Number, or e Day
On Apr 05, 1875 (1875/5/4), we can celebrate Planck Charge, or q_P Day
On May 04, 1875 (1875/5/4), we can celebrate Planck Charge, or q_P Day
On Mar 14, 2015 (3/14/15), we can celebrate Archimedes' Constant, or Pi Day
On Jan 02, 2016 (16/02/1), we can celebrate Elementary Charge, or e Day
On Jan 04, 2016 (1/4/16), we can celebrate Planck Temperature, or T_P Day
On Jan 06, 2016 (1/6/16), we can celebrate Planck Length, or l_P Day
On Jan 06, 2016 (16/1/6), we can celebrate Planck Length, or l_P Day
On Jan 08, 2016 (16/1/8), we can celebrate Golden Ratio, or Phi Day
On Jan 16, 2016 (16/16/1), we can celebrate Planck Length, or l_P Day
On Feb 01, 2016 (16/02/1), we can celebrate Elementary Charge, or e Day
On Feb 07, 2016 (16/7/2), we can celebrate Proton Mass, or m_p Day
On Feb 17, 2016 (16/02/17), we can celebrate Elementary Charge, or e Day
On Mar 18, 2016 (16/18/03), we can celebrate Golden Ratio, or Phi Day
On Apr 01, 2016 (1/4/16), we can celebrate Planck Temperature, or T_P Day
On May 06, 2016 (16/6/05), we can celebrate Atomic Mass Constant, or m_u Day
On Jun 01, 2016 (1/6/16), we can celebrate Planck Length, or l_P Day
On Jun 01, 2016 (16/1/6), we can celebrate Planck Length, or l_P Day
On Jun 05, 2016 (16/6/05), we can celebrate Atomic Mass Constant, or m_u Day
On Jul 02, 2016 (16/7/2), we can celebrate Proton Mass, or m_p Day
On Jul 26, 2016 (16/7/26), we can celebrate Proton Mass, or m_p Day
On Aug 01, 2016 (16/1/8), we can celebrate Golden Ratio, or Phi Day
On Feb 03, 2017 (17/3/2), we can celebrate Root 3 Day
On Feb 16, 2017 (16/02/17), we can celebrate Elementary Charge, or e Day
On Mar 02, 2017 (17/3/2), we can celebrate Root 3 Day
On Mar 20, 2017 (17/3/20), we can celebrate Root 3 Day
On May 29, 2017 (5/29/17), we can celebrate Bohr Radius, or a_0 Day
On Jan 06, 2018 (1/6/18), we can celebrate Golden Ratio, or Phi Day
On Feb 07, 2018 (2/7/18), we can celebrate Euler's Number, or e Day
On May 07, 2018 (18/7/5), we can celebrate Planck Charge, or q_P Day
On Jun 01, 2018 (1/6/18), we can celebrate Golden Ratio, or Phi Day
On Jul 02, 2018 (2/7/18), we can celebrate Euler's Number, or e Day
On Jul 05, 2018 (18/7/5), we can celebrate Planck Charge, or q_P Day
On Mar 17, 2020 (17/3/20), we can celebrate Root 3 Day
On Feb 06, 2021 (6/02/21), we can celebrate Avogadro's Number, or N_A Day
On Jun 02, 2021 (6/02/21), we can celebrate Avogadro's Number, or N_A Day
On Jun 07, 2021 (21/7/6), we can celebrate Planck Mass, or m_P Day
On Jul 06, 2021 (21/7/6), we can celebrate Planck Mass, or m_P Day
On Mar 06, 2022 (22/3/6), we can celebrate Root 5 Day
On Jun 03, 2022 (22/3/6), we can celebrate Root 5 Day
On Oct 13, 2025 (10/13/25), we can celebrate Standard Atmosphere, or atm Day
On Jun 06, 2026 (6/6/26), we can celebrate Planck Constant, or h Day
On Jun 06, 2026 (6/6/26), we can celebrate Planck Constant, or h Day
On Jul 16, 2026 (16/7/26), we can celebrate Proton Mass, or m_p Day
On Jan 08, 2027 (27/1/8), we can celebrate Euler's Number, or e Day
On Feb 18, 2027 (27/18/2), we can celebrate Euler's Number, or e Day
On Aug 01, 2027 (27/1/8), we can celebrate Euler's Number, or e Day
On Jul 09, 2029 (29/9/7), we can celebrate Speed of Light, or c Day
On Sep 07, 2029 (29/9/7), we can celebrate Speed of Light, or c Day
On Jan 04, 2031 (31/4/1), we can celebrate Archimedes' Constant, or Pi Day
On Apr 01, 2031 (31/4/1), we can celebrate Archimedes' Constant, or Pi Day
On Apr 15, 2031 (31/4/15), we can celebrate Archimedes' Constant, or Pi Day
On Jan 01, 2032 (1/01/32), we can celebrate Standard Atmosphere, or atm Day
On Jan 01, 2032 (1/01/32), we can celebrate Standard Atmosphere, or atm Day
On Jan 07, 2032 (1/7/32), we can celebrate Root 3 Day
On Jan 10, 2032 (10/1/32), we can celebrate Standard Atmosphere, or atm Day
On Jul 01, 2032 (1/7/32), we can celebrate Root 3 Day
On Oct 01, 2032 (10/1/32), we can celebrate Standard Atmosphere, or atm Day
On Feb 02, 2036 (2/2/36), we can celebrate Root 5 Day
On Feb 02, 2036 (2/2/36), we can celebrate Root 5 Day
On Jan 03, 2041 (3/1/41), we can celebrate Archimedes' Constant, or Pi Day
On Mar 01, 2041 (3/1/41), we can celebrate Archimedes' Constant, or Pi Day
On Jan 14, 2042 (14/1/42), we can celebrate Pythagoras' Constant, or Root 2 Day
On Aug 31, 2044 (8/31/44), we can celebrate Ideal Gas Constant, or R Day
On Jun 09, 2048 (9/6/48), we can celebrate Faraday Constant, or F Day
On Sep 06, 2048 (9/6/48), we can celebrate Faraday Constant, or F Day
On Mar 17, 2050 (17/3/2050), we can celebrate Root 3 Day
On Jan 09, 2052 (52/9/1), we can celebrate Bohr Radius, or a_0 Day
On Sep 01, 2052 (52/9/1), we can celebrate Bohr Radius, or a_0 Day
On Sep 17, 2052 (52/9/17), we can celebrate Bohr Radius, or a_0 Day
On Jan 09, 2053 (53/9/1), we can celebrate Planck Time, or t_P Day
On Sep 01, 2053 (53/9/1), we can celebrate Planck Time, or t_P Day
On Sep 10, 2053 (53/9/10), we can celebrate Planck Time, or t_P Day
On Oct 09, 2053 (53/9/10), we can celebrate Planck Time, or t_P Day
On Jul 18, 2055 (18/7/55), we can celebrate Planck Charge, or q_P Day
On Mar 07, 2056 (56/7/03), we can celebrate Stefan-Boltzmann Constant, or Theta Day
On Jul 03, 2056 (56/7/03), we can celebrate Stefan-Boltzmann Constant, or Theta Day
On Jan 06, 2060 (1/6/60), we can celebrate Atomic Mass Constant, or m_u Day
On Jan 22, 2060 (60/22/1), we can celebrate Avogadro's Number, or N_A Day
On Feb 02, 2060 (60/2/2), we can celebrate Avogadro's Number, or N_A Day
On Feb 02, 2060 (60/2/2), we can celebrate Avogadro's Number, or N_A Day
On Feb 21, 2060 (60/2/21), we can celebrate Avogadro's Number, or N_A Day
On Feb 23, 2060 (2/23/60), we can celebrate Root 5 Day
On Mar 22, 2060 (22/3/60), we can celebrate Root 5 Day
On Jun 01, 2060 (1/6/60), we can celebrate Atomic Mass Constant, or m_u Day
On Jan 16, 2061 (16/1/61), we can celebrate Planck Length, or l_P Day
On Feb 17, 2065 (2/17/65), we can celebrate Planck Mass, or m_P Day
On Jul 21, 2065 (21/7/65), we can celebrate Planck Mass, or m_P Day
On Feb 06, 2066 (66/2/6), we can celebrate Planck Constant, or h Day
On Mar 07, 2066 (66/7/3), we can celebrate Gravitational Constant, or G Day
On Jun 02, 2066 (66/2/6), we can celebrate Planck Constant, or h Day
On Jun 26, 2066 (66/26/06), we can celebrate Planck Constant, or h Day
On Jul 03, 2066 (66/7/3), we can celebrate Gravitational Constant, or G Day
On Jan 14, 2068 (14/1/68), we can celebrate Planck Temperature, or T_P Day
On May 06, 2070 (5/6/70), we can celebrate Stefan-Boltzmann Constant, or Theta Day
On Jun 05, 2070 (5/6/70), we can celebrate Stefan-Boltzmann Constant, or Theta Day
On Jan 06, 2072 (1/6/72), we can celebrate Proton Mass, or m_p Day
On Jun 01, 2072 (1/6/72), we can celebrate Proton Mass, or m_p Day
On Jun 06, 2073 (6/6/73), we can celebrate Gravitational Constant, or G Day
On Jun 06, 2073 (6/6/73), we can celebrate Gravitational Constant, or G Day
On Jan 08, 2075 (1/8/75), we can celebrate Planck Charge, or q_P Day
On Aug 01, 2075 (1/8/75), we can celebrate Planck Charge, or q_P Day
On Jan 02, 2076 (2/1/76), we can celebrate Planck Mass, or m_P Day
On Feb 01, 2076 (2/1/76), we can celebrate Planck Mass, or m_P Day
On Sep 29, 2079 (29/9/79), we can celebrate Speed of Light, or c Day
On Jan 03, 2080 (1/3/80), we can celebrate Boltzmann Constant, or k_B Day
On Jan 16, 2080 (16/1/80), we can celebrate Golden Ratio, or Phi Day
On Mar 01, 2080 (1/3/80), we can celebrate Boltzmann Constant, or k_B Day
On Jan 27, 2082 (27/1/82), we can celebrate Euler's Number, or e Day
On Jan 04, 2083 (83/1/4), we can celebrate Ideal Gas Constant, or R Day
On Apr 01, 2083 (83/1/4), we can celebrate Ideal Gas Constant, or R Day
On Apr 14, 2083 (83/14/4), we can celebrate Ideal Gas Constant, or R Day
On Aug 09, 2087 (8/9/87), we can celebrate Coulomb's Constant, or k_e Day
On Sep 08, 2087 (8/9/87), we can celebrate Coulomb's Constant, or k_e Day
On Jul 08, 2089 (89/8/7), we can celebrate Coulomb's Constant, or k_e Day
On Aug 07, 2089 (89/8/7), we can celebrate Coulomb's Constant, or k_e Day
On Feb 05, 2091 (5/2/91), we can celebrate Bohr Radius, or a_0 Day
On Mar 05, 2091 (5/3/91), we can celebrate Planck Time, or t_P Day
On Mar 09, 2091 (91/09/3), we can celebrate Electron Mass, or m_e Day
On May 02, 2091 (5/2/91), we can celebrate Bohr Radius, or a_0 Day
On May 03, 2091 (5/3/91), we can celebrate Planck Time, or t_P Day
On Sep 03, 2091 (91/09/3), we can celebrate Electron Mass, or m_e Day
On Sep 10, 2093 (9/10/93), we can celebrate Electron Mass, or m_e Day
On Oct 09, 2093 (9/10/93), we can celebrate Electron Mass, or m_e Day
On Apr 08, 2096 (96/4/8), we can celebrate Faraday Constant, or F Day
On Aug 04, 2096 (96/4/8), we can celebrate Faraday Constant, or F Day
On Feb 09, 2097 (2/9/97), we can celebrate Speed of Light, or c Day
On Sep 02, 2097 (2/9/97), we can celebrate Speed of Light, or c Day
On Jun 06, 2098 (98/06/6), we can celebrate Standard Gravity, or g Day
On Jun 06, 2098 (98/06/6), we can celebrate Standard Gravity, or g Day
On Jan 06, 2102 (1/6/02), we can celebrate Elementary Charge, or e Day
On Jun 01, 2102 (1/6/02), we can celebrate Elementary Charge, or e Day
On Jun 16, 2105 (16/6/05), we can celebrate Atomic Mass Constant, or m_u Day
On Aug 09, 2106 (9/8/06), we can celebrate Standard Gravity, or g Day
On Aug 13, 2106 (13/8/06), we can celebrate Boltzmann Constant, or k_B Day
On Sep 08, 2106 (9/8/06), we can celebrate Standard Gravity, or g Day
On Jan 09, 2109 (9/1/09), we can celebrate Electron Mass, or m_e Day
On Sep 01, 2109 (9/1/09), we can celebrate Electron Mass, or m_e Day
On Jan 03, 2110 (10/1/3), we can celebrate Standard Atmosphere, or atm Day
On Feb 13, 2110 (10/13/2), we can celebrate Standard Atmosphere, or atm Day
On Mar 01, 2110 (10/1/3), we can celebrate Standard Atmosphere, or atm Day
On Jun 08, 2113 (13/8/06), we can celebrate Boltzmann Constant, or k_B Day
On Aug 06, 2113 (13/8/06), we can celebrate Boltzmann Constant, or k_B Day
On Jan 04, 2114 (1/4/14), we can celebrate Pythagoras' Constant, or Root 2 Day
On Jan 04, 2114 (14/1/4), we can celebrate Pythagoras' Constant, or Root 2 Day
On Jan 06, 2114 (14/1/6), we can celebrate Planck Temperature, or T_P Day
On Feb 14, 2114 (14/14/2), we can celebrate Pythagoras' Constant, or Root 2 Day
On Mar 08, 2114 (8/3/14), we can celebrate Ideal Gas Constant, or R Day
On Apr 01, 2114 (1/4/14), we can celebrate Pythagoras' Constant, or Root 2 Day
On Apr 01, 2114 (14/1/4), we can celebrate Pythagoras' Constant, or Root 2 Day
On Jun 01, 2114 (14/1/6), we can celebrate Planck Temperature, or T_P Day
On Aug 03, 2114 (8/3/14), we can celebrate Ideal Gas Constant, or R Day
On Aug 16, 2114 (14/16/8), we can celebrate Planck Temperature, or T_P Day
On Feb 06, 2141 (6/02/2141), we can celebrate Avogadro's Number, or N_A Day
On Jun 02, 2141 (6/02/2141), we can celebrate Avogadro's Number, or N_A Day
On Jan 05, 2176 (2176/5/1), we can celebrate Planck Mass, or m_P Day
On May 01, 2176 (2176/5/1), we can celebrate Planck Mass, or m_P Day
On May 10, 2176 (2176/5/10), we can celebrate Planck Mass, or m_P Day
On Oct 05, 2176 (2176/5/10), we can celebrate Planck Mass, or m_P Day
On Jun 07, 2236 (2236/06/7), we can celebrate Root 5 Day
On Jul 06, 2236 (2236/06/7), we can celebrate Root 5 Day
On Oct 13, 2500 (10/13/2500), we can celebrate Standard Atmosphere, or atm Day
On Jun 06, 2606 (6/6/2606), we can celebrate Planck Constant, or h Day
On Jun 06, 2606 (6/6/2606), we can celebrate Planck Constant, or h Day
On Jul 16, 2621 (16/7/2621), we can celebrate Proton Mass, or m_p Day
On Jan 28, 2718 (2718/28/1), we can celebrate Euler's Number, or e Day
On Feb 08, 2718 (2718/2/8), we can celebrate Euler's Number, or e Day
On Aug 02, 2718 (2718/2/8), we can celebrate Euler's Number, or e Day
On Feb 09, 2997 (2997/9/2), we can celebrate Speed of Light, or c Day
On Sep 02, 2997 (2997/9/2), we can celebrate Speed of Light, or c Day
On Sep 24, 2997 (2997/9/24), we can celebrate Speed of Light, or c Day
On May 09, 3141 (3141/5/9), we can celebrate Archimedes' Constant, or Pi Day
On Sep 05, 3141 (3141/5/9), we can celebrate Archimedes' Constant, or Pi Day
On Jan 07, 3205 (1/7/3205), we can celebrate Root 3 Day
On Jul 01, 3205 (1/7/3205), we can celebrate Root 3 Day
On Jan 01, 3250 (1/01/3250), we can celebrate Standard Atmosphere, or atm Day
On Jan 01, 3250 (1/01/3250), we can celebrate Standard Atmosphere, or atm Day
On Jan 10, 3250 (10/1/3250), we can celebrate Standard Atmosphere, or atm Day
On Oct 01, 3250 (10/1/3250), we can celebrate Standard Atmosphere, or atm Day
On Feb 02, 3606 (2/2/3606), we can celebrate Root 5 Day
On Feb 02, 3606 (2/2/3606), we can celebrate Root 5 Day
On Jan 03, 4159 (3/1/4159), we can celebrate Archimedes' Constant, or Pi Day
On Mar 01, 4159 (3/1/4159), we can celebrate Archimedes' Constant, or Pi Day
On Jan 14, 4213 (14/1/4213), we can celebrate Pythagoras' Constant, or Root 2 Day
On Aug 31, 4462 (8/31/4462), we can celebrate Ideal Gas Constant, or R Day
On Jun 09, 4853 (9/6/4853), we can celebrate Faraday Constant, or F Day
On Sep 06, 4853 (9/6/4853), we can celebrate Faraday Constant, or F Day
On Jul 07, 5291 (5291/7/7), we can celebrate Bohr Radius, or a_0 Day
On Jul 07, 5291 (5291/7/7), we can celebrate Bohr Radius, or a_0 Day
On Jul 18, 5545 (18/7/5545), we can celebrate Planck Charge, or q_P Day
On Mar 07, 5670 (5670/3/7), we can celebrate Stefan-Boltzmann Constant, or Theta Day
On Jul 03, 5670 (5670/3/7), we can celebrate Stefan-Boltzmann Constant, or Theta Day
On Jan 04, 6022 (6022/1/4), we can celebrate Avogadro's Number, or N_A Day
On Jan 14, 6022 (6022/14/1), we can celebrate Avogadro's Number, or N_A Day
On Apr 01, 6022 (6022/1/4), we can celebrate Avogadro's Number, or N_A Day
On Dec 14, 6022 (6022/14/12), we can celebrate Avogadro's Number, or N_A Day
On Jan 06, 6053 (1/6/6053), we can celebrate Atomic Mass Constant, or m_u Day
On Jun 01, 6053 (1/6/6053), we can celebrate Atomic Mass Constant, or m_u Day
On Feb 23, 6067 (2/23/6067), we can celebrate Root 5 Day
On Mar 22, 6067 (22/3/6067), we can celebrate Root 5 Day
On Jan 16, 6199 (16/1/6199), we can celebrate Planck Length, or l_P Day
On Feb 17, 6510 (2/17/6510), we can celebrate Planck Mass, or m_P Day
On Jul 21, 6510 (21/7/6510), we can celebrate Planck Mass, or m_P Day
On Jun 09, 6626 (6626/06/9), we can celebrate Planck Constant, or h Day
On Sep 06, 6626 (6626/06/9), we can celebrate Planck Constant, or h Day
On Apr 08, 6673 (6673/8/4), we can celebrate Gravitational Constant, or G Day
On Aug 04, 6673 (6673/8/4), we can celebrate Gravitational Constant, or G Day
On Jan 14, 6833 (14/1/6833), we can celebrate Planck Temperature, or T_P Day
On May 06, 7037 (5/6/7037), we can celebrate Stefan-Boltzmann Constant, or Theta Day
On Jun 05, 7037 (5/6/7037), we can celebrate Stefan-Boltzmann Constant, or Theta Day
On Jan 06, 7262 (1/6/7262), we can celebrate Proton Mass, or m_p Day
On Jun 01, 7262 (1/6/7262), we can celebrate Proton Mass, or m_p Day
On Jun 06, 7384 (6/6/7384), we can celebrate Gravitational Constant, or G Day
On Jun 06, 7384 (6/6/7384), we can celebrate Gravitational Constant, or G Day
On Jan 08, 7554 (1/8/7554), we can celebrate Planck Charge, or q_P Day
On Aug 01, 7554 (1/8/7554), we can celebrate Planck Charge, or q_P Day
On Jan 02, 7651 (2/1/7651), we can celebrate Planck Mass, or m_P Day
On Feb 01, 7651 (2/1/7651), we can celebrate Planck Mass, or m_P Day
On Sep 29, 7924 (29/9/7924), we can celebrate Speed of Light, or c Day
On Jan 16, 8033 (16/1/8033), we can celebrate Golden Ratio, or Phi Day
On Jan 03, 8064 (1/3/8064), we can celebrate Boltzmann Constant, or k_B Day
On Mar 01, 8064 (1/3/8064), we can celebrate Boltzmann Constant, or k_B Day
On Jan 27, 8281 (27/1/8281), we can celebrate Euler's Number, or e Day
On Apr 06, 8314 (8314/4/6), we can celebrate Ideal Gas Constant, or R Day
On Jun 04, 8314 (8314/4/6), we can celebrate Ideal Gas Constant, or R Day
On Aug 09, 8755 (8/9/8755), we can celebrate Coulomb's Constant, or k_e Day
On Sep 08, 8755 (8/9/8755), we can celebrate Coulomb's Constant, or k_e Day
On May 05, 8987 (8987/5/5), we can celebrate Coulomb's Constant, or k_e Day
On May 05, 8987 (8987/5/5), we can celebrate Coulomb's Constant, or k_e Day
On Mar 05, 9106 (5/3/9106), we can celebrate Planck Time, or t_P Day
On May 03, 9106 (5/3/9106), we can celebrate Planck Time, or t_P Day
On Mar 08, 9109 (9109/3/8), we can celebrate Electron Mass, or m_e Day
On Aug 03, 9109 (9109/3/8), we can celebrate Electron Mass, or m_e Day
On Feb 05, 9177 (5/2/9177), we can celebrate Bohr Radius, or a_0 Day
On May 02, 9177 (5/2/9177), we can celebrate Bohr Radius, or a_0 Day
On Sep 10, 9382 (9/10/9382), we can celebrate Electron Mass, or m_e Day
On Oct 09, 9382 (9/10/9382), we can celebrate Electron Mass, or m_e Day
On Mar 05, 9648 (9648/5/3), we can celebrate Faraday Constant, or F Day
On May 03, 9648 (9648/5/3), we can celebrate Faraday Constant, or F Day
On Feb 09, 9792 (2/9/9792), we can celebrate Speed of Light, or c Day
On Sep 02, 9792 (2/9/9792), we can celebrate Speed of Light, or c Day
On May 06, 9806 (9806/6/5), we can celebrate Standard Gravity, or g Day
On Jun 05, 9806 (9806/6/5), we can celebrate Standard Gravity, or g Day

Now it's your turn

What other special numbers would you like to hold holidays for? Measurements of the Earth? Unit conversion factors? Interplanetary distances? The Fibonacci sequence? Try my function with some of your favorites and post them in the comments.

]]>An Ornamental Geometric Inequality
http://feedproxy.google.com/~r/mathworks/moler/~3/7MXX0DSDwOY/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/cleve/collatz_inequality_01.png"/></div><p>I came across this "ornamental geometric inequality" in a tribute to Lothar Collatz.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2015/03/16/an-ornamental-geometric-inequality/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1178Mon, 16 Mar 2015 17:00:14 +0000

I came across this "ornamental geometric inequality" in a tribute to Lothar Collatz.

I mentioned the German mathematician Lothar Collatz in my post in January on his 3n+1 Conjecture. Here is a beautiful inequality from his 1934 paper titled "Ornamental Geometric Inequalities". Consider the subregion of $-4 \le (x,y) \le 4$ for which

]]>Deploying P-code Files
http://feedproxy.google.com/~r/mathworks/pick/~3/ek7xBtAIYFg/
<p>
Jiro's pick this week is deploypcode by Sven.
Sven's no stranger to Pick of the Week. This time, he has created a tool that allows you to easily deploy your MATLAB files
as P-files (protected files). P-files... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/03/13/deploying-p-code-files/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5890Fri, 13 Mar 2015 13:00:09 +0000

Sven's no stranger to Pick of the Week. This time, he has created a tool that allows you to easily deploy your MATLAB files
as P-files (protected files). P-files are obfuscated versions of your MATLAB files and are good for hiding the contents of your MATLAB
source code. Note, however, that P-files are not encrypted and should not be considered secure.

Sven's deploypcode can take the contents of a particular folder and P-code all of the MATLAB files, including files in subfolders. The function
provides a number of additional features, as well. Two of my favorites are the ability to keep the help text and the ability
to copy additional files for easy deployment. For example, let's say that I have some MATLAB files and a data file.

]]>PicksMATLAB Data Types for Dates and Time, Part II
http://feedproxy.google.com/~r/mathworks/loren/~3/yXd6EWPPDYk/
<p>Once again we're going to hear from guest blogger Andrea Ho, who works for the MATLAB Documentation team here at MathWorks.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/03/11/matlab-data-types-for-dates-and-time-part-ii/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1109Wed, 11 Mar 2015 12:24:45 +0000

Once again we're going to hear from guest blogger Andrea Ho, who works for the MATLAB Documentation team here at MathWorks.

Last time, I introduced new data types for representing date and time data in MATLAB: datetime for representing specific points in time, duration for representing exact differences between points in time, and calendarDuration for representing flexible amounts of time such as months. Today, we'll look at how these data types help you to work with time zones, daylight saving time, and dates in languages other than English.

Interest Payments from a Swiss Bank

Last time, we looked at bank interest payments over the course of a year, without paying attention to the bank location. Let's do a similar analysis, but this time, our bank is located in Switzerland and we are in Phoenix, Arizona.

Suppose we opened a bank account on January 1, 2013, and the account pays interest monthly on the last day of each month. Create a sequence of dates in 2013 to represent when the bank account pays interest. The dateshift function can shift a date and create a sequence based on the initial date in one line.

p0 = datetime(2013,1,1);
p = dateshift(p0,'end','month',0:11)'

Suppose our records indicate that we received interest payments at either 3 AM or 4 AM, depending on the month of the year. We can input this information to our datetime vector p, by modifying its Hour property.

Our bank in Zurich, Switzerland, has sent us a file named interest_2013.txt that contains interest amounts for each month of the year.

For English-language dates and times, you can use the Import Tool to import your data graphically. However, the dates in our file are written in German:

Date,Interest
31. Januar 2013 12:00,1.93
28. Februar 2013 12:00,4.28
31. März 2013 12:00,4.82
30. April 2013 12:00,1.23
31. Mai 2013 12:00,5.89
30. Juni 2013 12:00,2.26
31. Juli 2013 12:00,3.84
31. August 2013 12:00,5.82
30. September 2013 12:00,2.51
31. Oktober 2013 12:00,2.99
30. November 2013 12:00,6.17
31. Dezember 2013 12:00,2.65

For non-English language dates, we must import the data using functions. We will use readtable to read the file data into a table. Use the %D specifier to read the first column of data as datetime values and specify the format of the date strings within curly braces. The 'DateLocale' name-value pair argument tells MATLAB how to interpret the date strings in the file. The value of DateLocale is a combination of a lowercase language code and an uppercase country code. For example, 'de_CH' indicates German-language dates in Switzerland.

T = readtable('interest_2013.txt',...'Format','%{dd.MMMM yyyy HH:mm}D %f','DateLocale','de_CH')

It appears that the bank's payments were made at noon each month, while our own records indicate otherwise. What happened? Our original times are based in Phoenix, whereas the imported times from the bank are based in Zurich.

By default, datetime values are not associated with any time zone. That is, they are "unzoned". Because we want to compare dates and times across geographic locations and therefore time zones, we should associate each datetime array with a time zone. You can learn more about time zones here.

Let's set the time zone for our Phoenix-based datetime vector, p. Valid values for the TimeZone property include names of time zone regions from the IANA Time Zone Database. Then, adjust the display format to show the UTC offset for each value.

The third payment date (March 31) falls on a Sunday. Suppose the bank shifted this payment date earlier by two days so that it would not fall on a weekend.

What happens if we subtract a duration of 2 standard, 24-hour long days to the corresponding payment date?

new_date = z(3) - days(2)

new_date =
29-Mar-2013 11:00 +0100

Notice that the new payment time is now 11 AM instead of noon because a daylight saving time shift occurred in Zurich on March 30, 2014.

To account for a daylight saving time shift in a calendar calculation, we should subtract 2 calendar days from the original payment date.

new_date = z(3) - caldays(2)

new_date =
29-Mar-2013 12:00 +0100

Now the payment time is at noon, and is consistent with the payment times during the rest of the year. If we calculate the exact duration between the original payment date and the new payment date, we will see that the difference is not 48 hours (two 24-hour long days). Rather, the difference is 47 hours due to the daylight saving time change.

new_date - z(3)

ans =
-47:00:00

Here's the big takeaway: When a datetime value is associated with a time zone that observes daylight saving time, it is affected by daylight saving time changes. To correctly account for such time changes, make sure you use calendar durations instead of durations in calendar calculations involving days, months, or years. calendarDuration values account for non-constant lengths of time such as the length of a day, which is not always equal to 24 hours in some time zones.

Quarterly Statistics

Now we will calculate the mean of the interest payments for each quarter. Just like how the day function can extract day names for each datetime value in an array, the quarter function lets us extract the quarter number.

q = quarter(z)

q =
1
1
1
2
2
2
3
3
3
4
4
4

Now we can calculate statistics for the interest payments associated with each unique value of of q.

for ii = 1:4
X = sprintf('Quarter %d',ii);
disp(X)
tf = q==ii;
m_quarter = mean(T.Interest(tf))
% T.Interest is a column containing interest values from% the text file we imported earlierend

Now we've seen how to create datetime values and how to import date and time data as datetime values that can account for time zones and daylight saving time. But what if you currently have data in the form of serial date numbers? Fear not, there's an easy way to convert your existing data to more convenient datetime values. Let's start with an array of serial date numbers:

dn = (735600:31:735755)'

dn =
735600
735631
735662
735693
735724
735755

Use the datetime function to easily convert the date numbers to datetime values.

An array of serial date numbers cannot account for a time zone, but you can add time zone information to the datetime array.

p.TimeZone = 'America/New_York';

You can even export the datetime array in a different language. Let's translate the dates into German and then export them to a text file.

C = cellstr(p,'dd. MMMM yyyy','de_DE')

C =
'01. Januar 2014'
'01. Februar 2014'
'04. März 2014'
'04. April 2014'
'05. Mai 2014'
'05. Juni 2014'

T = table(C,rand(6,1),'VariableNames',{'Date','InterestRate'})

T =
Date InterestRate
__________________ ____________
'01. Januar 2014' 0.2785
'01. Februar 2014' 0.54688
'04. März 2014' 0.95751
'04. April 2014' 0.96489
'05. Mai 2014' 0.15761
'05. Juni 2014' 0.97059

writetable(T,'myfile.txt')

Your Thoughts?

Have you tried using the datetime, duration, and calendarDuration data types? Let us know what you think by leaving a comment here.

]]>Why “Where’s Waldo”?
http://feedproxy.google.com/~r/mathworks/loren/~3/yNVCTM6c4QQ/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/findingWaldo.gif"/></div><p>Today I'd like to introduce guest blogger Seth DeLand who works for the MATLAB Product Marketing team here at MathWorks. Today, Seth will be discussing some different optimization techniques that he used for developing a strategy for the puzzle-books "Where's Waldo?".... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/03/06/why-wheres-waldo/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1125Fri, 06 Mar 2015 16:16:34 +0000

Today I'd like to introduce guest blogger Seth DeLand who works for the MATLAB Product Marketing team here at MathWorks. Today, Seth will be discussing some different optimization techniques that he used for developing a strategy for the puzzle-books "Where's Waldo?".

There was a nice post on Randal Olson’s blog last week about techniques for the popular puzzle-books “Where's Waldo?” (originally “Where’' Wally?” in the UK).

As a kid, I spent hours going through my collection of Where's Waldo? books, so I was inspired by Randal's post to do some of my own analysis. It would be great if I could show that in the last 20 years, I’ve learned enough to have a little better strategy for finding Waldo.

I started by grabbing the data for the Waldo locations in every book, and plotting the points:

The next thing I wanted to do was try to reproduce some of Randal’s results. The traveling salesman problem (What is the shortest path that goes through all of the points?) is always fun to dig into. I used the same start and end points for my path as in the original post, but took a different approach for finding the solution. Instead of using a genetic algorithm, I chose to formulate this as a binary programming problem, similar to what is done in this example. This approach has 2 advantages:

It's really fast for problems of this size (on the order of 100 stops)

We can prove we’ve actually found the shortest possible path

Let's look at the results:

[order,pathLength,output] = shortestPathAToB([pts.X pts.Y],1,67);
hold on
plot(pts.X(order),pts.Y(order))

Examining the output from the INTLINPROG optimization solver, we see that the Absolute Gap is 0, meaning we have found the shortest possible path.

output.absolutegap

ans =
0

A different take on the shortest path

Solving the traveling salesman problem is nice, but when you’re actually searching for Waldo you’re looking at a region instead of a point. It’s kind of like you’re applying a circular mask to the page, and you can see everything within that circle. So this got me thinking: what if we tried to find the shortest path, such that if you translated a circle of radius “r” along that path, you would cover all of the Waldo locations? This makes the problem a bit more complex, since we have to compute distances from all of the Waldo locations to arbitrary line segments making up the path.

I chose to formulate it as a nonlinear programming problem, where again the objective is to minimize the total distance, and we add a nonlinear constraint that says all of the Waldo points must be within “r” of some point on the path. The solution from the last problem serves as a good starting point. We then have enough to turn FMINCON loose on the problem. To make things interesting, I solved for 4 different values of “r”.

r = 0.25:0.25:1;
pts = pts(order,:);
npts = length(pts.X);
x0 = [pts.X; pts.Y]; % Initial guess at solution
lb = zeros(2*npts,1); % All variables >= 0
ub = [13*ones(npts,1); 8*ones(npts,1)]; % x-variables are <= 13, y-variables are <= 8
options = optimoptions('fmincon',...'Algorithm','sqp',...'Display','none',...'TolFun',1E-3,...'MaxFunEvals',1E5);
xopt = zeros(length(x0),length(r));
fval = zeros(1,length(r));
parfor i = 1:length(r)
noncon = @(x) nonlcon(x,[pts.X,pts.Y],r(i));
[xopt(:,i),fval(i)] = fmincon(@objfcn,x0,[],[],[],[],lb,ub,noncon,options);
end
figure;
for i = 1:length(r)
subplot(2,2,i);
drawPath(xopt(:,i),pts,r(i));
title({['View Radius = ' num2str(r(i))]; ...
['Path Length = ' num2str(fval(i))]});
end

Analyzing the solution

As expected, as the value of “r” increases, the length of the path decreases. Here’s an animation of what it would look like traversing this path for a view radius of 1:

The path length shrinks from a length of 58.8 when we explicitly visit each stop, to 31.2 when we have a radius of 1. So we have managed to shorten our search path by quite a bit. Looking at the path, I would describe it as a Z then an uppercase Lambda, or:

$$Z\Lambda$$

That makes it easy enough to remember. All of the code I used for this analysis is available in this GitHub repository.

Who has a different idea for finding Waldo? Have any of you image processing gurus tried tackling this problem? Let us know here.

]]>Plot Subfunctions
http://feedproxy.google.com/~r/mathworks/pick/~3/OpiJFZxXdXE/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/pick/Richard/main_FunctionPlotter/main_FunctionPlotter_01.png"/></div><p>
Posted by Richard Ruff , March 6, 2015
Richard is an Application Engineer at MathWorks focused on the Embedded Coder product for code generation, primarily in the
... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/03/06/plot-subfunctions/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5876Fri, 06 Mar 2015 14:00:58 +0000

Posted by Richard Ruff , March 6, 2015

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

My pick this week is the utility function plot_subfun for visualizing the subroutines in a file via a node map.

Have you ever been given a MATLAB file only to spend many hours reading through it to understand the calling structure? If
so, this function can save a lot of time by providing a snapshot of the calling structure, thus allowing you to quickly see
the interaction between the various subroutines. As the saying goes, a picture is worth a thousand words and this function
gives you the picture from the thousand words or five thousand words or ...

The File Exchange submission includes an extensive demo script, 'demo_subfun.m', that provides very detailed examples for
using plot_subfun. A basic example is shown here:

plot_subfun('plot_subfun');

There are a number of options that allow the user to customize the results generated by plot_subfun. These options are:

-ends: hide functions with one parent and no kids

-ext: show calls to external functions in a separate list

-hide: hide named functions and their children

-kids: same as -hide, but only hides the children

-list: display list of all functions and their calls to command line

-nest: hide nested functions

-noplot: do not plot the figure

-trim: trims the output to just the visible functions

-unhide: hidden functions are shown in grey instead

-unused: also show functions that are not called

-unusedonly: only show functions that are not called

I find plot_subfun particularly useful when developing custom UIs in MATLAB as I can quickly see what callbacks are being
used and any utility routines that are shared. In addition, it can help point out where missing connections are such as needing
to call a specific routine if the user selects an option such as switching between tabs on a UI and updating the data being
displayed.

While there is an option to list the extenal functions that are called, I would like to see the option to include a set of
files to analyze in order to see the dependency between the functions in multiple files. For Example, I have two files mainUI.m
and fillUIPanel.m, I’d like to see this syntax work:

plot_subfun({'mainUI.m', 'fillUIPanel.m'});

Comments

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

]]>PicksWelcome R2015a!
http://feedproxy.google.com/~r/SethOnSimulink/~3/YQmqcvBNohA/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/seth/2015Q1/timeScope.png"/></div><p>MATLAB R2015a just went live, let's explore some of the new features.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/03/05/welcome-r2015a/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4460Fri, 06 Mar 2015 04:39:44 +0000MATLAB R2015a just went live, let's explore some of the new features.

Dashboard Library

In R2015a, you will find a new Dashboard section in the Simulink Library Browser with many controls and displays:

When you drop one of those in your model, double-click on it to associate it with any signal. Click play and begin to use those new controls and displays.

Area

I have to admit, this feature should have been implemented a long time ago! I have seen many models where Annotations or empty Subsystems were used to group blocks in a model. This is what the new Area does.

Algebraic loop highlighting

Until now, it could be difficult to identify blocks in a algebraic loop. In R2015a, if there is an algebraic loop in your model, use Simulink.BlockDiagram.getAlgebraicLoops(bdroot) to highlight the blocks and signals involved in the loop.

Resettable Subsystem

For the regular readers, you might remember a post from last year titled What about a Resettable Subsystem?. Well... we made it! Woohoo!

This block can be very useful to reset the states of blocks that do not have reset ports, like the transfer function:

Simplification of Sample Time Specification

If you open the dialog of many basic math blocks, you will notice that something is gone when compared to previous releases. We hope this will make managing sample times simpler. If you want to specify multiple rates in your model, use the Rate Transition block.

When you add a Scope to your model, click the Try Time Scope button:

The Scope will be replaced by a new Time Scope very similar to the Time Scope that has been available in the DSP System Toolbox for many releases. This new scope has many features like cursors, statistics and triggers.

Now it's your turn

Take a few minutes to review the release notes and let us know what is your favorite feature in R2015a by leaving a comment here.

]]>Direct Access to Seismological Data using MATLAB
http://feedproxy.google.com/~r/mathworks/loren/~3/aYz9qHHsaBI/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/iris03sm.png"/></div><p>I've recently been engaged in several events of potential interest to the geophysics community. In December 2014, a few of us from MathWorks attended the AGU meeting (<a rel="nofollow" target="_blank" href="http://www.agu.org">American Geophysical Union</a>).... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/03/03/direct-access-to-seismological-data-using-matlab/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1130Tue, 03 Mar 2015 20:47:12 +0000

I've recently been engaged in several events of potential interest to the geophysics community. In December 2014, a few of us from MathWorks attended the AGU meeting (American Geophysical Union).

And last week, I was invited to give a webinar to the IRIS (Incorporated Research Institutions for Seismology) community. The video for MATLAB for Analyzing & Visualizing Geospatial Data is here.

Next, I would like to pass on some timely information from Chad Trabant and Robert Weekly, from IRIS, about accessing seismological data via IRIS sites and services.

Traditionally, seismologists spend considerable time and effort collecting and organizing the data they wish to use for a study. In an effort to reduce this prep work, the IRIS Data Management Center (DMC) are constantly trying to remove barriers associated with data collection and formatting to deliver data to users that is easy for them to use and will give researchers more time to actually perform their research.

The DMC operates one of the largest repositories of openly available seismological time series data. Data managed by the DMC is retrieved from stations around the world, with international partners contributing data, in addition to data collected by U.S. institutions. These data are primarily from passive source sensors and are what most folks are familiar with when they see earthquake signals.

One of the primary charges of the DMC is to collect, archive, and distribute data to not only university researchers and teachers, but to anyone with internet access. A significant amount of work at the DMC is dedicated to building software packages that enable people to access the data stored in the DMC archives through a variety of web-enabled clients.

Accessing data from the DMC in MATLAB

For MATLAB users, DMC has developed irisFetch.m which provides seamless access to data stored at the DMC, as well as other data centers that implement standard web services.

This bit of MATLAB code operates as an abstraction layer to the DMC IRIS-WS Java library, which in turn fetches data and information from a data center via web services. Importantly, a user does not need to know anything about Java or networking: syntax for the functions provided by irisFetch.m should be familiar to users of the MATLAB language, and all data is returned as MATLAB structures.

The irisFetch module provides access to the following:

Time series data, mostly seismological (ground motion)

Seismological station metadata, including instrument response

Event (earthquake) parameters, e.g., location, time, magnitude

Help getting started using irisFetch.m is available in the online manual for irisFetch. Examples are also included below.

A standard web services foundation, accessing data at multiple centers

The core web services for delivering time series, related metadata and event (earthquake) information used by irisFetch have been standardized by the FDSN2. DMC has purposely designed irisFetch, and the Java library it uses, to work with any services compliant with the specification. This means that irisFetch may be used to retrieve data from any data center supporting these services, however, not all data centers have implemented every FDSN-standardized service, so functionality with irisFetch is limited to what any given center supports. A list of FDSN data centers and the services they support is here.

A short aside on seismic data nomenclature and norms

In the world of passive source seismology (SEED format specifically) time series are identified with four identifiers: network, station, location and channel. Most of the data available from the DMC is recorded continuously, sometimes for decades, so a time range must also be specified to form a minimal request. In the irisFetch functions, time series are selected using these identifiers. In most cases, seismological data is returned in units of "digital counts", which is often proportional to velocity or acceleration within a certain frequency range.

Examples

The irisFetch module comprises three separate methods that can used to retrieve station metadata, time series, and event information. They are irisFetch.Stations, irisFetch.Traces and irisFetch.Events, respectively. Below are some examples of how each of these components of irisFetch could be used for common data retrieval tasks.

Example 1

Plot time series for the 2011 Tohoku-Oki earthquake. The example below uses irisFetch.Traces to retrieve time series data for global GSN stations to plot 1 hour of data (bandpass filtered from 1 - 5 Hz) following the Tohoku-Oni earthquake of March 2011.

sta = {'MAJO','WMQ','KAPI','CMB'};
for i=1:length(sta)
tohoku_tr(i) = irisFetch.Traces('_GSN',sta{i},'00','BHZ','2011-03-11 15:10:00',...
'2011-03-11 16:10:00','verbose');
end

colors = brighten(lines(numel(tohoku_tr)),-0.33);
figure(1)
for i=1:numel(tohoku_tr)
subplot(4,1,i)
tr = tohoku_tr(i);
data = (tr.data - mean(tr.data)) ./ tr.sensitivity;
wn1 = bandfilt_freq1/tr.sampleRate;
wn2 = bandfilt_freq2/tr.sampleRate;
[f1,f2] = butter(bandfilt_order,[wn1 wn2],'stop');
data = filter(f1,f2,data);
sampletimes = linspace(tr.startTime,tr.endTime,tr.sampleCount);
plot(sampletimes, data, 'color', colors(2,:),'LineWidth',1);
datetick;
if i==1
title('Tohoku-Oni Earthquake, 11-March-2011','FontSize',14)
end
if i~=numel(tohoku_tr)
set(gca,'xticklabel','');
else
xt = get(gca,'XTick'); xtl = get(gca,'XTickLabel');
xt = xt(1:2:end); xtl = xtl(1:2:end,:);
set(gca,'XTick',xt,'XTickLabel',xtl);
end
set(gca,'YTick',[],'YTickLabel','','TickDir','out',...
'TickLength',[.005 .015],'FontSize',12);
ylabel([tohoku_tr(i).network '.' tohoku_tr(i).station])
hold on;
end

Example 2

Plot a map of Transportable Array (TA) network stations and recent earthquakes in Oklahoma. This example uses irisFetch.Events and irisFetch.Stations to retrieve hypocenter data and station metadata, respectively, to make a simple plot in MATLAB.

figure(2)
plot([ok_ev.PreferredLongitude],[ok_ev.PreferredLatitude],'r.',...
'MarkerSize',10)
hold on
plot([ok_sta.Longitude],[ok_sta.Latitude],'b^','MarkerFaceColor','b')
deg2in=6.0/diff(xlim);
set(gca,'units','inches','pos',...
[1.5 1 diff(xlim)*cosd(35)*deg2in diff(ylim)*deg2in],...
'FontSize',12,'TickDir','out')
l = legend('Earthquake','TA station','Location','NorthEast');
set(l,'FontSize',12)
xlabel('Longitude')
ylabel('Latitude')
title(['M<' num2str(maxmag) ' -- ' datestr(starttime(1:10),1) ' to ' ...
datestr(endtime(1:10),1) ' -- ' ...
num2str(numel(ok_ev)) ' total events'],'FontSize',14)

Example 3

Retrieving data from other datacenters. The versatility of the irisFetch methods allow users to access data from web services hosted by other datacenters, as long as they conform to the specification set by the FDSN (as discussed above). This example uses data obtained from the USGS fdsn-event service and some simple commands from the MATLAB Mapping Toolbox.

figure(3)
% Here we use the MATLAB Mapping Toolbox
h = worldmap('world');
geoshow(h,'landareas.shp','FaceColor',[.3 .3 .3])
l = geoshow(h,[glob_ev.PreferredLatitude],[glob_ev.PreferredLongitude],'DisplayType','point');
set(l,'Marker','.')
mlabel('GLineWidth',.2,'FontSize',12,...
'MLabelLocation',90,'MLabelParallel','south');
title(['Global M>' num2str(minmag_glob) ' events since ' datestr(tstart(1:10),1)],...
'FontSize',14)

Trouble accessing your seismological data?

How do you get access to the seismological data you analyze – irisFetch – or something else? What are the stumbling blocks you encounter when accessing or working with seismological data? Let us know here.

In my previous post on iterative refinement I showed the need for an inner product routine that is more accurate than double precision. This post is about such a function.

Example with double precision

The example I am going to use is contrived so that the first and third terms in the inner product exactly cancel each other, leaving the much smaller second term to arise from the ashes.

format longe
x = [1 1/3 1]'
y = [1 3e-9 -1]'

x =
1.000000000000000e+00
3.333333333333333e-01
1.000000000000000e+00
y =
1.000000000000000e+00
3.000000000000000e-09
-1.000000000000000e+00

Computing the inner product with a conventional MATLAB statement shows the intended difficulty.

dot2p = x'*y

dot2p =
1.000000082740371e-09

The result should be 1.000000000000000e-09. We're getting only about half of the significant digits correct.

Of course, the problem is that the second intermediate sum is

s2 = 1 + 1/3*3e-9

s2 =
1.000000001000000e+00

That's OK in decimal, but not in binary. There is not enough room in one double precision floating point number to store the bits that are going to be needed when that leading one is cancelled by the third term in the sum.

Triple precision

I do not have a full blown triple precision arithmetic package by any means. I have just enough to do this one task. Here is the basic idea. A double precision floating point number has 14 hexadecimal digits in its fraction. I can use single and double to break a double into 6 high order hex digits and 8 low order hex digits, like this.

format hex
x = 1/3
xhi = double(single(x))
xlo = x - xhi

x =
3fd5555555555555
xhi =
3fd5555560000000
xlo =
be45555556000000

Two doubles with more than half of their low order bits equal to zero can be multiplied together with no roundoff error. For example

pihi = double(single(pi))
pio3hi = xhi*pihi

pihi =
400921fb60000000
pio3hi =
3ff0c1524860a920

That trailing zero in pio3hi is an indication that the result is exact.

Additions are not exact, when the two numbers differ by several orders of magnitude. This fact will eventually be the limiting factor of our inner product routine.

dot3p

My inner production routine is both accumulated, which means it uses extra precise arithmetic, and extended, which means an extra scalar is added to the result using the extra precision. You can download this function from here.

dbtype dot3p

1 function s = dot3p(x,y,s)
2 % DOT3P s = dot3p(x,y,s) Triple precision extended inner product.
3 % s = x'*y + s for vectors x and y and scalar s.
4
5 shi = double(single(s));
6 slo = s - shi;
7 for k = 1:length(x)
8 xhi = double(single(x(k)));
9 xlo = x(k) - xhi;
10 yhi = double(single(y(k)));
11 ylo = y(k) - yhi;
12 tmp = xhi*yhi;
13 zhi = double(single(tmp));
14 zlo = tmp - zhi + xhi*ylo + xlo*yhi + xlo*ylo;
15
16 tmp = shi + zhi;
17 del = tmp - shi - zhi;
18 shi = double(single(tmp));
19 slo = tmp - shi + slo + zlo - del;
20 end
21 s = shi + slo;

Example with triple precision

Let's run my example with dot3p in the debugger and look at some intermediate results.

% dbstop 16 dot3p% dbstop 20 dot3p% dot3p(x,y,0)

The variables shi and slo carry the sum in triple precision. The first time through the loop there are no roundoff errors, and shi and slo are set to 1.0 and 0.0.

Let's look at the second pass through the loop. Halfway through, at statement 16.

We can see that xhi is the first six hex digits of 1/3 and xlo is the remaining digits. Similarly, yhi is the first six hex digits of 3e-9 and ylo is the remaing digits. zhi is the first six hex digits of the xhi*yhi and zlo is a crucial correction term.

Stopping at the end of the second pass through the loop, at statement 20.

K>> format hex
K>> tmp, del, shi, slo
tmp =
3ff000000044b830
del =
0000000000000000
shi =
3ff0000000000000
slo =
3e112e0be826d694
K>>
K>> format long e,
K>> tmp, del, shi, slo
tmp =
1.000000001000000e+00
del =
0
shi =
1
slo =
9.999999999999999e-10

tmp is 1.0 plus some of the bits of 1.0e-10. del is zero because it is not needed in this example. It is involved when the terms vary over an even wider range. shi is exactly 1.0, which is the high order part of the evolving sum. And slo has become 1.0e-10 to full double precison accuracy.

On the third time through the loop there will be no roundoff errors, shi will be completely cancelled, and slo will bear the full responsibility for providing the final answer.

Of course, this example is contrived and unusual. Ordinarily, we can expect some cancellation (otherwise there would be no need for an accumulated inner product), with the high order part losing at least a few digits and the low order part filling them in.

residual3p

With a dot product in hand, it is easy to write the residual function. Here the extended feature is essential because we expect extreme, if not total, cancellation when the right hand side is subtracted from the dot product.

type residual3p

function r = residual3p(A,x,b)
% RESIDUAL3p Triple precision residual, A*x - b.
% r = residual3p(A,x,b) for matrix A, vectors x and b.
m = size(A,1);
r = zeros(m,1);
for k = 1:m
r(k) = dot3p(A(k,:),x,-b(k));
end
end

]]>About that dress …
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/dEdyQtUb-7Q/
<div class="overview-image"><img src="http://blogs.mathworks.com/steve/files/dress_small.jpg" class="img-responsive attachment-post-thumbnail wp-post-image" alt="dress_small"/></div><p>This afternoon my wife looked at me with an expression indicating she was convinced that I had finally gone completely and irrevocably crazy.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2015/03/01/about-that-dress/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1277Mon, 02 Mar 2015 04:29:56 +0000

This afternoon my wife looked at me with an expression indicating she was convinced that I had finally gone completely and irrevocably crazy.

I think it was because I was carrying from room to room a flashlight, a candle, a lighter, and a chessboard scattered with colored craft sticks and puff balls.

"It's about that dress," I said.

Yes, that dress. #TheDress, as it is known on Twitter.

Three nights ago I idly opened Twitter to see what was up. The Internet was melting down over: a) llamas, and b) a picture of a dress. (I have nothing more to say here about llamas.)

This is the picture as it was originally posted on Tumblr:

Some people see this dress as blue and black. Some see it as white and gold. Each group can't understand why the others see it differently.

By Friday afternoon, a myriad of explanations had popped up online and on various news outlets. Mostly, I found these initial attempts to be unsatisfying, although some better explanations have been published online since then.

Initially I didn't want to write a blog about this, because (as I often proclaim) color science makes my brain hurt. But I do know a little bit about how color scientists think, having worked with several, having read their papers and books, and having implemented their methods in software. So, here is my interpretation of this unusual visual phenomenon. It's in three parts:

The influence of illumination

The phenomenon of color constancy

How two different people could arrive at dramatically different conclusions about the color of that dress.

Let's start with the influence of illumination. Here is a small portion of a picture that I took today.

"Sage green," my wife said.

And here's a portion of a different picture.

"That's yellow," came the answer.

The truth: these two colors are from the same location of the same object. Here are the two original images with the locations marked.

Image A

Image B

The chessboard and other objects in these pictures are the same. The difference between the two images is caused entirely by the different light sources used for each one. Just for fun here are the colors of the puff ball on the upper right from two different pictures. (Remember, these are pixels from the exact same spot on the same object!)

The color of the light arriving at the camera depends not only on the color of the object, but also on the nature of the illumination. As you can see in the colored patches above, changing the illumination can make a big difference. So you cannot definitively determine the dress color solely from close examination of the digital image pixels.

Look at Image A again. What is the color of the index card?

Most people would call it white. If you look at just a chunk of pixels from the center of the card, though, it looks like a shade of green.

People have an amazing ability to compensate automatically and unconsciously for different light sources in a scene that they are viewing. If you looked at the same banana under a bright fluorescent light, and in candle light, and in the shade outdoors under a cloudy sky, you would see the banana as having the same color, yellow, each time. That is true even though the color spectrum of the light coming from the banana is actually significantly different in these three scenarios. Our ability to do this is called color constancy.

Our ability to compensate accurately for illumination depends on having familiar things in the scene we are viewing. It can be the sky, the pavement, the walls, the grass, the skin tones on a face. Almost always there is something in the scene that anchors our brain's mechanism that compensates for the illumination.

Now we come back to the photo of the dress. The photo completely lacks recognizable cues to help us compensate for illumination. Our brain tries to do it anyway, though in spite of the lack of cues. The different reactions of people around the world suggest that there are two dramatically different solutions to the problem.

Consider the diagram below. It illustrates two scenes. In the first scene, on the left, a hypothetical blue and black dress is illuminated by one source. In the second scene, on the right, a hypothetical white and gold dress is illuminated by a second source. By coincidence, both scenes produce the same light at the camera, and therefore the two photographs look the same.

Note: I tweaked the diagram above at 02-Mar-2015 14:50 UTC to clarify its interpretation.

For some people, their visual system is jumping to illumination 1, leading them to see a blue and black dress. For others, their visual system is jumping to illumination 2, leading them to see a white and gold dress.

On Friday afternoon, I heard a report that many people in the white and gold camp, when they see a version of the photograph that includes a woman's face, immediately change their perception of the dress to blue and black. This perceptual shift persists even when they view the original photograph again. This demonstrates that the presence of an object with a familiar color can significantly alter our perception of colors throughout the scene.

If you zoom way in and examine the pixels on the dress in the original image, you'll see that they are blue.

So what kind of illumination scenario could cause people to perceive this as white? I asked Toshia McCabe, a MathWorks writer who knows more than I do about color and color systems. She thinks the dress picture was edited from another one in which the dress was underexposed. As a result, the "blue coincidentally looks like a white that is in shadow (but daylight balanced)." In other words, light from a white object in shadowed daylight can arrive as blue light to the camera. So if your eye sees blue pixels, but your brain jumps to the conclusion that the original scene was taken in shaded daylight, then your brain might decide you are looking at a white dress.

For the record, my wife sees a white and gold dress on a computer monitor, but she sees blue and brown when it is printed. I see blue and brown.

]]>UncategorizedFor Each
http://feedproxy.google.com/~r/mathworks/pick/~3/SBJVeE6UKfc/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/pick/Sean/mainForEach/fileresults.png"/></div><p>
Sean's pick this week is For-Each by Jeremy Hughes.
My pick this week is a toolbox from our development organization for a for-each loop construct in MATLAB. A for-each... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/02/27/for-each/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5868Fri, 27 Feb 2015 14:00:15 +0000

My pick this week is a toolbox from our development organization for a for-each loop construct in MATLAB. A for-each loop allows you to iterate over a set or dimension without having to worry about indexing.

MATLAB actually already has this built in for single elements of numeric arrays:

% Display each element in xfor thisfile = [1 pi 8 -3]
disp(thisfile)
end

1
3.1416
8
-3

This could've course been written with an index like below.

% Display each element in x with an index
thisfile = [1 pi 8 -3];
for ii = 1:numel(thisfile)
disp(thisfile(ii))
end

1
3.1416
8
-3

But what if we want to traverse every column or slice of a three dimensional array? Or each element in a cell array? Or
each combination of two arrays? This is where Jeremy's for-each construct comes in handy since you won't need to worry about indexing.

Let's take a bunch of csv files and make copies of them with just the rows of data we care about. These files contain energy
outage information for energy outages in America. Being from New England, I just want the records that correspond to the
North East.

% Scrape the data directory for all csv files
datadir = '.\Data\';
files = dir([datadir '*.csv']);
filenames = {files.name};
disp(filenames)

Columns 1 through 4
'2002Data.csv' '2003Data.csv' '2004Data.csv' '2005Data.csv'
Columns 5 through 8
'2006Data.csv' '2007Data.csv' '2008Data.csv' '2009Data.csv'
Columns 9 through 12
'2010Data.csv' '2011Data.csv' '2012Data.csv' '2013Data.csv'

% Make a directory for output
outputdir = 'NorthEastData';
mkdir(outputdir)
% Loop over each filename, read, extract, writefor thisfile = each(filenames)
% Reach each file
Temp = readtable([datadir thisfile]);
Temp.Region = categorical(Temp.Region); % Convert to categorical% Identify and extract the northeast records
northeast = Temp.Region == 'NorthEast';
Output = Temp(northeast,:);
% Make new full file name
[~,name,ext] = fileparts(thisfile);
outputfile = fullfile(pwd,outputdir,[name '_NorthEast' ext]);
% Write out
writetable(Output,outputfile)
end

Now we can check to make sure it did what we expect.

Personally, I see this being really useful for working with sets of files or images, like I have done above, or for needle
in the haystack type problems where I want to find something in a bunch of combinations and stop as soon as I do.

Advanced Programming Challenge!

And now for a challenge:

I used the informal interface to this toolbox, i.e. the each function, to traverse each file. There is also a formal interface that provides the ability to define your own iterators so that you can traverse objects of your own class or anything else. I will give some MathWorks goodies to anyone who will
write their own eachFile iterator that I could use above to recreate what I have.

]]>PicksLet’s Code! Make a Cody Video
http://feedproxy.google.com/~r/mathworks/desktop/~3/l3-RiY6v3eE/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/community/files/solved_it.png"/></div><p>Ever seen a “Let’s Play” video, where somebody documents their own activity as they play a video game? They can be a lot of fun to watch, even when the author isn’t very good at the game. And they can be very instructive if you’re learning how to play that... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/community/2015/02/26/lets-code-make-a-cody-video/">read more >></a></p>http://blogs.mathworks.com/community/?p=3168Thu, 26 Feb 2015 23:21:26 +0000Ever seen a “Let’s Play” video, where somebody documents their own activity as they play a video game? They can be a lot of fun to watch, even when the author isn’t very good at the game. And they can be very instructive if you’re learning how to play that particular game. Effectively, you get to “watch over someone’s shoulder” as they play.

Learning by watching videos has become standard technique for students of everything from piano to surgery. And increasingly people are learning to code by watching other people code.

On this site we have lots of clever people solving MATLAB problems on Cody. I would love to see somebody post a video to YouTube that shows them actually in the process of solving a Cody problem. Or maybe just describing after-the-fact how they did solve it. Or even analyzing, sports-commentator style, the various techniques used to solve a particular problem. If you make a “How I Solved It” video for Cody leave a link here and tell us about it!

]]>UncategorizedSignoff
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/lBPyVRfS9AU/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2015/02/26/signoff/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201502/1590/62009828001_4084745180001_358-Signoff.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 01:15</span>
</div>
</a></div><p>After 14 awesome years at MathWorks, I am changing careers. It has been great working and learning with all of you. If you are new to this blog, go back and watch my other videos, they are still as useful today as they always were.
Keep watching here to see who... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2015/02/26/signoff/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1572Thu, 26 Feb 2015 13:34:51 +0000
Keep watching here to see who ends up taking over this space!

Thank you,

Doug

]]>Scaling Market Basket Analysis with MapReduce
http://feedproxy.google.com/~r/mathworks/loren/~3/6mQijLisENg/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/microsoftweb.png"/></div><p>In <a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/loren/?p=1089">an earlier post</a>, today's guest blogger Toshi Takeuchi gave us an introduction to Market Basket Analysis. This week, he will discuss how to scale this technique using <a rel="nofollow" target="_blank" href="http://www.mathworks.com/help/matlab/mapreduce.html"><tt>MapReduce</tt></a> to deal with larger data.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/02/25/scaling-market-basket-analysis-with-mapreduce/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1105Wed, 25 Feb 2015 12:52:15 +0000

In an earlier post, today's guest blogger Toshi Takeuchi gave us an introduction to Market Basket Analysis. This week, he will discuss how to scale this technique using MapReduce to deal with larger data.

R2014b was a major update to MATLAB core functionality and one of the several new exciting features to me was MapReduce. I was primarily interested in Market Basket Analysis to analyze clickstream data and I knew web usage data extracted from web server logs would be very large.

MapReduce was developed to process massive datasets in a distributed parallel computation, and it is one of the key technologies that enabled Big Data analytics.

MapReduce is made up of mappers and reducers. Mappers read data from file storage one chunk at a time and parse data to generate key-value pairs. Then reducers will receive those key-value pairs by key and process values associated with those values. Therefore what you need to do is:

Use datastore to designate data sources

Define mapper and reducer functions

Use mapreduce with datastore, the mapper and reducer to process data

Store the processed data for further analysis

Though mappers and reducers perform fairly simple operations, you can chain them together to handle more complex operations. In the Apriori algorithm, the most time-consuming steps are the generation of transactions and 1-itemset data. So let's use MapReduce to address these bottlenecks.

We start by setting up datastore. In this example, we are going to use a fairly small CSV file on a local drive, but datastore can

read a data that is too large to fit in a memory in a single computer, or

read files on multiple locations on a cluster, including those on a Hadoop distributed file system (HDFS), with appropriate add-on products.

It is important to start out with a small subset of the data to prototype and test your algorithm before you use it on really big data. MATLAB makes it really easy to prototype your algorithm on your local machine and then scale it up to the cluster or cloud later.

If you think of a visitor as a shopper, you can think of pages visited as items in the shopping cart (a transaction). A visitor can visit the same page multiple times, but such repeated visits are not factored in itemset counting.

One of the important considerations in designing a MapReduce algorithm is to minimize the number of keys you generate. For this reason, a good starting point would be to group items by transactions by using VisitorCookieID as the key, because we have a finite set of visitors but they can visit a larger number of pages.

type transactionMapper

function transactionMapper(data, info, intermKVStore)
tid = data.VisitorCookieID;
item = data.Page;
% get unique tids
u_tid = unique(tid);
% iterate over the data chunk to map multiple items to a unique tid
items = cell(size(u_tid));
for i = 1:length(tid)
row = ismember(u_tid,tid{i});
items{row}{end+1} = item{i};
end
% use addmulti to speed up the process
addmulti(intermKVStore, u_tid, items)
end

The mapper will then pass this key-value pair to the reducer.

type transactionReducer

function transactionReducer(key, value, outKVStore)
items = {};
% concatenate items from different mappers for the same key
while hasnext(value)
items = [items, getnext(value)];
end
% eliminate duplicates
u_items = unique(items);
% save the data to a key-value store
add(outKVStore, key, u_items);
end

The reducer receives key-value pairs by key, and merges multiple cell arrays with the same key into a single cell array and removes any duplicates. We can then store the result in a new datastore.

We now know that all items in a transaction are unique. All we need to do is to count the number of times an item appears in the transactions to see how many transactions contained that item. The pages were stored as a value in the previous step, so we simply need to retrieve just those and count their contents.

type oneItemsetMapper

function oneItemsetMapper(data, info, intermKVStore)
% keys are in a cell array
keys = data.Value{1};
% create a cell array of 1's
values = num2cell(ones(size(keys)));
% save the data to a key-value store
addmulti(intermKVStore, keys, values)
end

The mapper passes each instance of a page as 1.

type oneItemsetReducer

function oneItemsetReducer(key, value, outKVStore)
count = 0;
while hasnext(value)
count = count + getnext(value);
end
add(outKVStore, key, count);
end

The reducer then collects the counts of a given page and sums them. Now let's run this job and read the completed results into memory.

% Get 1-itemsets
oneItemset_ds = mapreduce(transaction_ds, @oneItemsetMapper, @oneItemsetReducer);
% Read the result into memory
oneItemsets = readall(oneItemset_ds);
disp(oneItemsets(655:659,:))

Now we are ready to feed the transactions and oneItemsets data to findFreqItemsets, which also accepts a table of 1-itemsets as an optional input. The code is available if you go to the earlier post.

Let's generate frequent itemsets based on a minimum support threshold 0.02, which means we see the same pattern among at least 2% of the visitors.

minSup = 0.02;
fprintf('Processing dataset with minimum support threshold = %.2f\n...\n', minSup)
[F,S,items] = findFreqItemsets(transactions,minSup,oneItemsets);
fprintf('Frequent Itemsets Found: %d\n', sum(arrayfun(@(x) size(x.freqSets,1), F)))
fprintf('Max Level : k = %d\n', length(F))
fprintf('Number of Support Data : %d\n\n', length(S))

Processing dataset with minimum support threshold = 0.02
...
Frequent Itemsets Found: 151
Max Level : k = 4
Number of Support Data : 4107

Generate Rules

This step is no different from the example in the earlier post. Let's use a minimum confidence threshold 0.6.

Now we can visualize the rules we generated. This sample dataset is very tiny, and the number of rules are limited.

conf = arrayfun(@(x) x.Conf, rules); % get conf as a vector
lift = arrayfun(@(x) x.Lift, rules); % get lift as a vector
sup = arrayfun(@(x) x.Sup, rules); % get support as a vector
colormap cool
scatter(sup,conf,lift*5, lift, 'filled')
xlabel('Support'); ylabel('Confidence')
t = colorbar('peer',gca);
set(get(t,'ylabel'),'String', 'Lift');
title('Sample Data')

Microsoft Web Data

Unfortunately we can’t share this sample data for you to try, but it is easy to adapt this code using publicly available dataset such as Anonymous Microsoft Web Data Data Set from UCI Machine Learning Repository. This data is preprossed from raw clickstream logs, and we need to process it back to the raw format we used in sample data. You can see the code used to process this dataset below.

When you apply Market Basket Analysis to this dataset, you should get something like this.

Rule #1
{Windows 95} => {Windows Family of OSs}
Conf: 0.91, Lift: 6.46

Rule #18
{Windows Family of OSs, isapi} => {Windows95 Support}
Conf: 0.64, Lift: 11.66

The sample code also includes MapReduce steps on a reduced dataset. I reduced the dataset because my MapReduce code is not optimal for this dataset and runs extremely slowly. My code assumes that there are more pages than visitors, and we have more visitors than pages in this Microsoft dataset.

If you would like to use MapReduce on the full Microsoft dataset, I strongly recommend that you wrote your own MapReduce code that is more optimized for this dataset.

Summary and the Next Step

Now we prototyped the algorithm and tested it. I can see if I am going to get the insight I am looking for with this process. Once I am happy with it, I still need to figure out where to store the larger data on a cluster or cloud. It's more of a business challenge than a technical one, and I am still in the middle of that.

Programming MapReduce in MATLAB is very straight forward and easy. The source data was local in this example, but you can also use it for a larger dataset that sits across multiple locations, such as HDFS files. You can prototype your MapReduce algorithm locally, then change the configuration to scale up to a larger dataset.

In this example, MapReduce was used only for the initial data processing, and the rest was still done in memory. As long as the processed result can fit in the memory, we can use this approach.

However, if the processed data gets larger, then we need to make more use of MapReduce in other steps in the Apriori alogrithm.

The key is to adapt the algorithm to parallel processing. In the current incarnation, you have a single thread in the candidate pruning stage.

Instead, we subdivide the dataset into several chunks and complete the entire process through rule generation in each. We need to adjust the minimum support to account for the reduction of transaction counts in subsets if we do so. Then we can combine the final output. This may not provide the same result as the single thread process, but it should be fairly close.

Your thoughts?

Do you see ways in which you might leverage MapReduce to handle larger data analysis projects? Let us know here.

Appendix - Process Microsoft Web Data

Here is the code I used to process the Microsoft dataset.

type processMicrosoftWebData.m

%% Let's load the dataset
% The UCI dataset is not in the raw log format, but in a special
% non-tabular ASCII format. We need to process the data into a format we
% can use. This is not a typical process you do when you work with actual
% raw clickstream data.
% clear everything
clearvars; close all; clc
% the URL of the source data
filename = 'anonymous-msweb.data';
% if te file doesn't exist, download of the UCI website
if exist(filename,'file') == 0
url = 'http://archive.ics.uci.edu/ml/machine-learning-databases/anonymous/anonymous-msweb.data';
filepath = websave('anonymous-msweb.data',url);
end
% accumulators
aid = []; % attribute id
vroot = {}; % vroot - page name
url = {}; % url of the page
vid = []; % visitor id
visits = []; % aid of vroots visited by a visitor
VisitorCookieID = {}; % same as vid but as string
% open the file
fid = fopen(filename);
% read the first line
tline = fgetl(fid);
% if the line contains a string
while ischar(tline)
% if the line contains attribute
if strcmp(tline(1),'A')
c = textscan(tline,'%*s%d%*d%q%q','Delimiter',',');
aid = [aid;c{1}];
vroot = [vroot;c{2}];
url = [url;c{3}];
% if the line contains case
elseif strcmp(tline(1),'C')
user = textscan(tline,'%*s%*q%d','Delimiter',',');
% if the line contains vote
elseif strcmp(tline(1),'V')
vid = [vid;user{1}];
vote = textscan(tline,'%*s%d%*d','Delimiter',',');
visits = [visits; vote{1}];
VisitorCookieID = [VisitorCookieID;['id',num2str(user{1})]];
end
tline = fgetl(fid);
end
% close the file
fclose(fid);
% sort attributes by aid
[~,idx] = sort(aid);
aid = aid(idx);
vroot = vroot(idx);
url = url(idx);
% populate |Page| with vroot based on aid
Page = cell(size(visits));
for i = 1:length(aid)
Page(visits == aid(i)) = vroot(i);
end
% create table and write it to disk if it doesn't exist
transactions = table(VisitorCookieID,Page);
if exist('msweb.transactions.csv','file') == 0
% just keep the first 318 rows to use as sample data
transactions(319:end,:) = []; % comment out to keep the whole thing
writetable(transactions,'msweb.transactions.csv')
end
%% Association Rule Mining without MapReduce
% Since we already have all necessary pieces of data in the workspace, we
% might as well do the analysis now.
% get unique uid
uniq_uid = unique(vid);
% create a cell array of visits that contains vroots visited
transactions = cell(size(uniq_uid));
for i = 1:length(uniq_uid)
transactions(i) = {visits(vid == uniq_uid(i))'};
end
% find frequent itemsets of vroots from the visits
minSup = 0.02;
fprintf('Processing dataset with minimum support threshold = %.2f\n...\n', minSup)
[F,S,items] = findFreqItemsets(transactions,minSup);
fprintf('Frequent Itemsets Found: %d\n', sum(arrayfun(@(x) size(x.freqSets,1), F)))
fprintf('Max Level : k = %d\n', length(F))
fprintf('Number of Support Data : %d\n\n', length(S))
% generate rules
minConf = 0.6;
rules = generateRules(F,S,minConf);
fprintf('Minimum Confidence : %.2f\n', minConf)
fprintf('Rules Found : %d\n\n', length(rules))
% plot the rules
conf = arrayfun(@(x) x.Conf, rules); % get conf as a vector
lift = arrayfun(@(x) x.Lift, rules); % get lift as a vector
sup = arrayfun(@(x) x.Sup, rules); % get support as a vector
colormap cool
scatter(sup,conf,lift*5, lift, 'filled')
xlabel('Support'); ylabel('Confidence')
t = colorbar('peer',gca);
set(get(t,'ylabel'),'String', 'Lift');
title('Microsoft Web Data')
% display the selected rules
selected = [1,2,3,5,18];
for i = 1:length(selected)
fprintf('Rule #%d\n', selected(i))
lenAnte = length(rules(selected(i)).Ante);
if lenAnte == 1
fprintf('{%s} => {%s}\nConf: %.2f, Lift: %.2f\n\n',...
vroot{rules(selected(i)).Ante(1)},vroot{rules(selected(i)).Conseq},...
rules(selected(i)).Conf,rules(selected(i)).Lift)
elseif lenAnte == 2
fprintf('{%s, %s} => {%s}\nConf: %.2f, Lift: %.2f\n\n',...
vroot{rules(selected(i)).Ante(1)},vroot{rules(selected(i)).Ante(2)},...
vroot{rules(selected(i)).Conseq},rules(selected(i)).Conf,rules(selected(i)).Lift)
end
end
%% MapReduce
% My MapReduce code is not well suited for this dataset and runs extremely
% slow if you use the whole dataset. I will just use just a subset for
% demonstration purpose. If you want to try this on the full dataset, you
% should write your own MapReduce code optimized for it.
% clear everything
clearvars; close all; clc
% set up source datastore
source_ds = datastore('msweb.transactions.csv');
disp(preview(source_ds))
% step 1: Group items by transaction
transaction_ds = mapreduce(source_ds, @transactionMapper, @transactionReducer);
transactions = readall(transaction_ds);
% step 2: Generate 1-itemsets
oneItemset_ds = mapreduce(transaction_ds, @oneItemsetMapper, @oneItemsetReducer);
oneItemsets = readall(oneItemset_ds);

]]>Downsampling polygons (Part 2)
http://feedproxy.google.com/~r/mathworks/pick/~3/buxCUdcL7es/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/pick/jiro/potw_decimatepoly/decimate_performance.png"/></div><p>
Jiro's pick this week is Decimate Polygon by Anton Semechko.
Just a little over a year ago, I wrote a blog post on downsampling polygons. I highlighted Peter's Polygon Simplification. His entry helped me out in one of my side projects I... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/02/20/downsampling-polygons-part-2/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5859Fri, 20 Feb 2015 14:00:22 +0000

Just a little over a year ago, I wrote a blog post on downsampling polygons. I highlighted Peter's Polygon Simplification. His entry helped me out in one of my side projects I was working on at the time. Soon after, I received an email from a
File Exchange contributor letting me know about another entry for accomplishing the task of downsampling polygons. I appreciate
people letting me know about these entries I've overlooked. When I select a file to highlight in my posts I try to do a search
on similar entries, but I don't always catch them all.

Let's see how DecimatePoly works. I'm going to use the same image I used for my previous post.

im = imread('amoeba.png');
imshow(im)

As before, I will convert this to black and white and trace out the boundaries using bwboundaries from the Image Processing Toolbox. Then I will overlay the largest boundary onto the original image.

im2 = im2bw(im);
boundaries = bwboundaries(~im2,'noholes')
largest = boundaries{1};
hold on
plot(largest(:,2),largest(:,1),'r','LineWidth',2)
hold off

With DecimatePoly, I can downsample this polygon by specifying either the maximum acceptable offset from the original boundary or the fraction
of the total number of points to retain. For example, if I want to ensure the downsampled boundary to be off by maximum of
.1 pixels,

b1 = DecimatePoly(largest, [.1 1]);

# of verts perimeter area
in 1695 1948.91 27992.50
out 926 1948.91 27992.50
-----------------------------------------------------
change -45.37% -0.00% 0.00 %

I can see that I have about 45% fewer points while maintaining the same perimeter and area. If I loosen my constraint,

b2 = DecimatePoly(largest, [1 1]);

# of verts perimeter area
in 1695 1948.91 27992.50
out 473 1863.07 28046.00
-----------------------------------------------------
change -72.09% -4.40% 0.19 %

I have over 70% fewer points while not sacrificing much change in perimeter or area.

So why is having fewer points preferred? Well, Anton has included a nice demo with his entry that shows improved performance
in in-polygon tests. This is one of the plots from the demo.

He looked at how different number of polygon points affected the computation time and the accuracy of in-polygon tests. The
horizontal axis (Percent simplification) represents the degree of downsampling. The lower the percent simplification, the
greater the downsampling. The Dice coefficient (blue) shows how similar the results are between the downsampled polygon and the full polygon. The closer the value is to
1, the more similar the results. The right axis (red) shows the relative speed up of the computation time compared to the
full polygon. We can see that even at around 10% simplification (downsampled by 90%), we are able to get over 98% similarity
and 10x speed up.

Comments

Do you have a need to downsample polygons? What is your use case? Let us know here or leave a comment for Anton.

]]>PicksManaging Conflicts in Simulink Projects
http://feedproxy.google.com/~r/SethOnSimulink/~3/y2KoLG_Jzfc/
<div class="overview-image"><img src="http://blogs.mathworks.com/seth/files/feature_image/initialSetup.png" class="img-responsive attachment-post-thumbnail wp-post-image" alt="Merging Simulink model changes from two users"/></div><p>This week we are welcoming guest blogger Tim Hosey from the Simulink Project development team.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/02/19/managing-conflicts-in-simulink-projects/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4428Fri, 20 Feb 2015 00:25:03 +0000This week we are welcoming guest blogger Tim Hosey from the Simulink Project development team.

Introduction

In this blog post I will highlight how Simulink Project can be used to manage conflicts that occur when two users submit changes to the same Simulink model using a source control tool.

A conflict arises when two users make changes to the same file. The first user to submit their changes to the source control system experiences no problem, but the second user’s changes are considered to be in conflict with the first user’s changes; so the source control tool prevents the second user from committing their file until they have merged their changes into the first user’s version of the file.

Environment set up

In order to show how conflict resolution can be performed with Simulink Project I have set up a modified version of our Airframe demo.

I created two working copies called demo1 and demo2 which are checked out from an SVN repository created by running the above demo. User 1 will use the folder demo1 and user 2 will use the folder demo 2.

Here is a screenshot of what the project looks like when user 1 opens their copy of the project in the working copy demo1.

The project contains a model called vertical_channel.slx which looks like this:

User 1 modifies this model to change the sign of both the sum block and the feedback gain, so that it now looks like this:

User 1 then commits his changes using the modified files view.

This trivial change was easy to make and commit to the SVN repository, but as we are about to see when user 2 makes a change to his copy of this file, submitting the change is going to be much more difficult.

User 2 hasn’t performed any update so his project looks the same as user 1’s project before he made his change.

User 2 opens the vertical channel model and adds a scope to his model so that it looks like this:

When user 2 tries to commit they get the following error:

This error tells him that he doesn’t have the latest version of vertical_channel.slx so he can’t commit his changes to it. To get the latest version he performs an update. This causes a problem because the changes in his working copy conflict with the changes made between the version of the file that has and the latest version in the repository submitted by user 1. User 2’s change is now in conflict with the change submitted to the repository by user 1. The project marks this file as conflicted user 2’s project.

Because we registered .slx files as binary with SVN, the SVN update did not attempt change the content of the file vertical_channel.slx, it still contains user 2’s addition of the scope to the original model.

User 2 can use Simulink Project to resolve the conflict as follows:

First the user needs to view the conflict, he can do this by right clicking on the file and selecting ‘view conflicts’ in the Simulink Project’s file context menu. If Simulink Report Generator is installed a Simulink XML comparison will be used to show the conflicts between user 1’s version of the model and user 2’s version.

User 2 can then easily select and merge the changes from user 1 that they want to include in their version of the file.

Once user 2 is happy with the resultant merged file he can right click on the file in the project’s file context menu and select ‘mark conflict resolved’. Now the conflicted file is marked as modified and can be submitted to the SVN repository.

When user 1 performs an update he now gets the merged model in his working copy:

Now it's your turn

Do you often need to merge Simulink models? Let us know by leaving a comment here.

]]>Iterative Refinement for Solutions to Linear Systems
http://feedproxy.google.com/~r/mathworks/moler/~3/SudiZnQxZTk/
<p>Iterative refinement is a technique introduced by Wilkinson for reducing the roundoff error produced during the solution of simultaneous linear equations. Higher precision arithmetic is required for the calculation of the residuals.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2015/02/16/iterative-refinement-for-solutions-to-linear-systems/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1166Mon, 16 Feb 2015 17:00:28 +0000

Iterative refinement is a technique introduced by Wilkinson for reducing the roundoff error produced during the solution of simultaneous linear equations. Higher precision arithmetic is required for the calculation of the residuals.

The first research paper I ever published, in 1967, was titled "Iterative Refinement in Floating Point". It was an analysis of a technique introduced by J. H. Wilkinson almost 20 years earlier for a post processing cleanup that reduces the roundoff error generated when Gaussian elimination or a similar process is used to solve a system of simultaneous linear equations, $Ax = b$.

The iterative refinement algorithm is easily described.

Solve $Ax = b$, saving the triangular factors.

Compute the residuals, $r = Ax - b$.

Use the triangular factors to solve $Ad = r$.

Subtract the correction, $x = x - d$

Repeat the previous three steps if desired.

Complexity

Almost all of the work is in the first step, which can be thought of as producing triangular factors such as $L$ and $U$ so that $A = LU$ while solving $LUx = b$. For a matrix of order $n$ the computational complexity of this step is $O(n^3)$. Saving the factorization reduces the complexity of the remaining refinement steps to something much less, only $O(n^2)$.

The residual

By the early 1960s we had learned from Wilkinson that if a system of simultaneous linear equations is solved by a process like Gaussian elimination or Cholesky factorization, the residual will always be order roundoff error, relative to the matrix and the computed solution, even if the system is nearly singular. This is both good news and bad news.

The good news is that $Ax$ is always close to $b$. This says that the computed solution always comes close to solving the equations, even though $x$ might not be close to the theoretical correct solution, $A^{-1}b$. The pitcher always puts the ball where the batter can hit it, even though that might not be in the strike zone.

The bad news is that it is delicate to compute the residual accurately. If the same precision arithmetic is used to compute the residual that was used to solve the system, the roundoff error involved in computing $r$ will be almost comparable to the effect of the roundoff error present in $x$, so the correction has little chance of being useful.

Inner product

We need to use higher precision arithmetic while computing the residual. Each component of the residual involves a sum of products and then one final subtraction. The exact product of two numbers with a certain precision has twice that precision. With the computers that Wilkinson used, and that I used early in my career, we had access to the full results of multiplications. We were able to write inner product routines that accumulated the sums with twice the working precision.

But it is not easy to write the accumulated inner product routine in modern, portable, machine independent software. It was not easy in Fortran. It is not easy in MATLAB. The original specification of the BLAS, the Basic Linear Algebra Subroutines, was deliberately silent on the matter. Subsequent proposals for extensions of the BLAS have introduced mixed precision, but these extensions have not been widely adopted. So, the key tool we need to implement iterative refinement has not been available.

In my next blog post, I will describe two MATLAB functions residual3p and dot3p. They provide enough of what I call "triple precision" arithmetic to produce an accumulated inner product. It's a hack, but it works well enough to illustrate iterative refinement.

Example

My example involves perhaps the world's most famous badly conditioned matrix, the Hilbert matrix. I won't begin with the Hilbert matrix itself because its elements are the fractions

$$ h_{i,j} = \frac{1}{i+j-1} $$

Many of these fractions can't be represented exactly in floating point, so I would have roundoff error before even getting started with the elimination. Fortunately, the elements of the inverse Hilbert matrix are integers that can be readily generated. There is a function invhilb in the MATLAB elmat directory. I'll choose the 8-by-8. The elements are large, so I need a custom format to display the matrix.

n = 8;
A = invhilb(n);
disp(sprintf('%8.0f %11.0f %11.0f %11.0f %11.0f %11.0f %11.0f %11.0f \n',A))

I am going to try to compute the third column of the inverse of this inverse, which is a column of the Hilbert matrix itself. The right hand side b is a column of the identity matrix. I am hoping to get the fractions x = [1/3 1/4 ... 1/9 1/10].

b = zeros(n,1);
b(3) = 1
format compact
format longe
x = A\b

b =
0
0
1
0
0
0
0
0
x =
3.333333289789291e-01
2.499999961540004e-01
1.999999965600743e-01
1.666666635570877e-01
1.428571400209730e-01
1.249999973935827e-01
1.111111087003338e-01
9.999999775774569e-02

Since I know what x is supposed to look like, I can just eyeball the output and see that I have only about half of the digits correct.

(I used backslash to solve the system. My matrix happens to be symmetric and positive definite, so the elimination algorithm involves the Cholesky factorization. But I'm going to be extravagant, ignore the complexity considerations, and not save the triangular factor.)

Inaccurate residual

Here's my first crack at the residual. I won't do anything special about the precision this time; I'll just use an ordinary MATLAB statement.

r = A*x - b

r =
-9.094947017729282e-13
1.746229827404022e-10
4.656612873077393e-10
1.862645149230957e-08
-1.490116119384766e-08
-2.980232238769531e-08
-3.725290298461914e-08
-1.862645149230957e-08

It's important to look at the size of the residual relative to the sizes of the matrix and the solution.

relative_residual = norm(r)/(norm(A)*norm(x))

relative_residual =
1.147025634044834e-17

The elements in this computed residual are the right order of magnitude, that is roundoff error, but, since I didn't use any extra precision, they are not accurate enough to provide a useful correction.

Now I will use residual3p, which I intend to describe in my next blog and which employs "triple precision" accumulation of the inner products required for an accurate residual.

r = residual3p(A,x,b)

r =
-4.045319634826683e-12
1.523381421009162e-10
-9.919851606809971e-10
2.429459300401504e-09
8.826383179894037e-09
-2.260851861279889e-08
-1.332933052822227e-08
-6.369845095832716e-09

Superficially, this residual looks a lot like the previous one, but it's a lot more accurate. The resulting correction works very well.

d = A\r
x = x - d

d =
-4.354403560053519e-09
-3.845999016894392e-09
-3.439925156187715e-09
-3.109578484769736e-09
-2.836169428940436e-09
-2.606416977917484e-09
-2.410777025186154e-09
-2.242253997222573e-09
x =
3.333333333333327e-01
2.499999999999994e-01
1.999999999999995e-01
1.666666666666662e-01
1.428571428571425e-01
1.249999999999996e-01
1.111111111111108e-01
9.999999999999969e-02

I've now got about 14 digits correct. That's almost, but not quite, full double precision accuracy.

Iterate

Try it again.

r = residual3p(A,x,b)

r =
3.652078639504452e-12
-1.943885052924088e-10
2.523682596233812e-09
-1.359348900109580e-08
3.645651958095186e-08
-5.142027248439263e-08
3.649529745075597e-08
-1.027348206505963e-08

Notice that the residual r is just about the same size as the previous one, even though the solution x is several orders of magnitude more accurate.

The correction changed the solution, but didn't make it appreciably more accurate. I've reached the limits of my triple precision inner product.

More accurate residual

Bring in the big guns, the Symbolic Math Toolbox, to compute a very accurate residual. It is important to use either the 'f' or the 'd' option when converting x to a sym so that the conversion is done exactly.

% r = double(A*sym(x,'d') - b)
r = double(A*sym(x,'f') - b)

r =
3.652078639504452e-12
-1.943885052924088e-10
2.523682707256114e-09
-1.359348633656055e-08
3.645651780459502e-08
-5.142027803550775e-08
3.649529478622071e-08
-1.027348206505963e-08

The correction just nudges the last two digits.

d = A\r
x = x - d

d =
-6.846375178532078e-16
-5.828670755614817e-16
-5.162536953886886e-16
-4.533410583885285e-16
-3.965082139594863e-16
-3.747002624289523e-16
-3.392348053271079e-16
-3.136379972458518e-16
x =
3.333333333333333e-01
2.500000000000000e-01
2.000000000000000e-01
1.666666666666667e-01
1.428571428571429e-01
1.250000000000000e-01
1.111111111111111e-01
1.000000000000000e-01

Now, with a very accurate residual, the elements I get in x are the floating point numbers closest to the fractions in the Hilbert matrix. That's the best I can do.

Further reading

There is more to say about iterative refinement. See Nick Higham's SIAM book and, especially, the 2006 TOMS paper by Demmel, Kahan and their Berkeley colleagues. A preprint is available from Kahan's web site.

James Demmel, Yozo Hida, William Kahan, Xiaoye S. Li, Sonil Mukherjee, and E. Jason Riedy, Error Bounds from Extra-Precise Iterative Refinement, ACM Transactions on Mathematical Software 32, pp. 325-351, 2006, <http://www.cs.berkeley.edu/~wkahan/p325-demmel.pdf>

]]>What’s the reaction when we touch?
http://feedproxy.google.com/~r/mathworks/pick/~3/5P6Lko44Xo4/
<div class="overview-image"><img class="img-responsive"/></div><p>
Contact Modeling in SimMechanics
Greg's pick this week is Rolling Ball on Plane by Janne Salomäki.
Janne put together a nice example of modeling the interaction of a ball contacting a plane, rolling on the plane, and then
falling off the plane.
Janne uses SimMechanics to model the geometry and motion of the ball... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/02/13/whats-the-reaction-when-we-touch/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5831Fri, 13 Feb 2015 14:00:21 +0000Contact Modeling in SimMechanics

Greg's pick this week is Rolling Ball on Plane by Janne Salomäki.
Janne put together a nice example of modeling the interaction of a ball contacting a plane, rolling on the plane, and then
falling off the plane.

Janne uses SimMechanics to model the geometry and motion of the ball and the plane. And then applies C-code S-Function to model
the interaction between the ball and the plane. Contact modeling isn’t something SimMechanics can do, so I’m excited to see
an example that shows the basic concept. I also saw an opportunity to compare using C-code and MATLAB Code to solve this particular
problem.

SimMechanics enables six degree of freedom (6DOF) modeling of rigid bodies in Simulink. A nice benefit is it comes with an animation feature
where you can visualize the motion of the bodies.
Bodies are defined in terms of mass, inertia tensor, and location for the center of mass. Connections to other bodies are
defined by coordinate systems applied to the body, and placing constraints on the type of motion permitted between two coordinate
systems. For example, a revolute joint only has one degree of freedom, while a gimbal joint has three degrees of freedom.
In this particular case, the joint between the ball and the plane has six degrees of freedom, which means there is no restriction on the relative motion between
these two bodies. Using this joint enables measurement and actuation of the relative motions between the two bodies, which is essential for modeling the contact interaction.

What is an S-function?

An S-Function is a “system function”. It is a way to create customized blocks in Simulink that can extend its functionality. People use
S-Functions to for all sorts of capabilities:

Interface Simulink models with other software running on your computer

Compile an optimized block algorithm to improve performance

While S-functions can be a little tricky, especially if you are unfamiliar with the API, they interact more or less directly with the Simulink solver engine, which enables development of a wide variety of block
types and capabilities.
Janne uses the C-code S-function to implement the mathematical model for the interaction between the ball and the plane.

Are there other ways to do this?

There are other ways to develop the interaction model. You could of course use basic Simulink blocks, and develop the equations
using Add blocks, Gain blocks, etc.
Another option is to write the interaction as MATLAB-Code in the MATLAB Function block. This would enable you to test the
algorithm in the flexible MATLAB environment. The MATLAB Function block is not designed to support mathematical integration.
So it is often best to pass out the derivatives of the states as outputs, route them through an integrator block, and pass
the result back as inputs to the MATLAB Function block.

How do they compare?

C-code S-Function

MATLAB Function

MATLAB Function

(Debugging
Disabled)

2.15 sec

2.67 sec

2.30 sec

~

+24%

+6.8%

I thought I would try some basic benchmarking to compare some different implementations.
It’s interesting to me that the MATLAB Function block is only about 7% slower. Remember that by default, the MATLAB Function
block converts the MATLAB Code to C-code and then compiles it to essentially an S-Function.

MATLAB Debugging Disabled

Leaving the debugging on certainly provides for a much bigger performance hit, but it also means you can step through the
code using the MATLAB Debugger. This is a bit more difficult to do when you need to debug the C-code S-Function.
I’ll leave it to you to decide if the convenience of the MATLAB Code is worth the hit in performance.

What do I like about this entry?

First, it is a neat, clean, and clear model. Both the model and the code for the S-function are done in a nice straight
forward manner. I like the use of the macros in the C-code S-function to make the code easier to read by hiding some of the
weird API language. What I find so fascinating is that SimMechanics handles all of the appropriate coordinate transformations,
so the algorithm for the contact model is relatively simple without a lot of math to figure out the body relationships.
I realize this approach might not scale well to a more general implementation of contact modeling, but I think it’s a beautiful
and creative example of investigating rigid body dynamics.

Do you need contact modeling in SimMechanics?

Let us know how you would use it. Do you create S-Functions? What’s one thing you like about them? What’s one thing you dislike
about them? You can leave your comments here.

Published with MATLAB® 8.4
]]>Commenting out… in Stateflow!
http://feedproxy.google.com/~r/SethOnSimulink/~3/XqaRatetGRQ/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/seth/2015Q1/displayingNotes.png"/></div><p>I am so happy right now!... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/02/12/commenting-out-in-stateflow/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4432Thu, 12 Feb 2015 19:32:24 +0000I am so happy right now!

Since R2012b, it has been possible to comment out blocks in Simulink. As soon as this feature was released, it immediately became popular among all the users I know.

Soon after, in R2013b, the possibility to comment-through blocks with same numbers of inputs and outputs was introduced. Also very useful!

Commenting out a Stateflow object is as easy as it should be. Select one or more objects and hit Ctrl+Shift+X

Adding text to the commented objects

One useful feature associated with commented Stateflow objects is the possibility to add notes, for example to explain why the object is commented. Click on the % sign on the bottom left corner of the object and a dialog will popup.

The text entered in this dialog will be displayed when you hover the mouse over the object.

Now it's your turn

In my opinion, the ability to comment out Simulink blocks and Stateflow objects, has most significantly affected my model editing workflow in many years.

Have you also included the commenting in your workflows? Let us know by leaving a comment here.

]]>Introduction to the New MATLAB Data Types for Dates and Time
http://feedproxy.google.com/~r/mathworks/loren/~3/T_BS-IBXYRQ/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/datetimePart1_02.png"/></div><p>Today I'd like to introduce a guest blogger, Andrea Ho, who works for the MATLAB Documentation team here at MathWorks. Today, Andrea will discuss three new data types for representing dates and times in MATLAB in R2014b - <a rel="nofollow" target="_blank" href="http://www.mathworks.com/help/matlab/ref/datetime.html"><tt>datetime</tt></a>, <a rel="nofollow" target="_blank" href="http://www.mathworks.com/help/matlab/ref/duration.html"><tt>duration</tt></a>, and <a rel="nofollow" target="_blank" href="http://www.mathworks.com/help/matlab/ref/calendarduration.html"><tt>calendarDuration</tt></a>.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/02/12/introduction-to-the-new-matlab-data-types-for-dates-and-time/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1092Thu, 12 Feb 2015 13:24:42 +0000

Today I'd like to introduce a guest blogger, Andrea Ho, who works for the MATLAB Documentation team here at MathWorks. Today, Andrea will discuss three new data types for representing dates and times in MATLAB in R2014b - datetime, duration, and calendarDuration.

Dates and times are frequently a critical part of data analysis. If they are part of your work, you've likely encountered some of the following limitations of date numbers, date vectors, and date strings:

Serial date numbers are useful for numeric calculations, but they are difficult to interpret or debug.

Date vectors and strings clearly communicate dates and times in a human-readable format, but they are not useful for calculations and comparisons.

Date numbers, date vectors, and strings that represent both points in time and lengths of time (elapsed time) are confusing. For example, does the date vector [2015 1 1 0 0 0] represent January 1, 2015 at midnight, or an elapsed time of 2015 years, 1 month, and 1 day?

Starting in R2014b, datetime, duration, and calendarDuration are data types dedicated to storing date and time data, and are easily distinguishable from numeric values. The three types differentiate between points in time, elapsed time in units of constant length (such as hours, minutes, and seconds), and elapsed time in units of flexible length (such as weeks and months). datetime, duration, and calendarDuration represent dates and times in a way that is easy to read, is suitable for computations, has much higher precision (up to nanosecond precision), and is more memory-efficient. With the new data types, you won't have to maintain or convert between different representations of the same date and time value. Work with datetime, duration, and calendarDuration much like you do with numeric types like double, using standard array operations.

Let's see the new data types in action.

datetime for Points in Time

The datetime data type represents absolute points in time such as January 1, 2015 at midnight. The datetime function can create a datetime variable from the current date and time:

t0 = datetime('now')

t0 =
16-Jan-2015 16:24:30

or from numeric inputs:

t0 = datetime(2015,1,1,2,3,4)

t0 =
01-Jan-2015 02:03:04

whos t0

Name Size Bytes Class Attributes
t0 1x1 121 datetime

If you are interested only in the time portion of t0, you can display just that.

t0.Format = 'HH:mm'

t0 =
02:03

Changing the format of a datetime variable affects how it is displayed, but not its numeric value. Though you can no longer see the year, month and day values, that information is still stored in the variable. Let's turn the display of the date portion of t0 back on.

t0.Format = 'dd-MMM-yyyy HH:mm'

t0 =
01-Jan-2015 02:03

By changing the display format, you can view the value of a datetime variable with up to nanosecond precision!

t0.Format = 'dd-MMM-yyyy HH:mm:ss.SSSSSSSSS'

t0 =
01-Jan-2015 02:03:04.000000000

If you are a frequent user of date string formats, you'll recognize that some of the format specifiers for datetime are different from those for the datestr, datenum, and datevec functions.. The new format syntax allows for more options and is now consistent with the Unicode Locale Data Markup Language (LDML) standard.

Example: Calculate Driving Time

The duration data type represents the exact difference between two specific points in time. Suppose you drove overnight, leaving your home on January 1, 2015 at 11 PM, and arriving at your destination at 4:30 AM the next day. How much time did you spend on the road? (Ignore the possibility of crossing into another time zone; we'll see how datetime manages time zones in a future post.)

dt is a duration variable displayed as hours:minutes:seconds. Like datetime, a duration behaves mostly like a numeric value but displays in a format that shows time well. Let's change the format to view dt as a number of minutes.

dt.Format = 'm'

dt =
330 mins

We can perform common arithmetic operations on date and time variables. Even though dt is currently displayed as a number of minutes, we can add or subtract a duration in any unit, including hours and seconds.

dt2 = dt - hours(1) + seconds(5)

dt2 =
270.08 mins

Example: Verify Payment Schedule

Suppose you opened a bank account on January 1, 2014. The account pays interest monthly on the last day of each month. Let's find the first date on which interest was paid. Is there a solution that does not depend on the account opening date? The dateshift function is useful for shifting dates to the start or end of a unit of time. In our case, we will shift the account opening date to the end of the month.

Yikes, most of these dates don't fall on the last day of the month. A duration of 30 days, with each day being exactly 24 hours long, was not the correct step size to use. The correct step size depends on the month, since some months have 30 days, while others have 31 (or perhaps 28, or 29). What we really want to do is add a series of calendar months to t1.

The calmonths function creates a calendarDuration array containing calendar months. A calendarDuration value represents a flexible or non-constant amount of time such as 1 month. There is no way to know exactly how many days, hours, minutes, or seconds there are in 1 month without knowing which month of the year I am referring to. The length of a calendar month is unknown until you relate it to a specific point in time.

How many days are between each of the payment dates?

dt = caldiff(t,'days')

dt =
28d
31d
30d
31d
30d
31d
31d
30d
31d
30d
31d

Each incremental calendar month represents a different number of days.

Our bank sent us a file named interest.txt that contains interest amounts for each month of the year. We can import this data interactively using the Import Tool.

Alternatively, we can use the readtable and textscan functions to read the file programmatically. Use the %D specifier to read a date or time string and specify the string format within curly braces.

T = readtable('interest.txt','Format','%{dd MMMM yyyy HH:mm}D %f')

T =
Date Interest
_______________________ ________
31 January 2014 12:00 1.93
28 February 2014 12:00 4.28
31 March 2014 12:00 4.82
30 April 2014 12:00 1.23
31 May 2014 12:00 5.89
30 June 2014 12:00 2.26
31 July 2014 12:00 3.84
31 August 2014 12:00 5.82
30 September 2014 12:00 2.51
31 October 2014 12:00 2.99
30 November 2014 12:00 6.17
31 December 2014 12:00 2.65

Even though the text file contains strings, the values in the Date column of the table are not strings. They are datetime values.

If we compare our payment times with the bank's data, we'll see that they are the same.

isequal(t,T.Date)

ans =
1

Now let's calculate the cumulative interest at the end of each month and plot the data. When using plot, working with datetime is much easier than working with serial date numbers.

The date axis is automatically labeled with date strings that are easy to understand. You no longer have to call datetick to do this! If we want to change the format of the dates, we can specify the desired format of the axis ticks in the call to plot.

To create charts containing date and time values using other plotting functions, you still need to work with serial date numbers for the time being. We recognize that this is inconvenient and we are planning improvements. You can learn more here about how to plot datetime and duration values.

What Else?

Notice that in our analysis of travel times and interest payments, there was no need to use serial date numbers, date vectors, or strings at all. datetime, duration, and calendarDuration allow you to manipulate dates and times quickly and more intuitively. The new data types are intended as a replacement for datenum, datevec, and datestr. However, these older functions will continue to exist in MATLAB for the foreseeable future because we know that many people rely on them in existing code.

In a future post, I'll show how datetime, duration, and calendarDuration allow you to work with time zones, daylight saving time, and dates in languages other than English.

Your Thoughts?

Do you see yourself using the new date and time data types with your data? Let us know by leaving a comment here.

]]>CAD APPS
http://feedproxy.google.com/~r/mathworks/pick/~3/B11V5_cFfzg/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/pick/Sean/maincadapps/moment.PNG"/></div><p>
Sean's pick this week is CAD APPS by Larry Silverberg.
For this week's pick I'm going to combine my grad school life in structural engineering with the two recent major... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/02/06/cad-apps/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5825Fri, 06 Feb 2015 14:00:49 +0000

For this week's pick I'm going to combine my grad school life in structural engineering with the two recent major snow storms
we've had in New England over the last few weeks. My house has seen about 42 inches of snow. So let's make a quick model
of the snow load on my back deck.

I'm going to model this as a Frame. Although it's not an ideal frame by any means, the snow load is distributed on it and
not just applied at the ends.

Larry's app allows me to draw the deck interactively (in meters, feet would be nice!)

First I needed to adjust the axes limits:

Next I need to draw the structure. There needs to be a node everywhere I apply a load because there doesn't seem to be a
distributed load option. The deck is about 12ft wide with a 3ft cantilever resting about 8ft off the ground.

Then the hard part, figuring out what the load is in terms of Newtons per meter of deck. Being lazy, I only shoveled the
outer half of the deck where it was easy to throw the snow off.

My wood structural design books are somewhere in my attic at home so I have no idea what percentage of capacity the current
snow load is using. But with another 6:10in of snow coming Monday, I probably oughtta clear off the rest of it. Sigh.

Comments

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

]]>Surprising Results from Modeling Queues with SimEvents
http://feedproxy.google.com/~r/SethOnSimulink/~3/wJtSvW41z2k/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/seth/2014Q4/theoreticalSolution.png"/></div><p>Today I am very happy to welcome guest blogger Ramamurthy Mani who works in our development organization.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/02/04/surprising-results-from-modeling-queues-with-simevents/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4305Thu, 05 Feb 2015 00:58:20 +0000Today I am very happy to welcome guest blogger Ramamurthy Mani who works in our development organization.

Introduction

Have you ever been in a supermarket checkout and wondered why you are in the slowest line? When I saw this article in Wired magazine, the title immediately grabbed my attention. As I read the article, I knew I had to try to put the theory to the test in SimEvents -- a discrete-event simulation add-on to Simulink that is suitable for modeling queuing systems.

Two Types of Queues

I started by building a simple model of a four register supermarket counter. I built two parallel versions of the setup -- one which uses four separate queues and one with a single "serpentine" queue that feeds all four registers.

Let see how this is implemented in SimEvents.

To get the simulation going, I need to model random customers entering the checkout area. I use entities in SimEvents to represent customers and I can generate these entities at random time intervals following an exponential distribution. During generation I can also specify a random duration (also exponentially distributed) that a customer will take to be served at a register by assigning a special attribute to the corresponding entity. I chose the average service time (usually denoted by $1/\mu$) to be 2 mins and the average arrival time (usually denoted by $1/\lambda$) to be 1 min.

I then clone each customer I generate so that I can exercise the two different line configurations identically.

To model the case where four separate queues are feeding the four cash registers, I use a Switch that routes customers to the shortest of the four Queues. Each Queue then feeds a Server representing a checkout register. This Server holds the customer for the amount of time that was setup during generation.

To model the "serpentine" queue, I use a single Queue that feeds the four registers via a Switch that routes customers to a free register when one becomes available.

As you probably noticed, the two models use Timer blocks (the ones with the Stop Watch icons) that keep track of timing statistics for the entities/customers. This allows me to plot the average wait time of customers in the two configurations as well as the minimum and maximum wait times.

Simulation Results

Let us run the model and see what gets plotted.

As you can see, the configuration with the four queues on average results in longer wait times!

Also heartening is that if we compare the results from SimEvents to the theoretical M/M/c queue solution, we can see that this matches pretty well with a result of 2.173913 seconds for the "serpentine" model:

Note that there are some approximations (but no analytical results) for the average arrival time of the 'shortest queue strategy' as well, however I will not get into that for this blog.

Advanced Analysis

Once I had my model setup, also did some other fun experiments with it. I will not go into details in this post, but if you are interested you can download the model and associated files using the link at the end of this post and play with them. To pique your interest, I have summarized some of the interesting things I observed during my experiments:

Measuring "Luckiness:: By assigning an ID to each customer, I tried to measure how often a “lucky” customer who arrives after a specific customer is able to reach a cash register before that “unlucky” customer. With the arrival and service times used above, I was surprised that only 13% of customers were running into an “unlucky” situation. I feel like this happens to me a lot more often in real life!

Impact of system utilization: With an arrival time of 1 minute and service time of 2 minutes, the system utilization is 50%. If we drop the customer arrival time to ~30 seconds, then the system is utilization increases to ~95% and we see that 57% of customers become unlucky. Perhaps this explains what we saw during the recent Christmas season in the grocery store!

Let's try to be smart: : Up to now, customers in my model always picked the shortest queue. By measuring more statistics from SimEvents Queues and Servers, it is possible to always pick the “optimal” queue (where you have a perfect estimate of how long each customer will take at the register) and get the four queues system to perform exactly the same as the serpentine queue. Of course, no human can perfectly estimate the time it will take reach customer at a cash register. To represent this effect, we can add a bit of randomness to include the effect of human optimism bias to our simulation. When I do this, the serpentine queue wins out again.

The dreaded “Price Check on lane 4”: “Price checks” have happened to all of us and these are another factor in the unpredictability of service times. In SimEvents, we can enable a "pause" input to the Servers and also model this effect. Once again when you do this, you begin to see that you are better off using the serpentine queue.

I never thought I would say this, but, based on these experiments, perhaps it is time to give airlines and banks their due – they do setup a single line for transactions.

]]>Origins of Colormaps
http://feedproxy.google.com/~r/mathworks/moler/~3/SMb4ThVkPC8/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/cleve/colormaps_blog_01.png"/></div><p>Steve Eddins has recently posted a <a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/steve/category/colormap/">series in his blog</a> about colormaps. I want to add one more post. With release R2014b, we are retiring <tt>jet</tt> as the default <tt>colormap</tt>, after many years of faithful service. But did you ever wonder where <tt>jet</tt> originated, and how it came to be the default? And did you ever come across colormaps like <tt>pink</tt> and <tt>bone</tt>?... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2015/02/02/origins-of-colormaps/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1162Mon, 02 Feb 2015 17:00:41 +0000

Steve Eddins has recently posted a series in his blog about colormaps. I want to add one more post. With release R2014b, we are retiring jet as the default colormap, after many years of faithful service. But did you ever wonder where jet originated, and how it came to be the default? And did you ever come across colormaps like pink and bone?

MATLAB first began to fully support color and 3D graphics with release 4.0 in 1992. Graphics hardware at the time had a limited amount of memory and so could only display a limited number of colors in any one figure. All color graphics went through color maps. Even the colors in full color photographs had to be quantized into a few bits.

Pseudocolor applies to the display of quantities, such as mathematical functions, that do not have any inherent color. Here the colormap interprets properties such as height, density or velocity as color. For example, the MathWorks logo, the L-shaped membrane, is the solution to a partial differential equation. As such, it does not have any color. A colormap, together with lighting and shading, makes it visually interesting.

jet

The image behind the jetcolormap comes from astrophysical fluid dynamics calculations made at NCSA, the National Center for Supercomputer Applications at the University of Illinois. I found the image on the Internet in the very early days of the Internet, possibly even before the availability of the Mosaic web browser, coincidentally also developed at NCSA in 1993. Astrophysicists Larry Smarr and Mike Norman were among the prominent researchers at NCSA.

The image is available in the MATLAB demos directory.

clear
load flujet
whos
snapnow
caption

Name Size Bytes Class Attributes
X 400x300 960000 double
caption 2x32 128 char
map 64x3 1536 double

caption =
Astrophysical jet from NCSA.
One fluid collides with another.

The 400-by-300 array X contains values of fluid density that have been quantized to "flints" (floating point integers) in the range from to 1 to 64. These serve as indices into the 64-by-3 array map, which is our jetcolormap. The doubles in map range from 0.0 to 1.0 and specify the brightness of the red, green, and blue components of each pixel.

image(X)
colormap(jet)
axis image
axis off

Even though jet was originally given as an array of 64-by-3 numbers, it is actually a piecewise linear function of row index with breaks at 1/8, 3/8, 5/8, and 7/8 of the length. Jet is a variant of the "rainbow" color maps whose faults Steve describes.

rgbploter(@jet)

hot

Here is another simple piecewise linear color map associated with black body radiation. It accurately describes the color obtained when a chisel made from high-carbon tool steel is heated with an oxy-acetylene torch. See the web page of Don Dougan, Sculptor. I learned about this map many years ago from Neil Ostlund, a chemistry professor at the University of Waterloo who started Hypercube, Inc., a company that did molecular modeling on the Intel iPSC.

clf
image(X)
colormap(hot(64))
axis image
axis off

rgbploter(@hot)

pink

Shortly after we introduced colormaps, I was sitting at a terminal demonstrating them to someone and I realized that the maps were just matrices. I could do array operations on them. I typed something like this.

p = sqrt((2*gray(64) + hot(64))/3);

Taking the square root of a color map! Wow, that's pretty strange. But the result was useful. A better name than "pink" would be "sepia". Here's an old photo I've used many times, tinted to make it look old.

clf
load gatlin
image(X)
colormap(pink(64))
axis image
axis off

rgbploter(@pink)

bone

This spine X-ray was high resolution back in the day. To get a colormap that looks like old X-ray file, interchange the blue and red columns of hot and then moderate it with gray.

b = (7*gray(64) + fliplr(hot(64)))/8;
clf
load spine
image(X)
colormap(bone(64))
axis image
axis off

rgbploter(@bone)

spring

We have colormaps named for all four of the seasons. They're all simple linear functions of the row index, like this spring.

The peaks function has proved to be remarkably popular for testing and demonstrating our 3D graphics. It's a modification of something I first saw back in my hypercube days from folks whose names I don't remember at Tektronix in Beaverton, Oregon. They were thinking of nearby Mount St. Helens.

Let's see how the old astrophysical jet looks with the new default colormap. For discussion of parula, including many reader's comments, see Steve's blog.

clf reset
load flujet
image(X)
axis image
axis off

rgbploter(@parula)

By the way, if you want to see a list of all the colormaps, type

]]>Fast, programmatic string search in MATLAB files
http://feedproxy.google.com/~r/mathworks/pick/~3/pW9OAgnobrs/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/pick/jiro/potw_findinm/potw_findinm_screen.png"/></div><p>
Jiro's pick this week is findInM by our very own Brett Shoelson.
If you know Brett, which you probably do if you've spent any time on MATLAB Central, you know about all the useful File Exchange entries he has contributed. It's no... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/01/30/fast-programmatic-string-search-in-matlab-files/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5816Fri, 30 Jan 2015 14:00:41 +0000

If you know Brett, which you probably do if you've spent any time on MATLAB Central, you know about all the useful File Exchange entries he has contributed. It's no surprise he's ranked number 8. Aside from many of the entries related to image processing, he
has uploaded a number of utility functions, and findInM is a must-have. The title of the File Exchange entry starts with "FAST, PROGRAMMATIC string searching...", and that's exactly
what it is. It searches for a string of text in MATLAB files, but it's fast and it's programmatic. There are already a number
of entries in File Exchange that searches for a text within files, including mfilegrep, mgrep, and grep. There is also an interactive way of searching from the toolstrip.

From the description of Brett's entry, findInM "can be much faster than any other method [he's] seen." The way Brett accomplishes this efficient search is by first creating
an index of the MATLAB files in the folders. This step takes some time, but afterwards, the search happens on the index file
and is very efficient.

In this example, I first created an index of my /toolbox folder (and its subfolders) of my MATLAB installation. Then searching
for some text from over 20000 files took less than 10 seconds.

tic
s = findInM('graph theory','toolbox')
toc

SORTED BY DATE, MOST RECENT ON TOP:
s =
'C:\Program Files\MATLAB\R2014b\toolbox\bioinfo\bioinfo\Contents.m'
'C:\Program Files\MATLAB\R2014b\toolbox\bioinfo\biodemos\graphtheorydemo.m'
'C:\Program Files\MATLAB\R2014b\toolbox\matlab\sparfun\Contents.m'
'C:\Program Files\MATLAB\R2014b\toolbox\matlab\sparfun\gplot.m'
'C:\Program Files\MATLAB\R2014b\toolbox\bioinfo\graphtheory\Contents.m'
Elapsed time is 8.334972 seconds.

As a comparison, the interactive "Find Files" tool in MATLAB took over 5 minutes to do the same search.

Thanks for this great tool, Brett! I do have a couple of suggestions for improvement.

Every 7 days, the function prompts the user to see if he/she wants to re-generate the index file. Perhaps this could be somewhat
automated if the indexing process captured the state of the files (file sizes, modified dates). It could automatically recommend
re-generating the index if it notices a change in the state.

The index file is a DOC file. It's easy to open/edit a DOC file. It might be better to use a non-standard extension, so that
it can't be accidentally opened and is easily distinguished from a regular DOC file. For example, in Windows, some folders
with images have a thumbnail database called "Thumbs.db". Perhaps findInM can create a file called "mIndex.mi" or something like that.

Comments

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

]]>PicksIntroduction to Market Basket Analysis
http://feedproxy.google.com/~r/mathworks/loren/~3/HmE_a57cFgc/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/gname.png"/></div><p>You probably heard about the "beer and diapers" story as the often quoted example of what data mining can achieve. It goes like this: some supermarket placed beer next to diapers and got more business because they mined their sales data and found that men often bought those two items together.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/01/29/introduction-to-market-basket-analysis/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1089Thu, 29 Jan 2015 13:20:54 +0000

You probably heard about the "beer and diapers" story as the often quoted example of what data mining can achieve. It goes like this: some supermarket placed beer next to diapers and got more business because they mined their sales data and found that men often bought those two items together.

Today's guest blogger, Toshi Takeuchi, who works in the web marketing team here at MathWorks, gives you a quick introduction to how such analysis is actually done, and will follow up with how you can scale it for larger dataset with MapReduce (new feature in R2014b) in a future post.

I have been interested in Market Basket Analysis not because I work at a supermarket but because it can be used for web usage pattern mining among many applications. Fortunately, I came across a good introduction in Chapter 6 (sample chapter available for free download) of Introduction to Data Mining. I obviously had to give it a try with MATLAB!

Let's start by loading the example dataset used in the textbook. Raw transaction data is turned into a binary matrix representation with 1 meaning the given item was in the given transaction, but otherwise 0.

transactions = {{'Bread','Milk'};...
{'Bread','Diapers','Beer','Eggs'};...
{'Milk','Diapers','Beer','Cola'};...
{'Bread','Milk','Diapers','Beer'};...
{'Bread','Milk','Diapers','Cola'}};
items = unique([transactions{:}]);
T = zeros(size(transactions,1), length(items));
for i = 1:size(transactions,1)
T(i,ismember(items,transactions{i,:})) = 1;
end
disp(array2table(T,'VariableNames',items,'RowNames',{'T1','T2','T3','T4','T5'}))

Market Basket Analysis is also known as Association Analysis or Frequent Itemset Mining. Let’s go over some basic concepts: itemsets, support, confidence, and lift.

The column headers of the table shows all the items in this tiny dataset. A subset of those items in any combination is an itemset. An itemset can contain anywhere from zero items to all the items in the dataset.

A transaction can contain multiple itemsets as its subsets. For example, an itemset {Bread, Diapers} is contained in T2, but not {Bread, Milk}.

Support count is the count of how often a given itemset appears across all the transactions. Frequency of its appearance is given by a metric called support.

Support = itemset support count / number of transactions

itemset = {'Beer','Diapers','Milk'}
% get col indices of items in itemset
cols = ismember(items,itemset);
N = size(T,1); fprintf('Number of transactions = %d\n',N)
% count rows that include all items
supportCount = sum(all(T(:,cols),2));
fprintf('Support count for this itemset = %d\n',supportCount)
itemSetSupport = supportCount/size(T,1);
fprintf('Support = %.2f\n',itemSetSupport)

itemset =
'Beer' 'Diapers' 'Milk'
Number of transactions = 5
Support count for this itemset = 2
Support = 0.40

Association Rules, Confidence and Lift

Association rules are made up of antecedents (ante) and consequents (conseq) and take the following form:

{ante} => {conseq}

A k-itemset where k > 1 can be randomly divided into ante and conseq to form such a rule. Here is an example:

ante = {'Diapers','Milk'}; % pick some subset itemset as |ante|% get items not in |ante|
conseq = setdiff(itemset,ante);
fprintf('Itemset: {%s, %s, %s}\n', itemset{:})
fprintf('Ante : {%s, %s}\n',ante{:})
fprintf('Conseq : {%s}\n',conseq{:})
fprintf('Rule : {%s, %s} => {%s}\n', ante{:},conseq{:})

You can think of this rule as "when diapers and milk appear in the same transaction, you see beer in the same transaction as well at some frequency". How do you determine which randomly generated association rules are actually meaningful?

The most basic measure of rules is confidence, which tells us how often a given rule applies within the transactions that contain the ante.

A given rule applies when all items from both antes and conseqs are present in a transaction, so it is the same thing as an itemset that contains the same items. We can use the support metric for the itemset to compute confidence of the rule.

Confidence = itemset support / ante support

cols = ismember(items,ante);
anteCount = sum(all(T(:,cols),2));
fprintf('Support Count for Ante = %d\n',anteCount)
anteSupport = anteCount/N;
fprintf('Support for Ante = %.2f\n',anteSupport)
confidence = itemSetSupport/anteSupport;
fprintf('Confidence = %.2f (= itemset support / ante support)\n',confidence)

Support Count for Ante = 3
Support for Ante = 0.60
Confidence = 0.67 (= itemset support / ante support)

Another measure of rules is lift. It compares the probability of ante and conseq happening together independently to the observed frequency of such a combination, as a measure of interestingness. If lift is 1, then the probabilities of ante and conseq occurring together is independent and there is no special relationship. If it is larger than 1, then lift tells us how strongly ante and conseq are dependent to each other. We can use respective support metrics to make this comparison.

Lift = itemset support / (ante support x conseq support)

cols = ismember(items,conseq);
conseqCount = sum(all(T(:,cols),2));
fprintf('Support Count for Conseq = %d\n',conseqCount)
conseqSupport = conseqCount/N;
fprintf('Support for Conseq = %.2f\n',conseqSupport)
lift = itemSetSupport/(anteSupport*conseqSupport);
fprintf('Lift = %.2f (= itemset support / (ante support x conseq support))\n',lift)

Support Count for Conseq = 3
Support for Conseq = 0.60
Lift = 1.11 (= itemset support / (ante support x conseq support))

Apriori Algorithm

Now that we know the basic concepts, we can define the goal of our analysis: finding association rules with a sufficient level of support (happens often enough) and confidence (association is strong). Lift can be a secondary measure to score the interestingness of the rules found.

Generate frequent itemsets that clear the minimum support threshold recursively from 1-itemsets to higher level itemsets, pruning candidates along the way

Generate rules that clear the minimum confidence threshold in a similar way

In a brute force method, you would calculate the support and confidence of all possible itemset combinations, but that would be computationally expensive, because the number of candidates grows exponentially.

The apriori algorithm addresses this issue by generating candidates selectively. To get an intuition, think about the frequency of an itemset that contains some infrequent items. That itemset will never be more frequent than the least frequent item it contains. So if you construct your candidates by combining the frequent itemsets only, starting from 1-itemset and continue to higher levels, then you avoid creating useless candidates.

Let's start with generating frequent itemsets and get their support measures. I wrote a couple of MATLAB functions to implement this algorithm.

minSup = 0.6; % minimum support threshold 0.6
[F,S] = findFreqItemsets(transactions,minSup);
fprintf('Minimum Support : %.2f\n', minSup)
fprintf('Frequent Itemsets Found: %d\n', sum(arrayfun(@(x) length(x.freqSets), F)))
fprintf('Max Level Reached : %d-itemsets\n', length(F))
fprintf('Number of Support Data : %d\n', length(S))

Minimum Support : 0.60
Frequent Itemsets Found: 8
Max Level Reached : 2-itemsets
Number of Support Data : 13

When we computed support for each itemset we evaluated, we stored the result in a Map objectS. This is used for rule generation in order to avoid recalculating support as part of the confidence computation. You can now retrieve support for a given itemset with this object, by supplying the string representation of the itemset as the key. Let's try [2,4,6]:

itemset = [2,4,6];
fprintf('Support for the itemset {%s %s %s}: %.2f\n',items{itemset(:)},S(num2str(itemset)))

Support for the itemset {Bread Diapers Milk}: 0.40

This itemset clearly didn't meet the minimum support criteria of 0.60.

Rule Generation Algorithm

We saw earlier that you can generate rule candidates from frequent itemsets by splitting their contents into antes and conseqs, and computing their confidence.

If we generate every possible candidates by such brute force method, it will be very time consuming. Apriori algorithm is also used to generate rules selectively. Let's say that this rule has low confidence.

{Beer, Diapers} => {Milk}

Then any other rules generated from this itemset that contain {Milk} in rule consequent will have low confidence.

Why? because support for those antes will be always greater than the initial ante {Beer, Diapers}, while the support for the itemset (hence also for the rule) remains the same, and confidence is based on the ratio of support between the rule and the ante.

We can take advantage of this intuition by first generating rules with only one item in conseq and drop those that don't meet the minimum criteria, and then merge those conseqs to generate rules with two items in conseq, and so forth.

Now we can generate association rules from the frequent itemsets we generated in the previous step.

minConf = 0.8; % minimum confidence threshold 0.8
rules = generateRules(F,S,minConf);
fprintf('Minimum Confidence : %.2f\n', minConf)
fprintf('Rules Found : %d\n\n', length(rules))
for i = 1:length(rules)
disp([sprintf('{%s}',items{rules(i).Ante}),' => ',...
sprintf('{%s}', items{rules(i).Conseq}),...
sprintf(' Conf: %.2f ',rules(i).Conf),...
sprintf('Lift: %.2f ',rules(i).Lift),...
sprintf('Sup: %.2f',rules(i).Sup)])
end

With a minimum support of 0.6 and a minimum confidence 0.8, we found only one rule that clears those thresholds: {Beer} => {Diapers}. The confidence is 1.00, which means we see diapers in all transactions that include beer, and the lift is 1.25, so this rule is fairly interesting.

Test Example: Congressional Voting Records

As a non-supermarket use case of Market Basket Analysis, the textbook uses the 1984 Congressional Voting Records dataset from the UCI Machine Learning Repository. We can test our code against the result in the textbook.

For the purpose of Market Basket Analysis, you can think of individual members of congress as transactions and bills they voted for and their party affiliation as items.

Now we generate frequent itemsets with a minimum support threshold of 0.3, but we need to use a slightly lower threshold because MATLAB uses double precision numbers.

minSup = 0.2988;
[F,S] = findFreqItemsets(votes,minSup);
fprintf('Minimum Support : %.2f\n', minSup)
fprintf('Frequent Itemsets Found: %d\n', sum(arrayfun(@(x) length(x.freqSets), F)))
fprintf('Max Level Reached : %d-itemsets\n', length(F))
fprintf('Number of Support Data : %d\n', length(S))

Minimum Support : 0.30
Frequent Itemsets Found: 1026
Max Level Reached : 7-itemsets
Number of Support Data : 2530

Finally, let's compare our results with the results from the textbook.

testAntes = [7,21,27;5,11,23;6,15,16;22,31,32];
testConseqs = [2,1,2,1];
testConf = [0.91,0.975,0.935,1];
% Get the rules with 1-item conseq with party classification
dem = rules(arrayfun(@(x) isequal(x.Conseq,1),rules));
rep = rules(arrayfun(@(x) isequal(x.Conseq,2),rules));
% Compare the resultsfor i = 1:size(testAntes,1)
rec = dem(arrayfun(@(x) isequal(x.Ante,testAntes(i,:)),dem));
rec = [rec,rep(arrayfun(@(x) isequal(x.Ante,testAntes(i,:)),rep))];
disp(['{', sprintf('%s, %s, %s',vars{rec.Ante}),'} => ',...
sprintf('{%s}', vars{rec.Conseq})])
fprintf(' Conf: %.2f Lift: %.2f Sup: %.2f\n',rec.Conf,rec.Lift,rec.Sup)
if isequal(rec.Conseq,testConseqs(i));
fprintf(' Correct! Expected Conf %.2f\n\n',testConf(i))
endend

{El Salvador = Yes, Budget Resolution = No, Mx Missile = No} => {Republican}
Conf: 0.91 Lift: 2.36 Sup: 0.30
Correct! Expected Conf 0.91
{Budget Resolution = Yes, Mx Missile = Yes, El Salvador = No} => {Democrat}
Conf: 0.97 Lift: 1.59 Sup: 0.36
Correct! Expected Conf 0.97
{Physician Fee Freeze = Yes, Right To Sue = Yes, Crime = Yes} => {Republican}
Conf: 0.94 Lift: 2.42 Sup: 0.30
Correct! Expected Conf 0.94
{Physician Fee Freeze = No, Right To Sue = No, Crime = No} => {Democrat}
Conf: 1.00 Lift: 1.63 Sup: 0.31
Correct! Expected Conf 1.00

Visualization of Metrics

MATLAB has strong data visualization capabilities. Let's visualize conf, lift and support.

conf = arrayfun(@(x) x.Conf, rules); % get conf as a vector
lift = arrayfun(@(x) x.Lift, rules); % get lift as a vector
sup = arrayfun(@(x) x.Sup, rules); % get support as a vector
colormap cool
scatter(sup,conf,lift*30, lift, 'filled')
xlabel('Support'); ylabel('Confidence')
t = colorbar('peer',gca);
set(get(t,'ylabel'),'String', 'Lift');
title('1984 Congressional Voting Records')

You see that rules we generated match with those in the textbook. The plot of support, confidence and lift is useful to identify rules that are high support, high confidence (upper right region of the plot) and high lift (more magenta). If you type gname into MATLAB prompt, you can interactively identify the indices of those points, and hit Enter to end the interactive session.

We used a small sample dataset in order to learn the basics of Market Basket Analysis. However, for real applications such as web usage pattern mining, we will typically be dealing with quite a large dataset. In a later post, I would like to figure out how to extend what we learned here with the help of MapReduce, introduced in R2014b. Any thoughts about your creative use of "Market Basket Analysis" with your data? Let us know here.

Appendix: Code used for this post

Here is the code used for this post.

findFreqItemsets() generates frequent itemsets using the Apriori method.

type findFreqItemsets.m

function [F,S,items] = findFreqItemsets(transactions,minSup,oneItemsets)
%FINDFREQITEMSETS generates frequent itemsets using Apriori method
% |transactions| is a nested cell array of transaction data or a table
% with |Value| column that contains such cell array. Each row contains a
% nested cell array of items in a single transaction. If supplied as a
% table, it also needs |Key| column with a cell array of transaction ids.
%
% |minSup| is a scalar that represents the minimum support threshold.
% Itemsets that meet this criteria are 'frequent' itemsets.
%
% |oneItemSets| is an optional table that list all single items that
% appear in |transactions| along with their frequencies. Items are listed
% in |Key| column and corresponding counts in |Value| column.
%
% |items| is a cell array of unique items.
%
% |F| is a structure array of frequent itemsets that meet that
% criteria. Items are represented as indices of cell arrat |items|.
%
% |S| is a Map object that maps itemsets to their support values. Items
% are represented as indices of cell arrat |items|.
%
% To learn more about the underlying alogrithm itself, please consult
% with Ch6 of http://www-users.cs.umn.edu/~kumar/dmbook/index.php
% check the number of input arguments
narginchk(2, 3)
if iscell(transactions)
transactions = table(num2cell(1:length(transactions))',transactions,'VariableNames',{'Key','Value'});
end
if nargin == 2
items = transactions.Value;
[uniqItems,~,idx] = unique([items{:}]');
oneItemsets = table(uniqItems,num2cell(accumarray(idx,1)),'VariableNames',{'Key','Value'});
end
% get the total number of transactions
N = height(transactions);
% sort the tables
transactions = sortrows(transactions,'Key');
oneItemsets = sortrows(oneItemsets,'Key');
% get all frequent 1-itemsets
[F,S,items] = getFreqOneItemsets(oneItemsets,N,minSup);
if isempty(F.freqSets)
fprintf('No frequent itemset found at minSup = %.2f\n',minSup)
return
end
% get all frequent k-itemsets where k >= 2
k = 2;
while true
% generate candidate itemsets
Ck = aprioriGen(F(k-1).freqSets, k);
% prune candidates below minimum support threshold
[Fk, support] = pruneCandidates(transactions,Ck,N,items,minSup);
% update support data; if empty, exit the loop
if ~isempty(support)
% create a map object to store suppor data
mapS = containers.Map();
% convert vectors to chars for use as keys
for i = 1:length(support)
mapS(num2str(Ck(i,:))) = support(i);
end
% update S
S = [S; mapS];
else
break;
end
% store the frequent itemsets above minSup
% if empty, exit the loop
if ~isempty(Fk)
F(k).freqSets = Fk;
k = k + 1;
else
break;
end
end
function [F1,S,C1]= getFreqOneItemsets(T,N,minSup)
% GETFREQ1ITEMSETS geneates all frequent 1-itemsets
% 1-items are generated from 1-itemset table |T| and pruned with the
% minimum support threshold |minSup|.
% get 1-itemset candidates
C1 = T.Key;
% get their count
count = cell2mat(T.Value);
% calculate support for all candidates
sup = count ./N;
% create a map object and store support data
S = containers.Map();
for j = 1:length(C1)
S(num2str(j)) = sup(j);
end
% prune candidates by minSup
freqSet = find(sup >= minSup);
% store result in a structure array
F1 = struct('freqSets',freqSet);
end
function [Fk,support] = pruneCandidates(T,Ck,N,items,minSup)
%PRUNECANDIDATES returns frequent k-itemsets
% Compute support for each candidndate in |Ck| by scanning
% transactions table |T| to identify itemsets that clear |minSup|
% threshold
% calculate support count for all candidates
support = zeros(size(Ck,1),1);
% for each transaction
for l = 1:N
% get the item idices
t = find(ismember(items, T.Value{l}));
% increment the support count
support(all(ismember(Ck,t),2)) = support(all(ismember(Ck,t),2)) + 1;
end
% calculate support
support = support./N;
% return the candidates that meet the criteria
Fk = Ck(support >= minSup,:);
end
end

aprioriGen() is a helper function to findFreqItemsets() and generates candidate k-itemsets using the Apriori algorithm.

type aprioriGen.m

function Ck = aprioriGen(freqSets, k)
% APRIORIGEN generates candidate k-itemsets using Apriori algorithm
% This function implements F_k-1 x F_k-1 method, which merges two pairs
% (k-1)-itemsets to generate new k-itemsets if the first (k-2) items are
% identical between the pair.
%
% To learn more about the underlying alogrithm itself, please consult
% with Ch6 of http://www-users.cs.umn.edu/~kumar/dmbook/index.php
% generate candidate 2-itemsets
if k == 2
Ck = combnk(freqSets,2);
else
% generate candidate k-itemsets (k > 2)
Ck = [];
numSets = size(freqSets,1);
% generate two pairs of frequent itemsets to merge
for i = 1:numSets
for j = i+1:numSets
% compare the first to k-2 items
pair1 = sort(freqSets(i,1:k-2));
pair2 = sort(freqSets(j,1:k-2));
% if they are the same, merge
if isequal(pair1,pair2)
Ck = [Ck; union(freqSets(i,:),freqSets(j,:))];
end
end
end
end
end

generateRules() returns the association rules that meets the minium confidence criteria. It also uses aprioriGen() to generate rule candidates.

type generateRules.m

function rules = generateRules(F,S,minConf)
%GENERATERULES returns the association rules found from the frequent
%itemsets based on the minimum confidence threshold |minConf|.
%Association rules are expressed as {ante} => {conseq}.
% |F| is a structure array of frequent itemsets
% |S| is a Map object that maps itemsets to their support values
% |rules| is a structure array of association rules that meet |minConf|
% criteria, along with confidence, lift and support metrics.
%
% To learn more about the underlying alogrithm itself, please consult
% with Ch6 of http://www-users.cs.umn.edu/~kumar/dmbook/index.php
rules = struct('Ante',{},'Conseq',{},'Conf',{},'Lift',{},'Sup',{});
% iterate over itemset levels |k| where k >= 2
for k = 2:length(F)
% iterate over k-itemsets
for n = 2:size(F(k).freqSets,1)
% get one k-itemset
freqSet = F(k).freqSets(n,:);
% get 1-item consequents from the k-itemset
H1 = freqSet';
% if the itemset contains more than 3 items
if k > 2
% go to ap_genrules()
rules = ap_genrules(freqSet,H1,S,rules,minConf);
else
% go to pruneRules()
[~,rules] = pruneRules(freqSet,H1,S,rules,minConf);
end
end
end
function rules = ap_genrules(Fk,H,S,rules,minConf)
%AP_GENRULES generate candidate rules from rule consequent
% |Fk| is a row vector representing a frequent itemset
% |H| is a matrix that contains a rule consequents per row
% |S| is a Map object that stores support values
% |rules| is a structure array that stores the rules
% |minConf| is the threshold to prune the rules
m = size(H,2); % size of rule consequent
% if frequent itemset is longer than consequent by more than 1
if length(Fk) > m+1
% prune 1-item consequents by |minConf|
if m == 1
[~,rules] = pruneRules(Fk,H,S,rules,minConf);
end
% use aprioriGen to generate longer consequents
Hm1 = aprioriGen(H,m+1);
% prune consequents by |minConf|
[Hm1,rules] = pruneRules(Fk,Hm1,S,rules,minConf);
% if we have consequents that meet the criteria, recurse until
% we run out of such candidates
if ~isempty(Hm1)
rules = ap_genrules(Fk,Hm1,S,rules,minConf);
end
end
end
function [prunedH,rules] = pruneRules(Fk,H,S,rules,minConf)
%PRUNERULES calculates confidence of given rules and drops rule
%candiates that don't meet the |minConf| threshold
% |Fk| is a row vector representing a frequent itemset
% H| is a matrix that contains a rule consequents per row
% |S| is a Map object that stores support values
% |rules| is a structure array that stores the rules
% |minConf| is the threshold to prune the rules
% |prunedH| is a matrix of consequents that met |minConf|
% initialize a return variable
prunedH = [];
% iterate over the rows of H
for i = 1:size(H,1);
% a row in H is a conseq
conseq = H(i,:);
% ante = Fk - conseq
ante = setdiff(Fk, conseq);
% retrieve support for Fk
supFk =S(num2str(Fk));
% retrieve support for ante
supAnte =S(num2str(ante));
% retrieve support for conseq
supConseq =S(num2str(conseq));
% calculate confidence
conf = supFk / supAnte;
% calculate lift
lift = supFk/(supAnte*supConseq);
% if the threshold is met
if conf >= minConf
% append the conseq to prunedH
prunedH = [prunedH, conseq];
% generate a rule
rule = struct('Ante',ante,'Conseq',conseq,...
'Conf',conf,'Lift',lift,'Sup',supFk);
% append the rule
if isempty(rules)
rules = rule;
else
rules = [rules, rule];
end
end
end
end
end

getData is a script that fetches the congressional voting data using the RESTful API access function webread() that was introduced in R2014b, and preprocesses the data for Market Basket Analysis.

type getData.m

% Let's load the dataset
clearvars; close all; clc
% the URL of the source data
url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/voting-records/house-votes-84.data';
% define the function handle to process the data you read from the web
myreadtable = @(filename)readtable(filename,'FileType','text','Delimiter',',','ReadVariableNames',false);
% set up the function handle as |ContentReader| for |webread| function
options = weboptions('ContentReader',myreadtable);
% get the data
votes = webread(url,options);
% get the variable names
votes.Properties.VariableNames = {'party',...
'handicapped_infants',...
'water_project_cost_sharing',...
'adoption_of_the_budget_resolution',...
'physician_fee_freeze',...
'el_salvador_aid',...
'religious_groups_in_schools',...
'anti_satellite_test_ban',...
'aid_to_nicaraguan_contras',...
'mx_missile',...
'immigration',...
'synfuels_corporation_cutback',...
'education_spending',...
'superfund_right_to_sue',...
'crime',...
'duty_free_exports',...
'export_administration_act_south_africa'};
% display the small portion of the data
disp(votes(1:5,1:4))
% turn table into matrix
C = table2array(votes);
% create a logical array of yes votes
arr = strcmp(C(:,2:end), 'y');
% append logical array of no votes
arr = [arr, strcmp(C(:,2:end), 'n')];
% create a logical arrays of party affiliation
dem = strcmp(C(:,1),'democrat');
rep = strcmp(C(:,1),'republican');
% combine them all
arr = [dem,rep,arr];
% create cell to hold the votes of each member of congress
votes = cell(size(dem));
% for each member, find the indices of the votes
for i = 1:length(dem)
votes{i} = find(arr(i,:));
end
% variable names that maps to the indices
vars = {'Democrat',...
'Republican',...
'Handicapped Infants = Yes',...
'Water Project = Yes',...
'Budget Resolution = Yes',...
'Physician Fee Freeze = Yes',...
'El Salvador = Yes',...
'Religious Groups = Yes',...
'Anti Satellite = Yes',...
'Nicaraguan Contras = Yes',...
'Mx Missile = Yes',...
'Immigration = Yes',...
'Synfuels Corp = Yes',...
'Education Spending = Yes',...
'Right To Sue = Yes',...
'Crime = Yes',...
'Duty Free = Yes',...
'South Africa = Yes',...
'Handicapped Infants = No'...,
'Water Project = No',...
'Budget Resolution = No',...
'Physician Fee Freeze = No',...
'El Salvador = No',...
'Religious Groups = No',...
'Anti Satellite = No',...
'Nicaraguan Contras = No',...
'Mx Missile = No',...
'Immigration = No',...
'Synfuels Corp = No',...
'Education Spending = No',...
'Right To Sue = No',...
'Crime = No',...
'Duty Free = No',...
'South Africa = No'};

]]>Big DataSimuflate: A time-based simulation of the Patriots’ footballs
http://feedproxy.google.com/~r/SethOnSimulink/~3/kAG3q98ufA4/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/seth/2015Q1/footbal_results.png"/></div><p>Today I want to share a discussion I had with a few fellow guest bloggers: Corey Lagunowich, Zack Peters, Matt Brauer and Seth Popinchalk.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/01/27/simuflate/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4397Wed, 28 Jan 2015 04:00:19 +0000Today I want to share a discussion I had with a few fellow guest bloggers: Corey Lagunowich, Zack Peters, Matt Brauer and Seth Popinchalk.

Deflate-gate

Corey: Hey guys, did you hear about the Patriots "Deflate-gate" story?

Zack: Yeah, they say that 11 out of 12 balls were inflated 2 PSI less than the minimum value!

Corey: There could be many reasons why the balls were measured at the prescribed 12.5 PSI before the game, but lose pressure by halftime.

Matt: Do you think they cheated, or that this situation can be explained by some math and physics?

Seth: We should make a Simulink model of the thermodynamics involved in that story.

Guy: My thermodynamics knowledge is pretty rusty, but I bet that with Simscape we can build up something in just a few minutes.

A few minutes later...

Guy: Et voilà! Here is what the model looks like:

Seth: Good job guys! Can you explain to us how you built that?

Zack: Sure! As you can see, the football is modeled using a Constant Volume Pneumatic Chamber. To the thermal port, we connected in series the convection inside, the conduction through the leather and the convection outside.

Guy: And the Compressor Control subsystem is there to first inflate the ball with air at a fixed temperature up to a desired pressure. A Stateflow chart turns off the inflating action for the rest of the simulation once the ball first gets up to the prescribed pressure.

Matt: Great, I heard a few theories that the ball being inflated with hot air could explain this pressure drop. Can we confirm that?

The Results

Corey: With this model, we are able to follow the evolution of the temperature and pressure inside the ball throughout the game day.

Zack: The story begins inside, with an ambient temperature of 73 degrees Fahrenheit. To test the hot air hypothesis, the ball is inflated with air at 150 degrees Fahrenheit. As you can see, this takes a bit less than a minute, and almost immediately after stopping to feed hot air, the air inside the ball cools down to ambient temperature.

Matt: I guess that rules out the hot air theory. I am curious though, I am sure I read somewhere that the football should take longer than that to cool down when taken from inside to outside.

Guy: We are getting there... One interesting thing I found is that the component with the most impact here is the thermal mass of the leather. This is why the hot air is immediately cooled down by the leather at ambient temperature.

Corey: The simulation times all correspond to what happened, or at least what's been reported. The balls are filled roughly 2 hours before game time, and not taken outside until just minutes before kickoff. As the simulation goes, we then move the ball outside, to about 50 degrees Fahrenheit, and then it takes roughly 20 minutes for the air inside to stabilize at the new ambient temperature.

Zack: The ideal gas law governs the magnitude of the pressure drop - roughly 1 PSI when the air temperature changes from 73 to 50 degrees Fahrenheit. But thanks to Simscape we can see the time it takes to get there.

Seth: And I can see that at halftime you simulate the ball being taken back inside and warming up.

Matt: We simulated the ambient temperature going back up to room temperature, assuming that the refs took the balls into their locker room to test them. If they stuck a gauge into the footballs right away, the pressure would still have read low. But wait a few more minutes…

Corey: When combined with other possible effects like accuracy of the gauge and leakage, it seems pretty obvious that the Patriots didn't cheat.

Guy: Corey, you are too much of a fan... The results of our simulation show how the ball pressure and temperature varied through that day, but in no way can they confirm what happened that day... too many unknowns.

Now it's your turn

If you are interested, you can download our simulation here.

Do you also use Simulink or Simscape to analyze controversial topics? Let us know by leaving a comment here.

And one last thing ... if you want to know just which way we might be biased on this issue, well, just check out where our corporate headquarters are located.

]]>Visualizing retention data over time in MATLAB
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/v77gJX4MN2g/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2015/01/26/visualizing-retention-data-over-time-in-matlab/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201501/631/62009828001_4021096362001_357-Retention.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 03:40</span>
</div>
</a></div><p>I have some data showing how many people come into and leave a group over time. I wanted to see the retention of members from one month to the next. Using area plots I was able to make a “stacked area plot” to get the visualization I wanted. This... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2015/01/26/visualizing-retention-data-over-time-in-matlab/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1564Mon, 26 Jan 2015 14:55:58 +0000

]]>Fourier series animation using phasor addition
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/2OhlQGxGyys/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/steve/files/harmonic-circles-screenshot.png"/></div><p>
I just saw a great animation illustrating the Fourier series decomposition of a square wave. Check it out. The video includes two different animations, so be sure to watch it all the way through to see the second one.
... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2015/01/26/fourier-series-animation-using-phasor-addition/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1268Mon, 26 Jan 2015 12:00:57 +0000
I just saw a great animation illustrating the Fourier series decomposition of a square wave. Check it out. The video includes two different animations, so be sure to watch it all the way through to see the second one.

]]>Fourier transformsWhat was Sheraton Trying to Tell Us?
http://feedproxy.google.com/~r/mathworks/loren/~3/gNLgiqXfTBk/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/Sher.png"/></div><p>Today's post is by Sean de Wolski and me. We spent a few days this week at training for the start of the year, held in a Sheraton hotel.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/01/23/what-was-sheraton-trying-to-tell-us/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1096Fri, 23 Jan 2015 20:31:29 +0000

What was Sheraton Trying to Tell Us?

Today's post is by Sean de Wolski and me. We spent a few days this week at training for the start of the year, held in a Sheraton hotel.

At the completion of our stay, we each got a note in our rooms. Here is a sample.

Next guess was how the text was encoded - 4 bits, 8 bits? We figured 8 was most likely.

newstr = reshape(str',8,[])';

Let's look at the first few lines.

snip = newstr(1:3,:)

snip =
01010100
01101000
01100001

Could these be text?

bin2dec(snip)

ans =
84
104
97

These could be ASCII characters!

Have you ever misread an O for a 0? Or an l for a 1? We have! So we decided to program defensively for those and any other weird swaps that might confuse the OCR.

]]>Introducing Community Profile Pages
http://feedproxy.google.com/~r/mathworks/desktop/~3/dt6XYgnZjOg/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/community/files/community-profile02.png"/></div><p>This week, guest-blogger David Wey welcomes the new community profile pages. David is a Senior Developer for MathWorks community sites.
Community Improvements: New Profile Pages
by David Wey
Have you ever wondered who created the file you’re looking at on the File Exchange, or who that nice person was that answered your... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/community/2015/01/23/introducing-community-profile-pages/">read more >></a></p>http://blogs.mathworks.com/community/?p=3150Fri, 23 Jan 2015 13:30:36 +0000This week, guest-blogger David Wey welcomes the new community profile pages. David is a Senior Developer for MathWorks community sites.

Community Improvements: New Profile Pages

by David Wey

Have you ever wondered who created the file you’re looking at on the File Exchange, or who that nice person was that answered your question on MATLAB Answers? Maybe you thought to yourself I wonder what else that person has done? As of this week answering these questions got a lot easier.

NEW PROFILE PAGES

Our latest update to MATLAB Central includes new profile pages. These pages aggregate content from across all the MATLAB Central areas, thus easing the pain of having separate profile pages for each area (Answers, Cody, File Exchange, and so on). On these profile pages you’ll be able to see how active the person has been and how recently they’ve participated in the community. You will see all their answers, files, Cody problems, trends, and links.

To view someone’s new profile page, click on any author link wherever you see it, whether in Cody or the File Exchange. The popup displays a summary of information related to that area of MATLAB Central (such as the File Exchange or Answers). Click again on the user name to view the profile page or click on the other labels to see their submissions as a filtered search.

We hope you’ll enjoy the new profile pages and we welcome your thoughts and reactions. Check out some of your favorite community member pages, as well as your own, and let us know what you think in the comments below.

]]>UncategorizedDivergent colormaps
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/s-GhhqD554Y/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/steve/2015/divergent_colormap_10.jpg"/></div><p>I promised earlier to post about divergent colormaps. A divergent colormap is usually constructed by concatenating two colormaps together that have different color schemes.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2015/01/20/divergent-colormaps/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1265Tue, 20 Jan 2015 23:35:44 +0000

I promised earlier to post about divergent colormaps. A divergent colormap is usually constructed by concatenating two colormaps together that have different color schemes.

Here is an example of a divergent colormap from www.ColorBrewer.org by Cynthia A. Brewer, Geography, Pennsylvania State University.

This colormap has the brightest color in this middle, which is typical for divergent colormaps.

A divergent colormap is used to compare data values to a reference value in a way that visually highlights whether values are above or below the reference.

Let me show you an example using sea surface temperatures. I learned today that 28 degrees Celsius has a specific meaning in hurricane modeling and prediction. It is considered to be the minimum temperature for hurricane formation.

Using the NASA EOSDIS Reverb tool, I found and downloaded sea surface temperature from the Aqua satellite for August 1, 2011. The data arrived as a netCDF file.

Now I want to make use of our divergent colormap to highlight the temperature of interest, 28 degrees. To do that, I'll use the CLim property of the axes object. The CLim property is a two-element vector. The first element tells us the data value that maps to the lowest colormap color, and the second element tells us the data value that maps to the highest colormap color. MATLAB computes these values automatically based on the range of the data.

get(gca,'CLim')

ans =
-1.9500 33.6000

To get the middle colormap color to represent 28 degrees, we need to pick our two CLim values so that they are equally far away from 28 degrees. After some experimentation, I settled on a range from 22 degrees to 34 degrees.

set(gca,'CLim',[22 34])

That's looking better, but it isn't helpful to have the missing data displayed using the bottom color of the colormap. I will fix that problem using this procedure.

Construct a gray mask image (in truecolor format).

Display it on top of our sea surface temperature image.

Set the AlphaData of the gray mask image so that the sea surface temperatures show through.

Here are just the first two steps.

mask = isnan(ssta);
graymask = 0.5 * mask;
graymask = repmat(graymask,1,1,3);
hold on
h = imshow(graymask);
hold off

The missing values are displayed as gray, but we can't see the sea surface temperatures. Fix that by using the NaN-mask as the overlay image's AlphaData.

h.AlphaData = mask; % This syntax requires MATLAB R2014b.

Temperatures close to the threshold are displayed using yellow, the brightest color in the colormap and the one right in the middle. Temperatures above the threshold are displayed using oranges and reds, while temperatures below are displayed using greens and blues.

Before I go, I want to give a shout to the new Developer Zone blog that just started today. If you are into advanced software developer topics, I encourage you to head over there and take a look.

]]>ColormapThe Three n Plus One Conjecture
http://feedproxy.google.com/~r/mathworks/moler/~3/yp1Y59gvRc0/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/cleve/threenplus1_7.png"/></div><p>If $n$ is odd, replace $n$ by $3n+1$, if not, replace $n$ by $n/2$. Repeat. A famous conjecture made by Lothar Collatz is that no matter what value of $n$ is chosen to start, the process eventually terminates at $n=1$. Do not expect a proof, or a counterexample, in this blog.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2015/01/19/the-three-n-plus-one-conjecture/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1157Mon, 19 Jan 2015 17:00:34 +0000

If $n$ is odd, replace $n$ by $3n+1$, if not, replace $n$ by $n/2$. Repeat. A famous conjecture made by Lothar Collatz is that no matter what value of $n$ is chosen to start, the process eventually terminates at $n=1$. Do not expect a proof, or a counterexample, in this blog.

You might want to download function threenplus1.m from my textbook Numerical Computing with MATLAB. This not only produces the graphs in this post, it also has a pair of uicontrols that allow you to increase and decrease the value of n.

Here is the core of threenplus1. The generated sequence is collected in a vector y. We don't know ahead of time how long y is going to be. In fact, that's the point of the computation. Does the while loop terminate and, if so, how long is y?. The ability of MATLAB to grow a vector an element at a time comes in very handy here.

dbtype 44:52threenplus1

44 y = n;
45 while n > 1
46 if rem(n,2)==0
47 n = n/2;
48 else
49 n = 3*n+1;
50 end
51 y = [y n];
52 end

n = 7

A good example is provided by starting with $n = 7$. Here is the sequence collected in y. Its length is 17 and its max is 52.

7 22 11 34 17 52 26 13 40 20 10 5 16 8 4 2 1

Here is the plot. The y-axis is logarithmic in anticipation of some fairly large values in later runs. Notice that as soon as the sequence hits a power of 2, it plunges to termination.

n = 108, 109, 110

It is interesting to see how the graphs produced by threenplus1 change from one value of $n$ to the next. Here are the first ten elements of y for $n$ = 108, 109, and 110.

After the first eight steps all three sequences reach the same value, so they are identical after that. The length is 114 and the max is 9232. Here are the three graphs, superimposed on one another. The similarity between the three sequences is very apparent when you use the increment n buttons provided by threenplus1.

n = 2^13-1

Prime numbers of the form $n = 2^p-1$ where $p$ itself is prime are known as Mersenne primes. But that's a whole 'nother story. For now, suffice it to say that such numbers are candidates to produce 3n+1 sequences with a high maximum value. Here is $n = 2^{13}-1$. The sequence reaches 6,810,136 before it eventually terminates in 159 steps.

I have to admit that this triggers a bug in some versions of MATLAB. If you see specious characters in the formatting of the legend on the y-axis, insert this line near the end of threenplus1.

set(gca,'yticklabel',int2str(get(gca,'ytick')'))

Orbits

Check out the two references below about orbits associated with the Collatz conjecture.

Background

When I Google "three n plus 1", the first link is to an excellent Wikipedia article on the "Collatz conjecture". According to Wikipedia, the famous German mathematican Lothar Collatz first made the conjecture, in 1937, that the process terminates for any starting value. I am already familiar with Collatz's name because of his work in numerical analysis and approximation theory, including finite difference methods for differential operators. My mentor John Todd had some contact with Collatz shortly after the end of World War II that could be the subject of a blog post some day.

Many others have also been interested in the 3n+1 problem. Stanislaw Ulam (see my previous blog post), Shiszuo Kakutani, Sir Bryn Thwaites, Helmut Hasse, Paul Erdos, John Horton Conway, Douglas Hofstadter, and Richard Guy. But Wikipedia quotes Paul Erdos as saying about the Collatz conjecture: "Mathematics may not be ready for its solution."

]]>Robot Game-Playing in MATLAB – Part 2
http://feedproxy.google.com/~r/mathworks/desktop/~3/Vur9W6mQ-Zs/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/community/files/reversi.gif"/></div><p>Last week we showed how to use a Monte Carlo approach to write a program capable of playing Tic Tac Toe. But that was just our warm-up exercise. Now let’s build another game that can take advantage of our game-playing harness.
The first game I thought of was Connect Four,... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/community/2015/01/19/robot-game-playing-in-matlab-part-2/">read more >></a></p>http://blogs.mathworks.com/community/?p=3124Mon, 19 Jan 2015 13:30:44 +0000Last week we showed how to use a Monte Carlo approach to write a program capable of playing Tic Tac Toe. But that was just our warm-up exercise. Now let’s build another game that can take advantage of our game-playing harness.

The first game I thought of was Connect Four, since it’s a straightforward extension of Tic Tac Toe. And sure enough, it’s not to hard to convert our Tic Tac Toe engine into a Connect Four engine.

CONNECT FOUR

In Connect Four, your board is a 6-by-7 grid, and your goal is to be the first to put four pieces in a line. The game board is oriented vertically, so gravity causes your pieces to drop to the bottom of the grid. This means that you have at most seven legal moves. This is good for us, since it constrains the number of possible moves, which in turn makes it quicker for us to play hundreds of games.

The game harness botMoves doesn’t change at all. We just make a new class, ConnectFour, that knows the rules and can draw the Connect Four game board. And when I say “knows the rules”, I mean exactly two things.

What legal moves can I make right now?

If the game is over, who won?

That’s it.

Finding legal moves is easy. We just look for any column that has zeroes left in it. But how do you determine if someone has won the game? We need to scan the board for four-in-a-row of either color in any of four directions: vertical, horizontal, and two diagonals. This seemed like it was going to call for some careful thinking, until I remembered that there was a Cody problem about this very topic! See Problem 90. Connect Four Win Checker.

I wired everything together, and it plays a pretty mean game. Here’s an animation of my Connect Four bot playing itself.

CONNECT FOUR VARIANTS

The Monte Carlo approach means that we haven’t had to put a lot of special game-specific smarts into the code. This gives us great freedom to explore variants on our basic Connect Four game. Suppose, instead of racing to get four in a row, players were trying to be the first to make an “L” shape. In order to make this work, I just had to change a few lines in the isGameOver method. And rather than clone and modify the entire ConnectFour class, I built a FourGameBase class from which I could quickly subclass variants. This base class knows how to draw the board and pick legal moves. The winning conditions are handled by the subclasses.

Here’s my ConnectEl game bot playing itself.

We can imagine a whole family of Connect Tetris variants: be the first to make a square, or a T-shape, or an S-shape. Not all of these variants are interesting as games. I built square and T-shape versions of the game, and in both cases defensive play was so easy that, as with Tic Tac Toe, all the games ended in ties. But this ability to rapidly create and play games leads to another observation. We can use our Monte Carlo code to mutate games in search of interesting variants. An interesting game is one in which either side has a reasonable chance of winning. An uninteresting game consistently ends in a tie or victory for a particular side. So we can effectively use bots to do our play-testing for us. The optimal strategies simply emerge from the Monte Carlo soup.

I want to mention one last Connect Four variant before I move on: Four Corners. In Four Corners, you try to force your opponent to be the first to put his pieces into the four corners of a rectangle.

In this case, I used code from a Cody answer by Alfonso Nieto-Castanon to test for the victory condition: Spot the rectangle.

Here is my bot slugging it out with itself in Four Corners.

This is a fun game to play against the bot. It’s surprisingly tricky to spot all the potential rectangles on the board.

REVERSI

For my last example, I wrote a program that does a passable job of playing Reversi. As before, I used Cody to diminish the coding effort by posting two problems.

If you, like me, are intrigued with the possibilities of Monte Carlo gaming, I encourage you to take a look at my code on GitHub. Maybe you can add some new games to the list that I’ve already created. There are plenty of games that would work with this harness. Or you can improve the code that’s already there, making it speedier or more efficient. Because of the brute force nature of the Monte Carlo technique, it would be an ideal place to experiment with parallel programming. Or you might make it more interactive, with a click-to-move interface.

I’ll stop here, but this is clearly a rich space to explore. It’s a powerful reminder that computers are very different from brains. Sometimes you can make up for not being clever by being stupid very fast.

]]>UncategorizedControlling the Location of the Generated Code and Temporary files
http://feedproxy.google.com/~r/SethOnSimulink/~3/aWVkk1MrjCs/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/seth/2014Q2/simulinkPrefs.png"/></div><p>Simulink often needs to generate files to simulate a model. An example of this is the model reference simulation target. Today I will describe a few option to control where those files are created.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/01/16/controlling-the-location-of-the-generated-code-and-temporary-files/">read more >></a></p>http://blogs.mathworks.com/seth/?p=3863Fri, 16 Jan 2015 19:37:46 +0000Simulink often needs to generate files to simulate a model. An example of this is the model reference simulation target. Today I will describe a few option to control where those files are created.

Simulation Target

If you click play in the following model:

You will notice that a folder named slprj with a lot of code and 4 MEX-files will be created:

Now, let's say I need to often switch between different projects, and each of those projects need to use different folders for code generation. In that case, it is possible to use Simulink.fileGenControl. For example, if I include Simulink Project in my workflow, the startup shortcut for my project could look like:

That way, all the code and MEX-files generated by the model will go in those folders and the current directory will remain clean.

Once you are done working on this project, you can restore the values stored in the Simulink Preferences using this line:

Simulink.fileGenControl('reset')

Don't Clean Up: Faster Initialization

Simulink checks the cached files at initialization to make sure they are up to date and match the model. This prevents re-generating files, and results in faster initialization. If you regularly work with the same models, keeping these derived files around can help you save time when switching between projects.

Now it's your turn

Let us know how you use this feature by leaving a comment here

]]>Process “Big Data” in MATLAB with mapreduce
http://feedproxy.google.com/~r/mathworks/loren/~3/hclgcbIrwtc/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/airlineExample_01.png"/></div><p>Today I’d like to introduce guest blogger Ken Atwell who works for the MATLAB Development team here at MathWorks. Today, Ken will be discussing with you the <a rel="nofollow" target="_blank" href="http://www.mathworks.com/help/matlab/mapreduce.html">MapReduce</a> programming technique now available in the R2014b release of MATLAB. MapReduce provides a way to process large amounts of file-based data on a single computer in MATLAB. For very large data sets, the same MATLAB code written using MapReduce can also be run on the "big data" platform, Hadoop®.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/01/15/process-big-data-in-matlab-with-mapreduce/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1085Thu, 15 Jan 2015 13:57:34 +0000

Today I’d like to introduce guest blogger Ken Atwell who works for the MATLAB Development team here at MathWorks. Today, Ken will be discussing with you the MapReduce programming technique now available in the R2014b release of MATLAB. MapReduce provides a way to process large amounts of file-based data on a single computer in MATLAB. For very large data sets, the same MATLAB code written using MapReduce can also be run on the "big data" platform, Hadoop®.

The data set we will be using consists of records containing flight performance metrics for USA domestic airlines for the years 1987 through 2008. Each year consists of a separate file. If you have experimented with "big data" before, you may already be familiar with this data set. The full data set can be downloaded from here. A small subset of the data set, airlinesmall.csv, is also included with MATLAB® to allow you to run this and other examples without downloading the entire data set.

Introduction to mapreduce

MapReduce is a programming technique that is used to "divide and conquer" big data problems. In MATLAB, the mapreduce function requires three input arguments:

A datastore for reading data into the "map" function in a chunk-wise fashion.

A "map" function that operates on the individual chunks of data. The output of the map function is a partial calculation. mapreduce calls the map function one time for each chunk of data in the datastore, with each call operating independently from other map calls.

A "reduce" function that is given the aggregate outputs from the map function. The reduce function finishes the computation begun by the map function, and outputs the final answer.

This is an over-simplification to some extent, since the output of a call to the map function can be shuffled and combined in interesting ways before being passed to the reduce function. This will be examined later.

Use mapreduce to perform a computation

Let's look at a straightforward example to illustrate how mapreduce is used. In this example we want to find the longest flight time out of all the flights recorded in the entire airline data set. To do this we will:

Create a datastore object to reference our airline data set

Create a "map" function that computes the maximum flight time in each chunk of data in the datastore.

Create a "reduce" function that computes the maximum value among all of the maxima computed by the calls to the map function.

Create a datastore

datastore is used to access collections of tabular text files stored on a local disk, or the Hadoop® Distributed File System (HDFS™). It is also the mechanism for providing data in a chunk-wise manner to the map calls when using mapreduce. A previous blog post explained how datastore works and how it is used for reading collections of data that are too large to fit in memory.

Let's create our datastore and first preview the data set. This allows us to peek into the data set, and identify the format of the data and the columns that contain the data we are interested in. The preview will normally provide a small chunk of data that contains all the columns present in the data set, although I will only show a handful of columns to make our article more readable.

For the computation we are looking to perform, we only need to look at the "ActualElapsedTime" column: it contains the actual flight time. Let's configure our datastore to only provide this one column to our map function.

ds.SelectedVariableNames = {'ActualElapsedTime'};

Create a map function

Now we will write our map function (MaxTimeMapper.m). Three input arguments will be provided to our map function by mapreduce:

The input data, "ActualElapsedTime", which is provided as a MATLAB table by datastore.

A collection of configuration and contextual information, info. This can be ignored in most cases, as it is here.

An intermediate data storage object, where the results of the calculations from the map function are stored. Use the add function to add key/value pairs to this intermediate output. In this example, we have arbitrarily chosen the name of the key ('MaxElapsedTime').

Our map function finds the maximum value in the table 'data', and saves a single key ('MaxElapsedTime') and corresponding value to the intermediate data storage object. We will now go ahead and save the following map function (MaxTimeMapper.m) to our current folder.

function MaxTimeMapper(data, ~, intermediateValuesOut)
maxTime = max(data{:,:});
add(intermediateValuesOut, 'MaxElapsedTime', maxTime);
end

Create a reduce function

Next, we create the reduce function (MaxTimeReducer.m). Three input arguments will also be provided to our reduce function by mapreduce:

A set of input "keys". Keys will be discussed further below, but they can be ignored in some simple problems, as they are here.

An intermediate data input object that mapreduce passes to the reduce function. This data is the output of the map function and is in the form of key/value pairs. We will use the hasnext and getnext functions to iterate through the values.

A final output data storage object where the results of the reduce calculations are stored. Use the add and addmulti functions to add key/value pairs to the output.

Our reduce function will receive a list of values associated with the 'MaxElapsedTime' key which were generated by the calls to our map function. The reduce function will iterate through these values to find the maximum. We create the following reduce function (MaxTimeReducer.m) and save it to our current folder.

function MaxTimeReducer(~, intermediateValuesIn, finalValuesOut)
maxElapsedTime = -inf;
while hasnext(intermediateValuesIn)
maxElapsedTime = max(maxElapsedTime, getnext(intermediateValuesIn));
end
add(finalValuesOut, 'MaxElapsedTime', maxElapsedTime);
end

Run mapreduce

Once the map and reduce functions are written and saved in our current folder, we can call mapreduce and reference the datastore, map function, and reduce function to run our calculations on our data. The readall function is used here to display the results of the MapReduce algorithm.

result = mapreduce(ds, @MaxTimeMapper, @MaxTimeReducer);
readall(result)

Parallel mapreduce execution on the local cluster:
********************************
* MAPREDUCE PROGRESS *
********************************
Map 0% Reduce 0%
Map 100% Reduce 100%
ans =
Key Value
________________ ______
'MaxElapsedTime' [1650]

If Parallel Computing Toolbox is available, MATLAB will automatically start a pool and parallelize execution of the map functions. Since the number of calls to the map function by mapreduce corresponds to the number of chunks in the datastore, parallel execution of map calls may speed up overall execution.

Use of keys in mapreduce

The use of keys is an important feature of mapreduce. Each call to the map function can add intermediate results to one or more named "buckets", called keys.

If the map function adds values to multiple keys, this leads to multiple calls to the reduce function, with each reduce call working on only one key's intermediate values. The mapreduce function automatically manages this data movement between the map and reduce phases of the algorithm.

Calculating group-wise metrics with mapreduce

The behavior of the map function in this case is more complex and will illustrate the benefits to using multiple keys to store your intermediate data. For every airline found in the input data, use the add function to add a vector of values. This vector is a count of the number of flights for that carrier for each day in the 21+ years of data. The carrier code is the key for this vector of values. This ensures that all of the data for each carrier will be grouped together when mapreduce passes it to the reduce function.

Here's our new map function (CountFlightsMapper.m).

for i = 1:numel(airlineName)
dayTotals = accumarray(dayNumber(airlineIndex==i), 1, [daysSinceEpoch, 1]);
add(intermediateValuesOut, airlineName{i}, dayTotals);
endend

The reduce function is less complex. It simply iterates over the intermediate values and adds the vectors together. At completion, it outputs the values in this accumulated vector. Note that the reduce function does not need to sort or examine the intemediateKeysIn values; each call to the reduce function by mapreduce only passes the values for one airline.

Here's our new reduce function (CountFlightsReducer.m).

while hasnext(intermediateValuesIn)
dayArray = dayArray + getnext(intermediateValuesIn);
end
add(finalValuesOut, intermediateKeysIn, dayArray);
end

To run this new analysis, just reset the datastore and select the active variables of interest. In this case we want the date (year, month, day) and the airline carrier name. Once the map and reduce functions are written and saved in our current folder, mapreduce is called referencing the updated datastore, map function, and reduce function.

reset(ds);
ds.SelectedVariableNames = {'Year', 'Month', 'DayofMonth', 'UniqueCarrier'};
result = mapreduce(ds, @CountFlightsMapper, @CountFlightsReducer);
result = result.readall();

Parallel mapreduce execution on the local cluster:
********************************
* MAPREDUCE PROGRESS *
********************************
Map 0% Reduce 0%
Map 100% Reduce 50%
Map 100% Reduce 100%

Up until this point we have been only analyzing the sample data set (airlinesmall.csv). To look at the number of flights each day over the full dataset, let's load the results of running our new MapReduce algorithm on the entire data set.

load airlineResults

Visualizing the results

Before looking at the number of flights per day for the top 7 carriers, let us first apply a filter to the data to smooth out the effects of weekend travel. This would otherwise clutter the visualization.

lines = result.Value;
lines = horzcat(lines{:});
[~,sortOrder] = sort(sum(lines), 'descend');
lines = lines(:,sortOrder(1:7));
result = result(sortOrder(1:7),:);
lines(lines==0) = nan;
for carrier=1:size(lines,2)
lines(:,carrier) = filter(repmat(1/7, [7 1]), 1, lines(:,carrier));
end

figure('Position',[1 1 800 800]);
plot(datetime(1987,10,1):caldays(1):datetime(2008,12,31),lines)
title ('Domestic airline flights per day per carrier')
xlabel('Date')
ylabel('Flights per day (7-day moving average)')
try
carrierList = readtable('carriers.csv.txt', 'ReadRowNames', true);
result.Key = cellfun(@(x) carrierList(x, :).Description, result.Key);
catchend
legend(result.Key, 'Location', 'SouthOutside')

The interesting characteristic highlighted by the plot is the growth of Southwest Airlines (WN) during this time period.

Running mapreduce on Hadoop

The code we have just created was done so on our desktop computer and it allows us to analyze data that normally would not fit into the memory of our machine. But what if our data was stored on the "big data" platform, Hadoop® and the data set is too large to reasonably transfer to our desktop computer?

Our datastore is updated to reference the location of the data files in the Hadoop® Distributed File System (HDFS™).

The mapreducer object is updated to reference Hadoop as the runtime environment.

After these configuration changes, the algorithm can then be run on Hadoop.

MATLAB Compiler can be used to deploy applications created using MapReduce, and it can also generate Hadoop specific components for use in production settings.

Conclusion

Learn more about MapReduce in MATLAB with our algorithm examples and let us know how you might see MapReduce being used to analyze your big data here.

]]>Robot Game-Playing in MATLAB
http://feedproxy.google.com/~r/mathworks/desktop/~3/GMPbDTkAtcM/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/community/files/tictac.gif"/></div><p>A story about just-in-time expertise. Sometimes the best learning is no learning.
COMPUTERS, CHESS, AND GO
I read an article in IEEE Spectrum about computer programs that play Go (AIs Have Mastered Chess. Will Go Be Next?). If you review the history of game-playing computers, you’ll see that chess programs improved steadily... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/community/2015/01/09/robot-game-playing-in-matlab/">read more >></a></p>http://blogs.mathworks.com/community/?p=3109Fri, 09 Jan 2015 21:15:02 +0000A story about just-in-time expertise. Sometimes the best learning is no learning.

COMPUTERS, CHESS, AND GO
I read an article in IEEE Spectrum about computer programs that play Go (AIs Have Mastered Chess. Will Go Be Next?). If you review the history of game-playing computers, you’ll see that chess programs improved steadily until eventually they could beat the best human players. Go programs, on the other hand, have been stuck at a level of play that was nowhere close to the best human. Why is that?

The basic element of a game-playing program is the look-ahead. Essentially, the program says “If I move here, is that better or worse than if I move there?” In chess, this is straightforward to evaluate. But in Go, this basic look-ahead strategy doesn’t work so well. It’s much harder to evaluate whether one board position is stronger than another.

But recently, Go programs have started to get much better. What happened?

TWO IDIOTS FINISH THE GAME
Go programs have improved by applying a Monte Carlo technique. It’s nothing like how a human plays, but it works remarkably well. And it only works because we can ask the computer to do a lot of dumb stuff very quickly. I call it “Two Idiots Finish the Game”.

Consider the following situation. You’ve reached a critical point in the game. We’ll call it position X. You’re considering
move A and move B. Which one should you make? Now instead of looking just one move ahead, play the game all the way to completion. But there’s an obvious problem with this. If you’re not smart enough to figure out your next move, how can you play an entire game? Simple: just ask two idiots to make random (but legal) moves until one of them wins. Then return the game to position X and have them play again. And again. And again and again. Sometimes they start with move A, and sometimes B. After your speedy but not-so-clever friends have played a few thousand games, examine the record. Is an idiot (with an idiot for an opponent) more likely to win with move A or move B? Those simulated games will give you the answer. Here’s the amazing thing: the idiot’s best move is your best move too. Don’t ask one clever mouse to solve the maze. Release ten thousand stupid mice and follow the lucky ones. This is what cheap computation buys you.

What’s beautiful about this approach is that it’s completely free of strategy. You don’t need to build up special knowledge structures about any particular game. You just need to know what moves are legal and how the game ends.

TIC TAC TOE
As soon as I read about this technique, I wanted to try it in MATLAB. So let’s make a program that can play Tic Tac Toe (also known as Naughts and Crosses). I’ve written Tic Tac Toe programs in MATLAB before. I’ve tried to make them clever and I’ve tried to make them learn. It’s not that hard. What’s fun about this Monte Carlo approach is that, with minimal effort I can teach it a new game. In fact, it makes playing lots of games easy. With a little object-oriented programming, you can write a generic game-playing harness. Then you just need to plug in some code that knows a few rules, and presto! You’ve got an instant game-playing program.

Here’s what I did. I made a class called TicTacToe that knows the rules of the game and how to draw the board. Then I wrote a function called botMoves that can look at the game object and make the next move. The separation is very clean. All of the Monte Carlo logic mentioned above lives in botMoves.

I only need a short script to have the bot play itself.

game = TicTacToe;
nSimulatedGames = 1000;
while ~game.isGameOver
botMoves(game,nSimulatedGames);
end

The variable nSimulatedGames refers to the number of simulated games we’ll ask our idiot friends to play for each potential move. Here’s an animation of what it looks like in action.

As it happens, the computer always ties itself. That’s actually good news, since Tic Tac Toe is unwinnable if your opponent is the least bit clever. So our bot is smart enough to prevent itself from winning. A little play-testing shows that it’s smart enough to avoid losing to a human too. But if we prefer, we can make the program less competitive by lowering the number of simulated games it plays. If I only let it run ten simulated games for each possible move, I can beat it easily.

I haven’t displayed much of my code here in the blog, but you can get your hands on it at this GitHub repository: Monte-Carlo-Games. Here is the TicTacToe class, and here is the botMoves function.

NEXT WEEK
This is the first of a two-part post. Next time we’ll show how quickly we can adapt our simple Tic Tac Toe harness for other games. We’ll also bring a community element into our programming. We’ll use Cody to source some of the tricky parts of our coding effort!

]]>UncategorizedVisualizing a Simple Saddle Point Algorithm in MATLAB
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/MhuDnGJuEaE/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2015/01/06/visualizing-a-simple-saddle-point-algorithm-in-matlab/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201501/896/62009828001_3993862324001_356-MinMax.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 03:50</span>
</div>
</a></div><p>There are lots of places in optimization and game theory where you would want to identify a saddle point on a two-dimensional surface. In this example we find and visualize the saddle point of a surface. We are trying to maximize the value of the surface by our choice in... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2015/01/06/visualizing-a-simple-saddle-point-algorithm-in-matlab/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1559Tue, 06 Jan 2015 14:34:34 +0000

]]>Prime Spiral
http://feedproxy.google.com/~r/mathworks/moler/~3/LJdQBY5pvtA/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/cleve/primespiral_10.png"/></div><p>The prime spiral was discovered by Stanislaw Ulam in 1963, and featured on the cover of Scientific American in March, 1964.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2015/01/05/prime-spiral/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1154Mon, 05 Jan 2015 17:00:51 +0000

The prime spiral was discovered by Stanislaw Ulam in 1963, and featured on the cover of Scientific American in March, 1964.

Stanislaw Ulam (1909-1984) was an eminent Polish-American mathematician who spent most of his career at Los Alamos, but who also had connections with several American universities, especially the University of Colorado. His many illustrious colleagues included von Neumann, Teller, Metropolis, Erdos, Richtmyer, Rota, and Kac.

Ulam had a wide range of mathematical interests. He was deeply involved in the early development of computers at Los Alamos. He is one of the primary developers of the Monte Carlo method. Later in his career he made contributions to mathematical biology and medicine. He wrote popular books, including "Adventures of a Mathematician" and "Analogies Between Analogies".

The Prime Spiral

As the story goes, Ulam was bored during a seminar in 1963 and began to doodle in his notebook. He numbered the integer lattice points in the plane in a spiral fashion and then highlighted the prime numbers. He was surprised to see that with this numbering the primes were not distributed randomly throughout the plane, but appeared to fall on diagonal segments.

10-by-10

Here is the 10-by-10 prime spiral, showing the spiral numbering and the 25 primes among the first 100 integers. The most prominent diagonal segment includes the primes 19, 7, 23, 47, and 79.

Piecewise Quadratics

Ulam's numbering at any point $(i,j)$ in the plane is a piecewise quadratic function of $i$ and $j$. For example, the values along the upper and lower halves of the antidiagonal, or their offsets by one, (depending upon whether the size is odd or even), are perfect squares. So there are certainly not any primes along these diagonals.

This reminds me of the quadratics that generate primes. The most famous is Euler's polynomial,

$$ p(n) = n^2-n+41 $$

This generates primes for $n$ from 1 through 40, but, obviously, not for $n$ equal to 41.

The Euler primes themselves don't show up in the Ulam prime spiral. But each of the diagonal segments corresponds to a quadratic polynomial, although probably one that generates fewer than 40 primes.

Improved Spiral Function

The MATLAB demos directory has a spiral function that generates Ulam's numbering. I'm afraid that the code is ancient, and pretty bad. It's not too hard to write better code. This generates successively larger spirals by rotating one of a given size 180 degrees and then appending another column at the right and another row at the bottom to produce the next larger size.

type improved_spiral

function S = improved_spiral(n)
% SPIRAL SPIRAL(n) is an n-by-n matrix with elements ranging
% from 1 to n^2 in a rectangular spiral pattern.
S = [];
for m = 1:n
S = rot90(S,2);
S(m,m) = 0;
p = m^2-m+1;
v = (m-1:-1:0);
S(:,m) = p-v';
S(m,:) = p+v;
end

Prime Spiral Is Spy

Once we have the spiral numbering, the prime spiral is our old friend, the spy plot.

type primespiral_core.m

function primespiral(n)
% PRIMESPIRAL PRIMESPIRAL(n) is Ulam's prime spiral.
P = zeros(n,n);
P(primes(n^2)) = 1;
S = improved_spiral(n);
P = P(S);
spy(P)

Increasing Sizes

Here are prime spirals of three different sizes.

Animated Prime Spiral

The function primespiral in my textbook Numerical Computing with MATLAB includes a provision for starting the numbering at a value larger than one. Animating this feature emphasizes the diagonal nature of the prime locations. The title on this plot is the starting value of the numbering.

]]>Season’s Greetings, 2014
http://feedproxy.google.com/~r/mathworks/moler/~3/6ODUiRvEqdU/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/cleve/greeting_5b.png"/></div><p>http://www.mathworks.com/moler/ncm/greetings.m... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2014/12/25/seasons-greetings-2014/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1150Thu, 25 Dec 2014 17:00:49 +0000

Season's Greetings, 2014

help greetings

% greetings Seasonal holiday fractal.% greetings(phi) generates a seasonal holiday fractal that depends% upon the parameter phi. The default value of phi is the golden ratio.% What happens for other values of phi? Try both simple fractions% and floating point approximations to irrational values.

]]>Editing an Existing Figure File in MATLAB
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/cUzKTKUFGnI/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2014/12/19/editing-an-existing-figure-file-in-matlab/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201501/2536/62009828001_3993752599001_355-Editing-Figure-File.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 02:08</span>
</div>
</a></div><p>The best way to modify a MATLAB figure is to just modify the code that generated it. I prefer this because it is better to be able to regenerate a figure from code if you want to keep modifying it later, and your workflow is visible for later inspection. Sometimes... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2014/12/19/editing-an-existing-figure-file-in-matlab/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1554Fri, 19 Dec 2014 20:02:24 +0000like this User on MATLAB Answers. I show a few techniques for modifying an existing figure when you do not have the original code that made it.

]]>Displaying Simulation Data
http://feedproxy.google.com/~r/SethOnSimulink/~3/Vu12vSFHYT0/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/seth/2014Q4/withPortValues.png"/></div><p>This week, I am working on a project with MathWorks consulting where I have been asked to implement in Simulink a simulator already existing in a different language.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2014/12/18/displaying-simulation-data/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4346Thu, 18 Dec 2014 13:24:04 +0000This week, I am working on a project with MathWorks consulting where I have been asked to implement in Simulink a simulator already existing in a different language.

As you can imagine, to debug and validate my work, I need to compare a lot of signals in my model with the original data I received.

As a first step, I am trying to validate the static behavior of the model. This means that I simulate the model for only one time step and validate that the outputs are as expected for different inputs.

If the output does not match, then I need to inspect intermediary signals to figure out where the difference originates.

Display blocks

The first, most obvious, way to see the value of a signal in Simulink is the Display block. In my case, using Display blocks quickly made my model look like:

As you can see, the model quickly becomes a mess. And since those displays are only temporary, I will need to remove them later.

Value Label Display

Another option that I often find more convenient to quickly view the value of a signal for debugging is the Value Label Display.

In the Display menu, I can decide if I want the label to be displayed when hovering over the block, or when clicking on it.

There are also options to control how often the labels are updated, and their accuracy. In my case, I chose the long format to see all the digits available.

Once this is enabled, you can simply click on blocks to enable the displays and simulate the model to see the values.

Once your debugging is completed, you can remove all the labels in one click from the Display menu.

Now it's your turn

What is your preferred method to to view signal values during debugging? Let us know by leaving a comment here.

]]>Jahnke and Emde, Revisited
http://feedproxy.google.com/~r/mathworks/moler/~3/tGANWXM1DEY/
<div class="overview-image"><img src="http://blogs.mathworks.com/cleve/files/feature_image/hankelbw.jpeg" class="img-responsive attachment-post-thumbnail wp-post-image" alt="hankelbw"/></div><p>An incredible book, published in several editions from 1909 to 1933, by German mathematicians Eugene Jahnke and Fritz Emde, contains definitions and formulas for mathematical functions, hand-calculated tables of function values, and meticulous hand-drawn 2- and 3-dimensional graphs. An English edition was published by Dover in 1945.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2014/12/15/jahnke-and-emde-revisited/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1140Mon, 15 Dec 2014 17:00:01 +0000

An incredible book, published in several editions from 1909 to 1933, by German mathematicians Eugene Jahnke and Fritz Emde, contains definitions and formulas for mathematical functions, hand-calculated tables of function values, and meticulous hand-drawn 2- and 3-dimensional graphs. An English edition was published by Dover in 1945.

I will include only two of their illustrations in this blog. Here is one of their graphs of the gamma function in the complex plane. The surface is the magnitude of the function, with its poles at the integers on the negative real axis. The contour lines are the modulus and phase. I can't imagine how they made this plot in 1933.

MATLAB Gamma Function

Here is my attempt to have MATLAB recreate the Jahnke and Emde plot. It is only partially successful. I do have little ellipses to cut off the spikes -- that was tricky. But I haven't begun to label the contours, or shade the edge faces.

Jahnke and Emde Bessel Function

I love this plot. It is the complex-valued Bessel function of order zero, also known as the Hankel function. There is a pole, a zero, and a saddle point. Again, imagine doing this by hand. Or, having some grad student do it.

MATLAB Bessel Function

Here we can use color, and our new parula color map, in place of 3d. But we can't get the nice feel of the hand drawn lettering.

[X,Y] = meshgrid(-4:1/128:2,-1.5:1/128:1.5);
H = besselh(0,X+i*Y);
s = 3.2;
contourf(X,Y,abs(H),0:0.2:s)
hold on
c2 = s*(angle(H)/(2*pi)+0.5);
r2 = s*((-180:10:180)/360+0.5);
contour(X,Y,c2,r2,'k-')
hold off
axis equal
axis tight
grid on

Reference

Eugene Jahnke and Fritz Emde, Tables of Functions with Formulas and Curves, (4th ed.), Dover, 1945

]]>A New Colormap for MATLAB – Part 4 – The Name
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/Aj6I9hLy3cM/
<div class="overview-image"><img src="http://blogs.mathworks.com/steve/files/tropical-parula-400x300.jpg" class="img-responsive attachment-post-thumbnail wp-post-image" alt="tropical-parula-400x300"/></div><p>In the summer of 2013 we were closing in a choice for the new MATLAB colormap. We were down to tweaking and fine-tuning.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2014/12/10/a-new-colormap-for-matlab-part-4-the-name/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1259Wed, 10 Dec 2014 23:06:02 +0000

In the summer of 2013 we were closing in a choice for the new MATLAB colormap. We were down to tweaking and fine-tuning.

But ... we needed a name!

For my many experiments, I had fallen back on an old graduate school habit of naming things after characters from Lord of the Rings. So I had filenames such as gandalf_20130623a.m and faramir_20130712b.m. (Faramir is my favorite LOTR minor character.) I certainly knew that wasn't going to work for the final name. How to choose one?

I looked over the names of the existing colormaps:

jet

hsv

hot

cool

spring

summer

autumn

winter

gray

bone

copper

pink

lines

colorcube

prism

flag

Well, we've got a few identifiable themes in there:

seasons (summer, autumn, winter, spring)

temperatures (cool, hot)

materials (bone, copper)

hues (pink, gray)

colorspaces (hsv, colorcube - sort of)

And there are some oddball names, such as jet, lines, prism, and flag. None of that seemed inspirational for naming a new colormap.

I decided to look for something descriptive. But descriptive of what?

To remind you, here's what the new colormap looks like:

showColormap(parula,'bar')

I picked the main colors (to my eye, these are blue, green, orange, and yellow) in the new colormap and started doing searches using these color names and different kinds of objects. Animals seemed obvious. I actually started with fish, but that got nowhere fast.

Then I tried birds, and up popped the tropical parula:

As Wikipedia describes it, the tropical parula "has mainly blue-grey upperparts, with a greenish back patch and two white wingbars. The underparts are yellow, becoming orange on the breast."

Perfect! I didn't really think the team would go for it, though. I sent around an email with the name parula and a picture of the bird to a small group of people working on the visual appearance changes for the new MATLAB graphics system. Somewhat to my surprise, everyone said they liked it. Later, a larger group of senior MATLAB designers reviewed it, and they also liked it. So the name stuck.

That left us with one big problem, however. How is parula pronounced?

I don't know. In trying to find a definitive answer, I have only managed to confuse myself. As a result, I don't even always pronounce it the same way myself. In looking at various references for pronunciation of bird names, I have seen all of these variations:

pah-ROO-lə

PAIR-yə-lə

PAIR-ə-lə

PAR-yə-lə

I know two amateur birders at MathWorkers who told me definitively how to pronounce it. Of course, they each gave me a different answer.

I guess that most American English speakers would choose the first variation, which has the accent on the second syllable. That has certainly been the case at MathWorks headquarters in Massachusetts.

So I'll tell you what I tell MathWorkers here: you can pronounce it however you like!

]]>ColormapSelectively Changing an Image in MATLAB with Logical Indexing
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/qwOCtoHyq-M/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2014/12/05/selectively-changing-an-image-in-matlab-with-logical-indexing/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201501/3456/62009828001_3993807956001_354-Applied-Binary-Index.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 02:51</span>
</div>
</a></div><p>A MATLAB user on MATLAB Answers wanted to know how to change all the pixels in an image that were below a certain threshold. The answer was correct, but very terse. I wanted to expand upon the question with an explanation of logical indexing in MATLAB.
... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2014/12/05/selectively-changing-an-image-in-matlab-with-logical-indexing/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1547Fri, 05 Dec 2014 14:34:09 +0000MATLAB user on MATLAB Answers wanted to know how to change all the pixels in an image that were below a certain threshold. The answer was correct, but very terse. I wanted to expand upon the question with an explanation of logical indexing in MATLAB.

]]>MathWorks Logo, Part Five, Evolution of the Logo
http://feedproxy.google.com/~r/mathworks/moler/~3/m5ZthCvSbg8/
<div class="overview-image"><img src="http://blogs.mathworks.com/cleve/files/feature_image/montage.jpeg" class="img-responsive attachment-post-thumbnail wp-post-image" alt="montage"/></div><p>Our plots of the first eigenfunction of the L-shaped membrane have changed several times over the last fifty years.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2014/12/01/mathworks-logo-part-five-evolution-of-the-logo/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1129Mon, 01 Dec 2014 17:00:12 +0000

Our plots of the first eigenfunction of the L-shaped membrane have changed several times over the last fifty years.

Ph.D. thesis. Two dimensional contour plot. Calcomp plotter with ball-point pen drawing on paper mounted on a rotating drum.

1968, University of Michigan

Calcomp plotter. Math department individual study student project. Original 3d perspective and hidden line algorithms. Unfortunately, I don't remember the student's name.

1985, MATLAB 1.0

The first MathWorks documentation and the first surf plot. Apple Macintosh daisy wheel electric typewriter with perhaps 72 dots-per-inch resolution.

1987, MATLAB 3.5

Laser printer with much better resolution.

1992, MATLAB 4.0

Sun workstation with CRT display and color. Hot colormap, proportional to height.

1994, MATLAB 4.2

A hand made lighting model. I have to admit this looks pretty phony today.

1996, MATLAB 5.0

Good lighting and shading. But at the time this was just a nice graphics demonstration. It was a few more years before it became the official company logo.

1990, Nantucket Beach

A prototype 3d printer. Pure granular silicon, a few dozen grains per inch. (Thanks to Robin Nelson and David Eells for curating photos.)

]]>MATLAB Onramp
http://feedproxy.google.com/~r/mathworks/desktop/~3/vaEJrpjnudU/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/community/files/blogpost.gif"/></div><p>In the last post, Wendy Fullam told us about MATLAB Examples, where you could see working examples and pick apart the code to see how they worked. But suppose you wanted more? Suppose you wanted an environment that could teach you, step by step, how to get started coding in... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/community/2014/12/01/3085/">read more >></a></p>http://blogs.mathworks.com/community/?p=3085Mon, 01 Dec 2014 13:00:05 +0000In the last post, Wendy Fullam told us about MATLAB Examples, where you could see working examples and pick apart the code to see how they worked. But suppose you wanted more? Suppose you wanted an environment that could teach you, step by step, how to get started coding in MATLAB? If that’s the case, then Chaitanya Chitale has just what you need. Chaitanya is a MATLAB Training Content Developer at MathWorks.

Introducing MATLAB Onramp

by Chaitanya Chitale

Are you new to MATLAB and looking for a way to get started? If so, check out the new MATLAB Onramp course.

MATLAB Onramp is an online course that provides a brief introduction to the MATLAB language. The course gives you hands-on MATLAB experience via the use of an integrated, web-based version of MATLAB, as shown below.

As you work through the course you can try out different solutions to the problems asked and get immediate feedback.

Course content and duration

The MATLAB Onramp course covers the basics of importing data, manipulating arrays, creating visualizations, and much more.

The course takes approximately 2 hours to complete. However, you can take the course at your own pace. After starting it, you can leave and come back to it at any time.

Ready to get started?

The MATLAB Onramp course is complimentary with your purchase of MATLAB. To get started, head over to the MATLAB Academy site: https://matlabacademy.mathworks.com.

You can also access the course from MATLAB R2014b by selecting Help → MATLAB Academy from the Home tab.

]]>UncategorizedTime zone terminology
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/xvvewGL408E/
<p>An email thread about a date-related design question had me thinking about terminology this week. (Note: There's nothing about image processing in today's post. Sorry!)... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2014/11/21/time-zone-terminology/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1251Fri, 21 Nov 2014 18:32:26 +0000

An email thread about a date-related design question had me thinking about terminology this week. (Note: There's nothing about image processing in today's post. Sorry!)

Specifically, I was thinking about the term time zone. In the U.S. you might see a date written this way:

November 21, 2004 1:02 PM (EST)

EST means Eastern Standard Time in the U.S. Or you might see a date written like this:

11/21/2004 8:30 (-04:00)

Many people refer to EST or -04:00 as the "time zone." However, because I work across the hall from the chief designer and implementer (Hi, P^2) of the new datetime suite of functions in MATLAB R2014b, I know that these are more properly referred to as time zone offsets.

A time zone is a much more complex concept than a simple offset from UTC. For example, a time zone encompasses both current and historical rules for using Daylight Saving Time.

The time zone offset for a given time depends not only on the time zone but also both the year and the day-of-year.

The dependence on the day-of-year is pretty obvious. In the fall season where I live, the time zone offset depends on whether the day is before or after the first Sunday in November.

The dependence of time zone offset on the year is somewhat less obvious. It turns out that Daylight Saving Time rules change occasionally. In the U.S. the rules were last changed in 2007. Let's compare the time zone offset for noon on Halloween for 2007 and 2006.

Now for just a little fun with MATLAB R2014b date calculations. Most of the U.S. changes from Daylight Saving time to Standard Time at 2:00 AM on the first Sunday of November. Given the year, how can we compute that?

Here's just one way. Start with 2:00 AM on November 1 of the specified year.

year = 2014;
d = datetime(year,11,1,2,0,0)

d =
01-Nov-2014 02:00:00

Now use dateshift to move to the first Sunday.

dateshift(d,'dayofweek','Sunday') % You can also use 1 instead of 'Sunday'.

ans =
02-Nov-2014 02:00:00

My horoscope for today said "Your next blog will include an example of an anonymous function for no good reason." So here it is. Let's make an anonymous function that returns the exact time for the transition from Daylight Saving Time to Standard time in the U.S.

f = @(year) dateshift(datetime(year,11,1,2,0,0),'dayofweek','Sunday')

f =
@(year)dateshift(datetime(year,11,1,2,0,0),'dayofweek','Sunday')

f(2014)

ans =
02-Nov-2014 02:00:00

f(2015)

ans =
01-Nov-2015 02:00:00

f(2016)

ans =
06-Nov-2016 02:00:00

Of course, this function only works for years between 2007 and whenever the U.S. Congress changes its mind again!

Next week, I promise, I'll get back to writing about image processing. Or maybe colormaps. Unless I write about something else.

]]>UncategorizedMATLAB: In an assignment A(I) = B, the number of elements in B and I must be the same.
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/yLPG6f3PA3s/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2014/11/20/matlab-in-an-assignment-ai-b-the-number-of-elements-in-b-and-i-must-be-the-same/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201411/3371/62009828001_3911528083001_353-scalar-vector-matrix-thumb.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 02:33</span>
</div>
</a></div><p>The very common MATLAB error: “In an assignment A(I) = B, the number of elements in B and I must be the same.” occurs when one side of an assignment has a dfferent dimension than the other. So for instance the left side of the assignment might be... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2014/11/20/matlab-in-an-assignment-ai-b-the-number-of-elements-in-b-and-i-must-be-the-same/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1538Thu, 20 Nov 2014 18:52:06 +0000

]]>A New Colormap for MATLAB – Part 3 – Some Reactions
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/tkfFGbGmT1w/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/steve/2014/colormaps_part_3_05.png"/></div><p>I hesitated for some time today about where to go next in my discussion about the new default colormap in MATLAB. I finally decided to spend a little time on the reactions so far. Those reactions have ranged from "Great, what took you so long?" to "I hate it; give... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2014/11/12/a-new-colormap-for-matlab-part-3-some-reactions/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1243Wed, 12 Nov 2014 16:39:50 +0000

I hesitated for some time today about where to go next in my discussion about the new default colormap in MATLAB. I finally decided to spend a little time on the reactions so far. Those reactions have ranged from "Great, what took you so long?" to "I hate it; give me back jet!" In the Twitterverse, one scientist proclaimed "End of an era!" Another simply said "yuck." (Jet is still in MATLAB, by the way.)

Today I want to dig into some of those reactions a little bit. Along the way I'll show some experiments with a test pattern to explain some of my comments about jet.

To remind you, last time I summarized 15 years or so of criticisms about rainbow colormaps:

Rainbow colormaps confuse viewers because there is no natural perceptual ordering of the spectral colors.

Rainbow colormaps obscure small details in the data.

Rainbow colormaps mislead viewers by suggesting data features that are not really there. These "phantom features" often take the form of false boundaries. This effect can falsely segment the data.

Rainbow colormaps lose critical information about high and low data values when displayed or printed on a gray-scale device.

Rainbow colormaps can be difficult to interpret for some color-impaired viewers.

I have seen a claim that that rainbow incorporates a universally understood order based on associating blue with cold and red with hot. There's no doubt that this order makes sense for some people (although blue flame is hotter than red flame?). I quibble with the "universally understood" part, though. The dominant evidence is that individuals do not universally order the spectral hues the same way unless they refer to spectral order mnemonic, such as "Roy G. Biv." It is well-established, however, that individuals can easily and consistent order colors by relative brightness.

Another commenter said that the rainbow colormap has higher contrast than parula. Yes, this is true. Rainbow goes all the way from dark to bright in the bottom third of the data range, and it goes all the way from bright to dark in the top third of the data range. For some parts of the data range, therefore, there is higher visual contrast than with parula, which uses the entire data range to go from dark to bright. A price paid for this increased contrast, however, is the visual illusions and confusion caused by the reversed roles of dark and bright in different parts of the colormap.

Matteo Niccoli gives the following nice example of a rainbow-based visual illusion. The arrow in the right-hand image below points to something that gives the strong but false impression of being a depression in the bone structure.

And I showed another visual confusion example last time with these visualizations of two different sinusoidal patterns:

These visualizations give the impression that the two sinusoidal patterns have opposite phase, but in reality they have exactly the same phase.

Another commenter said that, with sufficiently familiarity, the bright "stripes" that are characteristic of the rainbow colormap can be an advantage for certain kinds of data. That might be true in some particular cases, but I don't think it's true in general. Let me show an example, using a test pattern, that illustrates how much difficulty the cyan and yellow stripes can cause in interpreting what should be easily recognizable data.

url = 'http://blogs.mathworks.com/images/steve/2014/eia-2000-cropped.png';
I = imread(url);
imshow(I)
colormap(jet)

Let me label some portions of the visualization and then ask some quetions about it.

At periodic intervals, the dark blue vertical bars have a bright outline (see the points marked "A"). What is the data significance of these bright outlines? Especially when compared with the dark blue vertical bars that do not have a bright outline?

The vertical bars in between the "A" marks don't have bright outlines. Those bars alternate between dark blue and dark red, a low-contrast combination that is a little difficult to focus clearly on (a red-blue visual illusion called Chromostereopsis).

Above the "B", I see very thin vertical dark bars. When the bars are this thin, it's difficult to distinguish between the dark red and dark blue, which is unfortunate because those colors represent opposite ends of the data range. The bright outlines between the bars is broken up by a diagonal pattern of blue stripes.

The numbers around the border are oddly broken up and hard to decipher.

Let's investigate why some vertical bars are outlined and some not by looking at the data cross-section along the white path marked "C".

c = improfile(I,[77 160],[225 225]);
figure
plot(c,'*-')

I'll superimpose two lines that show the data values corresponding to the bright cyan and bright yellow in the jet colormap.

hold on
plot([1 84],[1 1]*0.375*255,'c','LineWidth',5)
plot([1 84],[1 1]*0.625*255,'y','LineWidth',5)
hold off

It is now apparent why some vertical bars are outlined but not others. In each transition from low to high and from high to low, there is one data value that is in between. When that one transitional data value happens to be close the data values mapping to cyan or yellow, the vertical bar is shown with a bright outline. Otherwise, there is no outline. Although the effect is visually dramatic and makes some regions of the data look quite different from others, the visual difference has no real data significance.

OK, now let's show this data using the parula colormap.

imshow(I)
colormap(parula)

This visualization is trivially easy to interpret. The vertical bars in the middle are all the same. The digits are easy to read. The visual difference between low data values and high data values has high contrast and is easy on the eyes. And the data above the point marked "B" is seen to have a completely different form than the way it appeared with jet. The stripes have a lower spatial frequency and are seen to be slightly fanning out in direction instead of appearing vertical.

A valid complaint about this example is that it is completely artificial, with data rather unlike what people are likely to be looking at with jet or parula. I'd like to suggesting thinking about it this way, though: If jet can take familiar, easily interpretable data and render it almost unrecognizable, what can jet do to your own data?

Finally, I'll mention that a blog reader commented that using jet made it easy to distinguish positive and negative values from each other in a data set whose values are centered around zero. And yes, jet can be better than parula for that. Jet is bright and relatively unsaturated in the middle, and it is dark and saturated at both ends. These are characteristics often shared by colormaps designed specifically for highlighting a central, specific value and distinguishing the values above and below that. These are called divergent colormaps. In a future post I'll show where to find some good divergent colormaps and how to use them in MATLAB.

]]>ColormapRotate Axes Labels in MATLAB
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/zNHK5V_xQy0/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2014/11/05/rotate-axes-labels-in-matlab/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201411/2371/62009828001_3911528086001_352-rotate-Tick-Labels-thumb.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 01:35</span>
</div>
</a></div><p>The new release of MATLAB R2014b brings a new graphics engine. This allows for some new capabilities such as rotating tick labels and using dot notation in referencing graphics properties. These are a couple of small but very useful parts of this release.
... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2014/11/05/rotate-axes-labels-in-matlab/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1523Wed, 05 Nov 2014 19:53:04 +0000

]]>Introducing MATLAB Examples
http://feedproxy.google.com/~r/mathworks/desktop/~3/Fp7Foc5EI8A/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/community/files/related.png"/></div><p>This week, guest-blogger Wendy Fullam is trumpeting the arrival of a new Support page feature. Wendy is the Technical Marketing Product Manager for MathWorks online support and website search program.
MATLAB Examples
by Wendy Fullam
Have you ever found yourself thinking, “I wish there was a place where I could go to see... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/community/2014/11/03/introducing-matlab-examples/">read more >></a></p>http://blogs.mathworks.com/community/?p=3046Mon, 03 Nov 2014 15:00:17 +0000This week, guest-blogger Wendy Fullam is trumpeting the arrival of a new Support page feature. Wendy is the Technical Marketing Product Manager for MathWorks online support and website search program.

MATLAB Examples

by Wendy Fullam

Have you ever found yourself thinking, “I wish there was a place where I could go to see a bunch of MATLAB examples…”?

MathWorks has just what you’re looking for. We just launched a big new feature on the Support page: MATLAB Examples. On it, you can discover thousands of code examples for MATLAB and Simulink: data manipulation, signal processing, machine learning, statistics, finance, computational biology, finance, you name it. Maybe you’ll even discover a few things that you didn’t know were possible.

MATLAB Examples is a single destination to find high-quality code examples, including those authored by both MathWorks staff and contributors to the File Exchange. Here’s what the main page looks like.

What can you do here?

You can browse by topic area:

You can search by keyword, and filter the results:

After finding an example that you want to use, you can quickly see what product contains that example:

or access the code:

You can also find related examples:

The Help Team at MathWorks hopes you enjoy the new interface and welcomes your thoughts and reactions. Check out the new area and then let us know what you think in the comments below!

]]>UncategorizedGetting my ICIP 2014 workshop files
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/PTBpa_5RAXc/
<p>For ... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2014/10/29/getting-my-icip-2014-workshop-files/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1237Wed, 29 Oct 2014 12:51:50 +0000here are the files.]]>UncategorizedUsing DBstack in MATLAB to get a stack trace
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/kHpiyEfW1h4/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2014/10/24/using-dbstack-in-matlab-to-get-a-stack-trace/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201411/1861/62009828001_3911528089001_351-dbStack-to-Debug-thumb.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 03:22</span>
</div>
</a></div><p>When you have timers and callbacks in your code, it can be difficult to find out when certain errors are occuring. Use DBstack to find the debug stack trace in MATLAB. This will allow you to find out what function called your current function, and what function called... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2014/10/24/using-dbstack-in-matlab-to-get-a-stack-trace/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1515Fri, 24 Oct 2014 17:19:01 +0000

]]>ICIP 2014
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/nyFjEuJQa6M/
<p>
Next week I'll be at ICIP 2014. That's the IEEE International Conference on Image Processing in Paris.
... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2014/10/24/icip-2014/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1232Fri, 24 Oct 2014 13:17:53 +0000
Next week I'll be at ICIP 2014. That's the IEEE International Conference on Image Processing in Paris.

If you're going to be there, come see my MATLAB Techniques for Image Processing talk on Wednesday, 12:45 - 13:45. Or stop by the MathWorks booth.

Now you'll have to excuse me while I finish preparing my talk.

]]>UncategorizedMATLAB and Git
http://feedproxy.google.com/~r/mathworks/desktop/~3/TfvkVQ-Xepc/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/community/files/branching-300x223.png"/></div><p>This week we hear from Toshi Takeuchi about how to take advantage of MATLAB’s recent improvements to Git integration. Toshi is a Senior Marketing Manager at MathWorks.
Quick Introduction to Git with MATLAB
by Toshi Takeuchi
One of the new R2014b features that deserves your attention is Git integration. Git is a source... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/community/2014/10/20/matlab-and-git/">read more >></a></p>http://blogs.mathworks.com/community/?p=3018Mon, 20 Oct 2014 13:00:31 +0000This week we hear from Toshi Takeuchi about how to take advantage of MATLAB’s recent improvements to Git integration. Toshi is a Senior Marketing Manager at MathWorks.

Quick Introduction to Git with MATLAB

by Toshi Takeuchi

One of the new R2014b features that deserves your attention is Git integration. Git is a source control system (also known as version control or source code management system) that enables collaborative software development. Why does that matter to you? Programming is an essential skill in many technical fields even outside computer science, and some universities now offer software carpentry workshops to enhance coding skills for researchers. Source control is one of those essential skills in software carpentry.

Until now, you may have tinkered alone with code you needed for your project. However, there are other people who may be working on similar problems and they may be writing similar programs. Source control enables you to work with other people so that you don’t have to do it all alone. Collaboration lets you be more productive in other aspects of your project.

Even if you don’t care about such collaboration, wouldn’t it be cool to share your personal project and see other people using it? They may even fix bugs and improve your code for you!

GitHub is one of the most popular websites that host Git repositories. The best place to share your MATLAB projects is File Exchange because of its popularity with the MATLAB user community. And guess what – File Exchange is integrated with GitHub! Now you see the connection?

Basic terminology

What is a Git repository? A repo (repository) is a directory that holds your source code and any associated files. Local repos are on your local drive, and the remote repos are on GitHub or other hosts, and you sync the local repos to remote repos as you write your code. You can start with either the local or remote repos, but in this example I am going to start with a remote repo.

The process looks like this for a single developer:

Create or fork a repo on GitHub

Clone the repo to your local drive – this is your local repo

Add your files to the local repo

Sync your local repo to remote repo

Repeat this process to keep your source code in sync as you write more code

Share your GitHub repo on File Exchange

What is forking?

When you talk about Git, you cannot go without mentioning “forking”. In the simplest terms forking means copying someone else’s public repo on the remote server, rather than starting a repo from scratch. In practice forking is used as a way to contribute to the existing projects or to start a new project using an existing project as a starting point. Once you make changes to your forked project, you can send a merge request to the original developer, and your changes may be accepted and incorporated into the main project.

Forking enables a flexible distributed style of collaboration and number of forks you have on your project acts as a measure of popularity – similar to the count of likes or followers on Facebook or Twitter. The social aspect of forking is an interesting topic on its own, but we need to skip it for this post.

Creating a repo on GitHub is very easy – just follow these instructions. From this point on I will assume you named your repo Hello-World and initialized it with a README file. Please note that you can only create a public repo with a free account.

Cloning the repo to your local drive with MATLAB

Until recently, you needed to use the command line tool for this step, but starting with R2014b we can just use MATLAB’s Current Folder window. No more Git commands like git init, git status, git add, or git commit!

Open your copy of MATLAB and create an empty folder. Right-clicking the empty space in the Current Folder window to bring up a contextual menu, and select Source Control > Manage Files.

This will open a new dialog box: Manage Files Using Source Control.

For Select control integration, choose Git

For Repository path, click Change

You now see a new dialog box: Select a Repository. Copy and paste the URL of the remote repo you just created. You can find the URL in the right sidebar of your new repo on GitHub.

You choose either SSH or HTTPS depending on how you setup your authentication on GitHub.

Click Validate. You may be asked for your login password for authentication. You can close the dialog box when your path is validated.

Back in Manage Files dialog box, the sandbox should be already set to your current folder. All you need to do now is hit Retrieve.

You have now successfully cloned the remote repo to your local drive. Check your Current Folder window. You should see just one file – README.md, but with a green circle next to it. This is just a text file but you can apply wiki-like syntax called markdown to make it appear like a regular web page when viewed on GitHub. README serves as the front page of your repo on GitHub.

Adding files to your local repo

Let’s add a new MATLAB script file helloworld.m. It will appear with a blank circle – it means it is not added to Git source control yet. To add it to Git, right-click on the file and select Add to Git. The empty circle changes to “+” symbol. When you modify a file already under source control, the symbol becomes a blue square.

Taking a snapshot with commit

You can continue editing files as you like, but at some point, you want to take a snapshot of the edits you made. That’s when you do a commit. You can select any empty space in the Current Folder window to bring up the contextual menu and select Commit to Git Repository. This will bring up a dialog box where you can add your comment about the changes you made since the last commit. Comments will be helpful to keep to track of your changes and revert back to earlier commits if necessary.

Synching your local repo to remote repo

When you commit, the snapshot is saved in the local repo, but it is also a good idea to mirror the changes to the remote repo as well. To do so, bring up the contextual menu by right-clicking an empty space in the Current Folder window and select Push. That will push your changes to the remote repo. You may need to enter your password.

Branching and Merging

The real power of source control comes from the ability to create multiple branches from your project. By default, you have a branch called “master” in your repo. You can create a new branch from the master branch, makes changes, and then merge those changes back to the master. This mechanism is used for working on new experimental features without affecting the working code on the master. You can branch and merge in MATLAB but the details are beyond the scope of this post.

Closing

If you have been curious about Git but put off by its complicated command line interface, Git integration in R2014b makes Git much more pleasant and approachable. I hope this quick introduction motivates you to take advantage of this new feature. When you do, please don’t forget to post your project to the File Exchange. To learn more about Git, it is actually helpful to start with reviewing the underlying concept about how Git works.

]]>UncategorizedMATLAB r2014bGraphics has launched
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/aGtt7DWFa-Y/
<p>
We are very excited to announce that MATLAB r2014b introduces a new graphics system. Graphics developer, Mike Garrity, is launching a graphics specific blog. Here is the landing page for this release of MATLAB and for the new graphics system in general.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2014/10/09/matlab-r2014bgraphics-has-launched/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1503Thu, 09 Oct 2014 14:02:06 +0000

We are very excited to announce that MATLAB r2014b introduces a new graphics system. Graphics developer, Mike Garrity, is launching a graphics specific blog. Here is the landing page for this release of MATLAB and for the new graphics system in general.]]>
Topic: Cool featureAcquire Data from Android Device Sensors with MATLAB Mobile
http://feedproxy.google.com/~r/mathworks/desktop/~3/Soozdn0ekkA/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/community/files/matlab_mobile_sensors_04.png"/></div><p>
With the new MATLAB® Support Package for Android™ Sensors, you can now use MATLAB Mobile™ to acquire data from the sensors on your Android device. This data can be sent to a MATLAB session running on your computer for further analysis and visualization.
Contents
What data, you ask?
Viewing Sensor Data
Analyze Data with... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/community/2014/10/06/acquire-data-from-device-sensors-with-matlab-mobile/">read more >></a></p>http://blogs.mathworks.com/community/?p=2905Mon, 06 Oct 2014 15:00:16 +0000
With the new MATLAB® Support Package for Android™ Sensors, you can now use MATLAB Mobile™ to acquire data from the sensors on your Android device. This data can be sent to a MATLAB session running on your computer for further analysis and visualization.

On Android devices, MATLAB Mobile supports data acquisition from motion sensors like the accelerometer as well as positional sensors like the GPS. A list of all sensors is shown below.

Viewing Sensor Data

You can access these sensors by selecting the Sensors option from the drop-down menu in MATLAB Mobile. You can tap on a sensor to enable it and view related measurements. The screenshot below is the result of turning on the Accelerometer and Magnetometer.

Analyze Data with MATLAB

Displaying this data is cool, but to make this truly useful, you will want to perform further analysis and processing. Fortunately, the MATLAB Support Package for Android Sensors helps you do just that! It enables you to send sensor data to a MATLAB session on your computer. To do this:

Connect MATLAB Mobile to your computer with the MATLAB Connector. For more information on how to do this, refer to the getting started instructions here. This feature is only supported on MATLAB R2014a and later, so make sure you are on a compatible version.

Install the MATLAB Support Package for Android Sensors. Choose Add-ons from the MATLAB Toolstrip, and then choose Get Hardware Support Packages. This will open the support package installer. Choose Android Sensors from the list and follow the instructions.

To establish communication between the sensors on your device and MATLAB, create a mobiledev object, as follows:

m = mobiledev;

Example: Counting Steps by Capturing Acceleration Data

The mobiledev object facilitates communication between the sensors on your Android device and the MATLAB session running on your computer. Let’s explore this workflow through an example that illustrates the collection of acceleration data and using it to count the number of steps taken.

Step 1: Turn on the Accelerometer Once you have completed the 3 steps from the above section, go to MATLAB Mobile and turn on the accelerometer. You should see something akin to this:

You can also enable the sensor directly from MATLAB, by executing the following command:

m.AccelerationSensorEnabled = 1;

Step 2: Send Data to MATLAB Did you notice the enabled Start Sending button towards the bottom of your screen? Tap on it, and voila! You are now sending data to MATLAB.
Alternatively, you can start sending data directly from MATLAB, through the following command:

m.Logging = 1;

You can verify this in MATLAB, note the Current Sensor Values in the result:

Step 3: Stop Acquiring Data and Retrieve Logs Walk around your campus/home/floor with your device. Once you are satisfied, stop sending this data to MATLAB. You can either tap on the Stop Sending button on MATLAB Mobile, or issue the following command in MATLAB:

m.Logging = 0;

To retrieve the data, use the accellog variable:

[a, t] = accellog(m);

Step 4: Plot Raw Sensor Data Once you have retrieved the logged acceleration data, you can plot it in MATLAB:

plot(t, a);
legend('X', 'Y', 'Z');
xlabel('Relative time (s)');
ylabel('Acceleration (m/s^2)');

Calculate the magnitude to convert your X, Y and Z vectors to scalar values. Then, plot it.

x = a(:,1);
y = a(:,2);
z = a(:,3);
% Calculate and plot magnitude.
mag = sqrt(sum(x.^2 + y.^2 + z.^2, 2));
plot(t, mag);
xlabel('Time (s)');
ylabel('Acceleration (m/s^2)');

To remove constant effects such as gravity, you can subtract the mean from this data.

The plotted data is now centered on zero, and shows peaks which correspond to a step taken while walking.

Step 5: Count Number of Steps Taken To determine the number of steps taken, you can use to FINDPEAKS function from Signal Processing Toolbox. In this example, we are treating only peaks with a minimum height above one standard deviation as a step. This threshold should be tuned experimentally to match a person’s level of movement while walking, hardness of floor surfaces etc.

% Use FINDPEAKS to determine the local maxima.
minPeakHeight = std(magNoG);
[pks, locs] = findpeaks(magNoG, 'MINPEAKHEIGHT', minPeakHeight);

The number of steps taken is the number of peaks:

numSteps = numel(pks)

numSteps =
15

Finally, you can also identify these locations on your plot of acceleration magnitude data:

hold on;
% Place a red marker on the locations that correspond to peaks.
plot(t(locs), pks, 'r', 'Marker', 'v', 'LineStyle', 'none');
title('Counting Steps');
xlabel('Time (s)');
ylabel('Acceleration Magnitude, Gravity Removed (m/s^2)');
hold off;

Step 6: Clean Up Once you are done, make sure you turn off the acceleration sensor and clear the mobiledev object.

m.AccelerationSensorEnabled = 0;
clear m;

Try it out!

To learn more about acquiring data from sensors on your mobile device, refer to the following links:

]]>MATLAB MobileWelcome Mike!
http://feedproxy.google.com/~r/mathworks/desktop/~3/y1m20zJjqYE/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/graphics/files/feature_image/welcome_teapot_03.png"/></div><p>I am very excited to welcome a new blogger to MATLAB Central: Mike Garrity is going to be writing about graphics on the appropriately named Mike on MATLAB Graphics blog.
Mike is a computer graphics guru, well informed on the state of the art and knee-deep in the operational knowledge of... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/community/2014/10/03/welcome-mike/">read more >></a></p>http://blogs.mathworks.com/community/?p=2916Fri, 03 Oct 2014 14:48:55 +0000I am very excited to welcome a new blogger to MATLAB Central: Mike Garrity is going to be writing about graphics on the appropriately named Mike on MATLAB Graphics blog.

Mike is a computer graphics guru, well informed on the state of the art and knee-deep in the operational knowledge of How It All Works. He helped develop many of the features he’ll be describing, and he is a thoroughly entertaining writer. Beyond all this goodness, MATLAB graphics have taken a tremendous leap forward with the latest release, and Mike’s blog will give you an ideal perch to learn all about it. I recommend heading over there and seeing what he has to say.

]]>UncategorizedReading parameters from a file in MATLAB
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/qEsf5BtZA8Q/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2014/09/29/reading-parameters-from-a-file-in-matlab/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201411/523/62009828001_3877414682001_350-ParamReader-thumb.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 03:11</span>
</div>
</a></div><p>In last week’s post, I showed how to generate code in MATLAB to automatically read in a file. In this video I show how to make a simple wrapper function around code generated in that way. This wrapper will allow MATLAB to quickly look up a parameter in... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2014/09/29/reading-parameters-from-a-file-in-matlab/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1493Mon, 29 Sep 2014 17:11:39 +0000how to generate code in MATLAB to automatically read in a file. In this video I show how to make a simple wrapper function around code generated in that way. This wrapper will allow MATLAB to quickly look up a parameter in a text file. This architecture was picked because I need to have a non-MATLAB user modify parameters when the application is deployed on their site.