MATLAB Central BlogsPipes Output
http://pipes.yahoo.com/pipes/pipe.info?_id=e72c9f53f10ba5cad505dbbed8d501cb
Wed, 01 Jul 2015 06:18:43 +0000http://pipes.yahoo.com/pipes/Dubrulle Creates A Faster Tridiagonal QR Algorithm
http://feedproxy.google.com/~r/mathworks/moler/~3/ssISXGoyGPo/
<div class="overview-image"><img src="http://blogs.mathworks.com/cleve/files/feature_image/eigsvdgui.jpg" class="img-responsive attachment-post-thumbnail wp-post-image" alt="eigsvdgui"/></div><p>Augustin (Austin) Dubrulle deserves to be better known in the numerical linear algebra community. His version of the implicit QR algorithm for computing the eigenvalues of a symmetric tridiagonal matrix that was published in a half-page paper in <i>Numerische Mathematik</i> in 1970 is faster than Wilkinson's version published earlier. It is still a core algorithm in MATLAB today.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2015/06/29/dubrulle-creates-a-faster-tridiagonal-qr-algorithm/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1221Mon, 29 Jun 2015 17:00:34 +0000

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

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

Austin Dubrulle

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

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

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

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

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

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

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

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

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

References

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

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

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

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

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

Handbook Alters History

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

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

Impact Today

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

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

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

eigsvdgui

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

]]>Getting short path names
http://feedproxy.google.com/~r/mathworks/pick/~3/J-DpeaxEKp8/
<p>
Jiro's pick this week is Short Path Name on Windows by Jerome Briot.Thanks, Jerome! Your File Exchange entry helped me big time in my project. I was creating a program that managed a bunch of files, and some of the tasks required me to copy files from one network location... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/06/26/getting-short-path-names/">read more >></a></p>http://blogs.mathworks.com/pick/?p=6065Fri, 26 Jun 2015 13:00:02 +0000

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

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

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

]]>PicksHow Do You Modify the Background of an Image?
http://feedproxy.google.com/~r/mathworks/loren/~3/zjDQQIq1few/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/LorenColorVisaImageWhiteBGInterpolated.png"/></div><p>Today I'd like to introduce guest blogger <a rel="nofollow" target="_blank" href="http://www.mathworks.com/matlabcentral/profile/authors/845693-brett-shoelson">Brett Shoelson</a>. Some of you may know Brett through <a rel="nofollow" target="_blank" href="http://www.mathworks.com/matlabcentral/profile/authors/845693-brett-shoelson?utf8=%E2%9C%93&detail=fileexchange">his File Exchange submissions</a>, or through his involvement with the <a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/pick/">Pick of the Week</a> blog, or from occasional guest posts on <a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/steve/">Steve’s blog on image processing</a>.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/06/24/how-do-you-modify-the-background-of-an-image/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1201Wed, 24 Jun 2015 15:36:01 +0000

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

The original photo

First question: how do we isolate the background?

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

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

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

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

Is the effort of automation justified?

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

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

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

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

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

Improving the mask

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

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

Now the mask looks more appropriate for my task:

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

Modifying the background

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

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

Two problems remain...

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

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

Planewise manipulations

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

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

Fixing the interface

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

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

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

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

A final note

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

imgout = planewise(fcnhandle,rgbimg,varargin)

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

]]>Image ProcessingNew image batch processor app in R2015a
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/nKVshOpkcr4/
<div class="overview-image"><img src="http://blogs.mathworks.com/steve/files/imagebatchprocessor_spec_fcn.png" class="img-responsive attachment-post-thumbnail wp-post-image" alt="imagebatchprocessor_spec_fcn"/></div><p>Nine years ago I wrote a blog post showing how to do batch processing of image files. That is, I showed how to use MATLAB to perform the same operation on a bunch of image files in a particular folder.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2015/06/23/new-image-batch-processor-app-in-r2015a/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1338Tue, 23 Jun 2015 21:23:31 +0000

Nine years ago I wrote a blog post showing how to do batch processing of image files. That is, I showed how to use MATLAB to perform the same operation on a bunch of image files in a particular folder.

The basic procedure is:

Get a list of filenames.

Determine the processing steps to follow for each file.

Put everything together in a for loop.

The processing steps for each file typically looked like this:

Read in the data from the file.

Process the data.

Construct the output filename.

Write out the processed data.

I showed a sample processing loop that cropped and resized a bunch of images the same way:

files = dir('*.JPG');
for k = 1:numel(files)
rgb = imread(files(k).name);
rgb = rgb(1:1800, 520:2000, :);
rgb = imresize(rgb, 0.2, 'bicubic');
imwrite(rgb, ['cropped\' files(k).name]);
end

You can still do it exactly this way today, of course, in 2015. But you really might want to use the new Image Batch Processor App instead:

This is a very nice app that was just added to the Image Processing Toolbox in R2015a (released earlier this year, in March). If you look at the tool strip from left to right, you can see the entire workflow laid out for you.

Load the images. You can specify the folder containing your data.

What function do you want to run on each image? You can specify an existing function, or you can let the app make a shell function for you, and then you can fill in the details.

Where do you want to put the output? Here you can say you just want to overwrite the input files. (Be careful!)

Do you want to use a parallel cluster?

Do you want to process all the files, or just a selection?

Do you want to generate some code so you can automate the entire process?

]]>UncategorizedBus Object Bus Creator
http://feedproxy.google.com/~r/mathworks/pick/~3/kw8aY4AfqHg/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/pick/Richard/main_BusObjectBusCreator/main_BusObjectBusCreator_04.png"/></div><p>
Posted by Richard Ruff , June 19, 2015Richard is an Application Engineer at MathWorks focused on the Embedded Coder product for code generation, primarily in the Aerospace industry.Richard's pick this week is Bus Object Bus Creator by Landon Wagner.PickMy pick this week is the Simulink library submission for automatically populating... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/06/19/bus-object-bus-creator/">read more >></a></p>http://blogs.mathworks.com/pick/?p=6056Fri, 19 Jun 2015 13:00:19 +0000

Posted by Richard Ruff , June 19, 2015

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

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

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

The description from the entry sums up the benefits nicely:

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

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

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

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

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

Comments

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

]]>PicksGetting Started with Kaggle Data Science Competitions
http://feedproxy.google.com/~r/mathworks/loren/~3/28aLOvKGRvw/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/Kaggle_Titanic_05.png"/></div><p>Have you been interested in data science competitions, but not sure where to begin? Today's guest blogger, Toshi Takeuchi, would like to give a quick tutorial on how to get started with Kaggle using MATLAB.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/06/18/getting-started-with-kaggle-data-science-competitions/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1195Thu, 18 Jun 2015 18:44:57 +0000

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

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

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

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

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

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

train.csv - the training data

test.csv - the test data

Data Import and Preview

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

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

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

1 - Survived
0 - Didn't survive

Establishing the Baseline

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

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

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

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

Back to Examining the Data

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

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

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

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

avgAge =
29.699

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

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

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

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

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

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

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

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

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

Feature Engineering

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

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

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

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

Your Secret Weapon - Classification Learner

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

classificationLearner

Click on Import Data

Select Train in Step 1 in Set Up Classification dialog box

In Step 2, change the "Import as" value for PassengerId to "Do not import", and Survived to "Response". All other variables should be already marked as Predictor.

In Step 3, just leave it as is to Cross Validation.

Random Forest and Boosted Trees

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

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

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

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

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

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

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

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

oobAccuracy =
0.82212

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

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

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

Model Evaluation

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

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

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

Create a Submission File

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

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

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

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

]]>UncategorizedRegion filling and Laplace’s equation
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/fh8TjSFI_6M/
<div class="overview-image"><img src="http://blogs.mathworks.com/steve/files/exploring_regionfill_06.png" class="img-responsive attachment-post-thumbnail wp-post-image" alt="Region to be filled"/></div><p>Today I want to show you how to use a linear system of 90,300 equations and 90,300 unknowns to get rid of a leaf.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2015/06/17/region-filling-and-laplaces-equation/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1333Wed, 17 Jun 2015 15:44:25 +0000

Today I want to show you how to use a linear system of 90,300 equations and 90,300 unknowns to get rid of a leaf.

Earlier this month I told you how the function roifill was recently renamed to regionfill to correct a functional design flaw. Now I want to follow up and talk about the algorithm underlying regionfill.

This was one of my earliest algorithm projects at MathWorks. And it was my earliest lesson the power of sparse linear algebra in MATLAB. About 21 years ago, Clay Thompson, one of the two developers of version 1.0 of Image Processing Toolbox, commented to me one day that maybe we could "erase" objects in an image by treating it as a "soap-film problem." (The other developer of that first toolbox version was Loren Shure.)

In other words, let's treat the pixel values surrounding the region we want to erase as heights of a closed loop of wire and then figure out the shape of a soap film bounded by that wire.

That sounded plausible, so I looked into it. (The World Wide Web was barely a year old at the time, and there was nothing like the modern search engine, so I sought out information using what we used to call "books.")

I found that, if you make some simplifying assumptions (such as no gravity), the soap film surface satisfies Laplace's equation:

where the height of the wire loop around the region gives us the boundary conditions for the PDE.

In a simple, discretized version of Laplace's equation, the value of every grid element in the interior of the region equals the average of its north, east, south, and west neighbors in the grid.

Before I explore that idea further, though, let's look at some pictures to illustrate what we're trying to accomplish.

Here is the "trees" image (in gray) with a region drawn around what looks like a small leaf.

[X,map] = imread('trees.tif');
I = im2double(ind2gray(X,map));
imshow(I)
x = [240 245 255 270 285 290 280 270 250 240];
y = [155 165 170 173 172 166 161 155 153 155];
hold on
plot(x,y,'y')
hold off

Here's a zoomed-in view.

axis([230 300 140 180])

We are going to "erase" that leaf by filling in the pixels inside that region. The function poly2mask is useful for finding the set of pixels inside a polygonal region.

mask = poly2mask(x,y,258,350);
imshow(mask)

Let's think about the gray-scale pixel values as heights and view the image as a surface.

h = surf(I);
h.EdgeColor = 'none';
h.FaceColor = 'interp';
axis ij
view(-10,80)

Here's a close-up of the region we're interested in.

xlim([230 300])
ylim([140 180])

We can use the region mask computed above to "cut out" the part of the surface inside the region.

h.FaceAlpha = 'interp';
h.AlphaData = ~mask;

A very simple approach to erasing the left would be to fill in the region with the mean value of pixels in the region. Here's how to do it. (I just love logical indexing with binary image masks!)

J = I;
J(mask) = mean(I(mask));
h.ZData = J;
h.AlphaData = true(size(J));
title('Region filled with mean value')

imshow(J)
title('Region filled with mean value')

axis([230 300 140 180])

Now let's work toward a solution based on Laplace's equation. I'm going to set up a linear system of equations, $Ax = b$, so that the solution gives me every pixel of the output image. I'll think of the pixels as being numbered in columnwise order from 1 to the total number of pixels. The trees image has 90,300 pixels, so this is going to be a 90,300 by 90,300 linear system! It's a good thing it'll be very sparse.

For every pixel outside the masked region, the equation is simply the identity: output_pixel equals input_pixel.

For every pixel inside the masked region, the pixel value should equal the average of the pixel values of its north, east, south, and west neighbors.

First, number the pixels inside the masked region.

u = find(mask);

Now number the pixels outside the region.

w = find(~mask);

Using a technique I called "neighbor indexing" years ago, we can find the north, east, and south, and west neighbors for all the pixels in the masked region this way:

M = size(mask,1);
u_north = u - 1;
u_east = u + M;
u_south = u + 1;
u_west = u - M;

(Note: for the purpose of this blog post, I'm omitting code needed to handle the case where some mask pixels lie on the border of the image.)

Next, construct a sparse matrix representing the linear system.

v = ones(size(u));
ijv_mask = [...
u u 1.00*v
u u_north -0.25*v
u u_east -0.25*v
u u_south -0.25*v
u u_west -0.25*v ];
ijv_nonmask = [w w 1.00*ones(size(w))];
ijv = [ijv_mask; ijv_nonmask];
A = sparse(ijv(:,1),ijv(:,2),ijv(:,3));

What does the sparse linear system look like?

spy(A)

The right-hand vector, b, contains either the original pixel values (for pixels outside the mask) or 0 (for pixels inside the mask).

Today, the process of filling image regions is often called "inpainting." Many different inpainting algorithms have been published.

I'll finish by confessing that I wrote an iterative solver for my first MATLAB implementation of this technique. When I showed my work to Cleve, he said, "Why don't you use sparse backslash?"

]]>UncategorizedALGOL 60, PL/0 and MATLAB
http://feedproxy.google.com/~r/mathworks/moler/~3/Fuudu6r_EEw/
<div class="overview-image"><img src="http://blogs.mathworks.com/cleve/files/feature_image/pl0.jpeg" class="img-responsive attachment-post-thumbnail wp-post-image" alt="pl0"/></div><p>The 1960 programming language ALGOL 60 influenced the design of many subsequent languages, including C and a miniature language called PL/0. I based the design of the first MATLAB on PL/0.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2015/06/15/algol-60-pl0-and-matlab/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1218Mon, 15 Jun 2015 16:01:22 +0000

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

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

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

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

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

JOVIAL

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

ALGOL 60

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

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

ACM Algorithms

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

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

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

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

Example

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

Numerische Mathematik

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

EISPACK

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

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

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

LINPACK

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

Niklaus Wirth

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

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

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

PL/0

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

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

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

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

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

]]>A Simple Backup Utility
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/j9lo6tp-scQ/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2015/06/15/a-simple-backup-utility/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201504/711/62009828001_4161185672001_simple-backup-thumbnail4.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 04:16</span>
</div>
</a></div><p>I have a number of scripts that I have been maintaining for many years. They are not in version control but I still want to make edits and be able to access the code in previous versions. So I use a simple utility that saves a version of the file... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2015/06/15/a-simple-backup-utility/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1663Mon, 15 Jun 2015 14:57:42 +0000I have a number of scripts that I have been maintaining for many years. They are not in version control but I still want to make edits and be able to access the code in previous versions. So I use a simple utility that saves a version of the file to a backup folder before I make major edits.

For more information about programmatically accessing the MATLAB Editor type:

>>help matlab.desktop.editor

]]>UncategorizedSpeeding up Initialization for Accelerator Mode
http://feedproxy.google.com/~r/SethOnSimulink/~3/Lbym7DpFHoM/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/seth/2015Q2/simNoChecksum.png"/></div><p>Some time ago, I wrote a post about Getting the most out of Rapid Accelerator mode. That post describes how to use the RapidAcceleratorUpToDateCheck = 'off' option to skip the initialization time.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/06/12/speeding-up-initialization-for-accelerator-mode/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4621Fri, 12 Jun 2015 16:48:20 +0000Some time ago, I wrote a post about Getting the most out of Rapid Accelerator mode. That post describes how to use the RapidAcceleratorUpToDateCheck = 'off' option to skip the initialization time.

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

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

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

The Problem

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

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

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

Let see how we can improve that.

The Solution: Accelerated Model Reference

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

When doing the same test as before, we get:

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

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

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

And run our test to measure the initialization time:

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

Now it's your turn

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

]]>Introduction to MuPAD® for Plant Modeling
http://feedproxy.google.com/~r/mathworks/pick/~3/7Bdv_M7NkVA/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/pick/jiro/potw_mupadintro/potw_mupaddocexample.png"/></div><p>
Jiro's pick this week is Introduction to MuPAD for Plant Modeling by JMAAB.
If you are affiliated with the automotive industry, you may have heard of MathWorks Automotive Advisory Board (MAAB). JMAAB (Japan MBD Automotive Advisory Board) is the Japanese counterpart.
... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/06/12/introduction-to-mupad-for-plant-modeling/">read more >></a></p>http://blogs.mathworks.com/pick/?p=6047Fri, 12 Jun 2015 13:00:53 +0000

This introductory document on MuPAD was originally published about 2 years ago by JMAAB in Japanese. It was a part of an initiative
for studying the workflow of using MuPAD for plant modeling. In fact, the overarching motivation was to promote the use of
formula manipulation techniques in the field of automotive control system development.

Why use MuPAD for plant modeling? There are a variety of MathWorks tools for doing plant modeling. One of them is Simscape, which allows you to create plant models using blocks that represent physical components.

However, there are cases where you may need to formulate and derive the equations that govern the physical system based on
first principles. MuPAD can greatly help in this area. This document introduces the MuPAD language and environment from the
point of view of plant modeling. It also shows how to output the MuPAD results in the Simscape language. This File Exchange
entry is an English translation of the original document.

While the initial target audience were people needing formula manipulation tools for plant modeling, this can be a nice introductory
tutorial for anyone wanting to use MuPAD. MuPAD is a full language in and of itself, so becoming completely proficient may
take some time. However, this document uses a couple of simple examples to demonstrate how one might go about deriving mathematical
formulations. The examples are shown in a two-column format with additional information on the right column.

In addition, there are some useful appendices discussing concepts that are somewhat confusing to MATLAB users, such as "Concept
of Assignment (:=) and Equal(=)" and "Differences and Relationship Between MuPAD and MATLAB".

Comments

Do you perform symbolic calculations as part of creating plant models? Let us know here or leave a comment for JMAAB.

]]>Advice for Making Prettier Plots
http://feedproxy.google.com/~r/mathworks/loren/~3/KImEYoXR2Cg/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/prettierPlots_05.png"/></div><p>A few years ago, <a rel="nofollow" target="_blank" href="http://www.mathworks.com/matlabcentral/profile/authors/869871-jiro-doke">Jiro</a> wrote a <a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/loren/2007/12/11/making-pretty-graphs/">popular post</a> for making pretty plots on this blog. We also host a <a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/graphics">blog specifically about graphics</a> by Mike. And with the R2014b release of MATLAB came an updated graphics system that Dave described last year in a 3 part series: <a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/loren/2014/10/03/matlab-r2014b-graphics-part-1-features-of-the-new-graphics-system/">Part 1</a>, <a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/loren/2014/10/14/matlab-r2014b-graphics-part-2-using-graphics-objects/">Part 2</a>, and <a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/loren/2014/11/05/matlab-r2014b-graphics-part-3-compatibility-considerations-in-the-new-graphics-system/">Part 3</a>.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/06/11/advice-for-making-prettier-plots/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1182Thu, 11 Jun 2015 13:27:29 +0000

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

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

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

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

Update the Y Tick Labels with Euros

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

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

Finally Annotate a Plot with Some Math

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

]]>UncategorizedXKCDIFY
http://feedproxy.google.com/~r/mathworks/pick/~3/3AmlFCRBAyI/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/pick/Sean/mainxkcdify/mainxkcdify_03.png"/></div><p>
Sean's pick this week is xkcdify by Stuart Layton.
MATLAB has lots of stock plotting routines and many more that come with the various toolboxes. You can see the... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/06/05/xkcdify/">read more >></a></p>http://blogs.mathworks.com/pick/?p=6044Fri, 05 Jun 2015 13:00:24 +0000

MATLAB has lots of stock plotting routines and many more that come with the various toolboxes. You can see the ones available
to you on the plots tab:

However, all of these plots are very formal. MATLAB doesn't have any stock options for informal plotting. The comic strip
xkcd does a very good job communicating data in an informal way and Stuart has given us xkcdify to take our existing MATLAB plots and xkcdify or informalize them.

We've had some unseasonably cold weather in Massachusetts over the last few days (surprise I know!). Let's take a look at
it with a control chart before and after xkcdification. I'll use the temperatures that I'm comfortable in a t-shirt and shorts
for as the limits.

If this isn't informal enough, we can re-xkcdify too!

xkcdify(gca)

There are a couple changes I made to xkcdify to make it work on more plot types.

To separate a line into multiple pieces, you can use NaN. For the length of line calculation, I removed these NaNs because
otherwise the line length was always NaN.

% Line 209:210
goodIdx = ~isnan(x) & ~isnan(y) & isfinite(x) & isfinite(y);
len = sum(hypot(diff(x(goodIdx)),diff(y(goodIdx))));

There are cases where a line can be empty so I put an isempty() check around this.

% Line 129:132if ~isempty(x)
x = x + smooth( generateNoise(n) .* rand(n,1) .* jx )';
y = y + smooth( generateNoise(n) .* rand(n,1) .* jy )';
end

Comments

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

]]>PicksImage region filling – an updated design
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/KenK1CEiMs8/
<div class="overview-image"><img src="http://blogs.mathworks.com/steve/files/roifill-mask-2.png" class="img-responsive attachment-post-thumbnail wp-post-image" alt="roifill-mask-2"/></div><p>About a year ago, I wrote a blog post criticizing one of my functional designs from the 1990s, roifill. Now I have an update based on the R2015a release of Image Processing Toolbox.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2015/06/04/image-region-filling-an-updated-design/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1326Thu, 04 Jun 2015 20:01:24 +0000

About a year ago, I wrote a blog post criticizing one of my functional designs from the 1990s, roifill. Now I have an update based on the R2015a release of Image Processing Toolbox.

Let me show you an example to refresh your memory about roifill.

I = imread('rice.png');
imshow(I)

You can use roifill to fill in pixels inside a region. (ROI in the function name stands for ''region of interest.'') You can specify the region interactively, using a mouse, or by providing the coordinates of a polygonal boundary, or by providing a binary mask. Here I'll provide the coordinates directly.

As I described in my earlier post, the form of roifill that takes a mask image suffered from a design flaw. Most people would reasonably assume that the mask image would specify the set of pixel values that are to be filled. And that's the way I should have designed it. Unfortunately, though, roifill actually only fills the interior pixels of the mask.

Suppose your fill mask looked like this:

The function roifill only replaces these interior pixels of the mask:

The pixels on the edge of the mask are used to establish boundary conditions for the fill equation. Those pixels stay the same, which is not what most people expect.

So how could we fix this problem? Changing the default behavior of roifill would be a significant incompatibility. But we have learned that optional behaviors often go undiscovered, so most people don't benefit from them.

The Image Processing Toolbox team decided to address the problem by introducing a new function, regionfill, with the desired behavior. The function roifill remains in the product, so existing code that uses it will continue to work. Presumably, anyone currently using roifill with the mask syntax has already worked around its awkward behavior.

At the bottom of the reference page for regionfill, there is a note that this function was introduced in R2015a.

And at the top of the reference page for roifill, there is a note encouraging the use of regionfill.

]]>UncategorizedFriedrich Bauer
http://feedproxy.google.com/~r/mathworks/moler/~3/C56Ar8IEdBI/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/cleve/fritzbauer.jpg"/></div><p>Fritz Bauer, eminent German computer scientist and last surviving member of the organizing committee of the 1964 Gatlingburg Conference on Numerical Algebra, passed away on March 26 at the age of 90.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2015/06/01/friedrich-bauer/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1208Mon, 01 Jun 2015 17:00:42 +0000

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

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

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

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

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

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

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

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

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

Stack

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

Computer Science and Software Engineering

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

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

This definition of software engineering is now universally quoted.

Algol

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

Community

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

Numerical Analysis

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

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

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

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

is the condition number of the eigenvector matrix.

There is a proof in Wikipedia.

Bavarian Alpinist

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

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

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

]]>Setting Initial MATLAB Working Folder to Last Folder Used
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/Nehls8_ES4A/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2015/06/01/setting-initial-matlab-working-folder-to-last-folder-used/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201504/320/62009828001_4162901277001_initial-working-folder-thumbnail7b.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 01:43</span>
</div>
</a></div><p>In release 2014b there is a new preference that lets you set MATLAB’s initial working folder to be the last folder used in the previous MATLAB session. Here I show when it can be useful and when you can’t use it.
... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2015/06/01/setting-initial-matlab-working-folder-to-last-folder-used/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1642Mon, 01 Jun 2015 15:42:33 +0000In release 2014b there is a new preference that lets you set MATLAB’s initial working folder to be the last folder used in the previous MATLAB session. Here I show when it can be useful and when you can’t use it.

]]>Format: VideoimageSet Viewer
http://feedproxy.google.com/~r/mathworks/pick/~3/ugqDuHrLbRU/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/pick/files/AviPOWT1_01-1024x514.png"/></div><p>Avi Nehemiah is the product marketing manager for computer vision applications.
Avi's pick of the week is imageSet viewer by Brett Shoelson.
My pick of the week is the imageSet Viewer user interface created by my friend Brett Shoelson. The biggest challenges I face when working with large sets of images(that are... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/05/29/imageset-viewer/">read more >></a></p>http://blogs.mathworks.com/pick/?p=6012Fri, 29 May 2015 13:00:05 +0000Avi Nehemiah is the product marketing manager for computer vision applications.
Avi's pick of the week is imageSet viewer by Brett Shoelson.
My pick of the week is the imageSet Viewer user interface created by my friend Brett Shoelson. The biggest challenges I face when working with large sets of images(that are normally stored in several different directories) is it was difficult to bring the data into MATLAB, and interactively visualize the images.
The imageSet functionality in the Computer Vision System Toolbox helps me bring the data into MATLAB. I can bring all the images in a folder into MATLAB, and maintain the hierarchical relationship between folders using the imageSet. The image data I am using is from Caltech 101, collected by Fei-Fei Li, Marco Andreetto, and Marc 'Aurelio Ranzato. imageSet functionality is new in R2014b, you can learn more about it in the documentation .

imageData = imageSet('Data','recursive');

Now let's take a look at some of properties of the imageData variable I just created.

This is where the imageSet Viewer comes in to help me interactively visualize this data

imageSetViewer(imageData)

Other options to navigate through your images using the user interface are

Click on the tabs at the top to view different folders or categories

Click on an image to view in a separate window

Right click on an image to save to workspace

Automatically creates imageSet if you pass in the directory name as an argument

Click on the name of an image to copy to clipboard

I'd recommend trying out the imageSet viewer anytime you are working with sets of images.
]]>PicksCreating Custom Valve in Simscape
http://feedproxy.google.com/~r/SethOnSimulink/~3/AapfwQxUJ4U/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/seth/2015Q2/pilotValveTimeResults.png"/></div><p>Today I want to share a technique I like to use when I need to create custom hydraulic components in Simscape.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/05/26/creating-custom-valve-in-simscape/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4585Wed, 27 May 2015 02:31:22 +0000Today I want to share a technique I like to use when I need to create custom hydraulic components in Simscape.

The Problem

I often get questions similar to the following from users:

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

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

Overview of the solution

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

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

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

Getting the Data

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

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

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

This gives us the data we need for our model.

Option 1: Using blocks

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

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

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

Option 2: Custom Simscape Component

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

The code looks like the following:

The Results

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

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

Now it's your turn

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

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

]]>Including External Code in Published Document
http://feedproxy.google.com/~r/mathworks/loren/~3/lwIyA3mqul4/
<p>When I wanted to show you a code snippet in the past, I displayed the code in the text of my blog post.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/05/22/including-external-code-in-published-document/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1176Fri, 22 May 2015 16:40:22 +0000

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

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

How I Do It Now

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

<include>fileOfInterest.m</include>

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

Let's Try It Out!

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

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

]]>Writing A Cell Array – A Comparison
http://feedproxy.google.com/~r/mathworks/pick/~3/Gp55T01hGLk/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/pick/Sean/maincell2csv/cell2csv.png"/></div><p>
Sean's pick this week is a comparison of a few cell-writing functions.
Contents
... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/05/22/writing-a-cell-array-a-comparison/">read more >></a></p>http://blogs.mathworks.com/pick/?p=6002Fri, 22 May 2015 13:00:13 +0000

Sean's pick this week is a comparison of a few cell-writing functions.

There are many high-level functions for writing common datatypes to text file formats such as writetable, csvwrite, and dlmwrite. However, you might have noticed that none of these work on cell arrays. Since cell arrays can and often do contain highly
heterogenous data, you'll often need to use lower-level functionality like fprintf to write them out given the knowledge of what the cell contains. My pick today is a comparison of a few files that have
tried to automate this while trying to do the "write" thing.

Selected Entries

Below, I've selected a few files to compare based on a quick search and the fact that they're all covered under the BSD license:

fwritecell was unable to resolve this cell array because it needs consistent formatting. Because it relies on cell2mat, this continues to pop up. The others pass.

Test 3: Speed

This test is for speed. I'm going to write a large cell array to see what happens.

C = repmat(C,1000,10);
disp('Size of C')
disp(size(C))

The winner here is going to be impcell2csv with dlmcell coming in a close second. impcell2csv wrote NA in the spot of the nested cell, dlmcell left it blank. The others gave non-descript errors with the file partially written.

Challenge Does anyone know of a file that will handle this? Perhaps padding with extra commas or providing other options?

Test 5: N-Dimensional Arrays

Here we will have a three-d cell array.

C = {pi, 'Sean', uint8(42); uint8(42), exp(1), exp(45); 'MaTLaB','Test 2', 1; -1.23 -exp(45) 19};
C = cat(3,C,C);
disp(C)

impcell2csv and cellwrite did the same thing here, append the second page to the right of the first page. I guess that's a tie, cell2csv and dlmcell both ignored the second page altogether.

Test 6: Help/Flexibilty

How's the help/flexibility?

help dlmcell

<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> %%
<><><><><> dlmcell - Write Cell Array to Text File <><><><><> %
<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> %
Version: 01.06.2010 %
(c) Roland Pfister %
roland_pfister@t-online.de %
...with many thanks to George Papazafeiropoulos %
for his corrections and improvements. %
1. Synopsis %
%
A single cell array is written to an output file. Cells may consist of %
any combination of (a) numbers, (b) letters, or (c) words. The inputs %
are as follows: %
%
- file The output filename (string). %
- cell_array The cell array to be written. %
- delimiter Delimiter symbol, e.g. ',' (optional; %
default: tab ('\t'}). %
- append '-a' for appending the content to the %
output file (optional). %
%
2. Example %
%
mycell = {'Numbers', 'Letters', 'Words','More Words'; ... %
1, 'A', 'Apple', {'Apricot'}; ... %
2, 'B', 'Banana', {'Blueberry'}; ... %
3, 'C', 'Cherry', {'Cranberry'}; }; %
dlmcell('mytext.txt',mycell); %
%
<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> %

help fwritecell

FWRITECELL writes formatted data from a cell array to a text file.
fwritecell(filename, format, data)
writes data using the C language conversion specifications in
variable "format" as with SPRINTF
Example: fwritecell('textfile1.txt','%2d %1d %21s %8.5f',C);
fwritecell(filename, data)
writes data using a fixed width format padded with whitespace. Note
that the decimal point for floating point numbers may not be aligned

help cell2csv

Writes cell array content into a *.csv file.
CELL2CSV(fileName, cellArray, separator, excelYear, decimal)
fileName = Name of the file to save. [ i.e. 'text.csv' ]
cellArray = Name of the Cell Array where the data is in
separator = sign separating the values (default = ';')
excelYear = depending on the Excel version, the cells are put into
quotes before they are written to the file. The separator
is set to semicolon (;)
decimal = defines the decimal separator (default = '.')
by Sylvain Fiedler, KA, 2004
updated by Sylvain Fiedler, Metz, 06
fixed the logical-bug, Kaiserslautern, 06/2008, S.Fiedler
added the choice of decimal separator, 11/2010, S.Fiedler

help impcell2csv

% Writes cell array content into a *.csv file.
%
% CELL2CSV(fileName, cellArray[, separator, excelYear, decimal])
%
% fileName = Name of the file to save. [ e.g. 'text.csv' ]
% cellArray = Name of the Cell Array where the data is in
%
% optional:
% separator = sign separating the values (default = ',')
% excelYear = depending on the Excel version, the cells are put into
% quotes before they are written to the file. The separator
% is set to semicolon (;) (default = 1997 which does not change separator to semicolon ;)
% decimal = defines the decimal separator (default = '.')
%
% by Sylvain Fiedler, KA, 2004
% updated by Sylvain Fiedler, Metz, 06
% fixed the logical-bug, Kaiserslautern, 06/2008, S.Fiedler
% added the choice of decimal separator, 11/2010, S.Fiedler
% modfiedy and optimized by Jerry Zhu, June, 2014, jerryzhujian9@gmail.com
% now works with empty cells, numeric, char, string, row vector, and logical cells.
% row vector such as [1 2 3] will be separated by two spaces, that is "1 2 3"
% One array can contain all of them, but only one value per cell.
% 2x times faster than Sylvain's codes (8.8s vs. 17.2s):
% tic;C={'te','tm';5,[1,2];true,{}};C=repmat(C,[10000,1]);cell2csv([datestr(now,'MMSS') '.csv'],C);toc;

help cellwrite

CELLWRITE Write a cell array to a comma separated value file.
CELLWRITE(FILENAME, C) writes cell array C into FILENAME as comma
separated values.
NOTE: This function is not completely compatible with CSVWRITE.
Offsets are not supported and 0 values are not omitted.
See also CSVWRITE, CSVREAD, DLMREAD, DLMWRITE, WK1READ, WK1WRITE.

The help for all of them is good. I like dlmcell's ability to append, and I like impcell2csv's description of the enhancements though the help is mostly borrowed from cell2csv.

Comments

Do you use cell arrays to store data? If so, do you need to write them out to text often? Have you ever written your own
fprintf wrapper or used one of the above files or another that I didn't see? Let us know here or leave a comment for one of the authors above on their File Exchange entry page.

]]>PicksWhat is the most efficient aircraft seating strategy?
http://feedproxy.google.com/~r/SethOnSimulink/~3/5-tpCzy6l10/
<div class="overview-image"><img src="http://blogs.mathworks.com/seth/files/feature_image/EntityGeneration.png" class="img-responsive attachment-post-thumbnail wp-post-image" alt="EntityGeneration"/></div><p>Today I am happy to welcome my colleague Ramamurthy Mani for another very interesting study using SimEvents.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/05/20/what-is-the-most-efficient-aircraft-seating-strategy/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4580Wed, 20 May 2015 15:14:34 +0000Today I am happy to welcome my colleague Ramamurthy Mani for another very interesting study using SimEvents.

Introduction

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

Generating a queue of passengers

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

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

The six strategies include:

frontToBackZone: Boarding from front of aircraft to back broken into 4 zones.

backToFrontZone: Boarding from back of aircraft to front broken into 4 zones.

random: Boarding happens across the aircraft randomly

alternateRowZone: Boarding (back to front) is done in four zones and each zone only includes alternate rows.

windowToAisleZone: Boarding broken into 3 zones as ‘window’, ‘middle’, and ‘aisle’

windowToAisleAlternateZone: Hybrid of 4 and 5 with 6 zones where alternate rows in window, middle, and aisle are treated as separate zones.

Run simulations

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

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

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

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

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

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

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

Results

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

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

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

Now it's your turn

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

]]>The Ziggurat Random Normal Generator
http://feedproxy.google.com/~r/mathworks/moler/~3/zg2PPlDHXJk/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/cleve/ziggurat_blog_01.png"/></div><p>This is the third in a multi-part series on the MATLAB random number generators. MATLAB has used variants of George Marsaglia's ziggurat algorithm to generate normally distributed random numbers for almost twenty years.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2015/05/18/the-ziggurat-random-normal-generator/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1201Mon, 18 May 2015 17:00:42 +0000

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

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

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

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

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

z = zigplot(6)

z =
0
0.8288
1.1713
1.4696
1.7819
2.1761

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

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

Initialization

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

Central algorithm

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

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

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

Accuracy

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

Variations

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

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

Underlying Uniform Generator

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

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

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

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

m =
25000000

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

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

tu =
0.3416
tn =
0.4520
ratio =
0.7558

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

Acknowledgements

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

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

References

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

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

]]>Control Excel using ActiveX
http://feedproxy.google.com/~r/mathworks/pick/~3/HIuylMnAWHQ/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/pick/jiro/potw_xcollection/xcollection_methods.png"/></div><p>
Jiro's pick this week is The X Collection by Yvan Lengwiler.
My day job is a customer training engineer. I get to take people who have never used MATLAB or Simulink to the point where
they feel... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/05/15/control-excel-using-activex/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5990Fri, 15 May 2015 13:00:59 +0000

My day job is a customer training engineer. I get to take people who have never used MATLAB or Simulink to the point where
they feel comfortable in applying the tools in their day-to-day work. It's certainly a rewarding role.

Recently during my MATLAB Fundamentals training course, I received a few questions regarding how to exchange data between
MATLAB and Microsoft® Excel®. In the training, we talk about functions like xlsread and xlswrite. The person wanted the ability to control the formatting of cells, which xlswrite does not provide. Specifically, he wanted to change the width of a column and some of the cell colors. I pointed him to the
function actxserver, which provides a way to connect to Windows® applications using ActiveX. Once connected to Excel, he can use the various APIs that are provided by Excel to customize. The reference page for actxserver has an example for connecting to Excel.

I started browsing the File Exchange in hopes for finding additional examples that I could point people to for using actxserver with Excel. Sure enough, I found quite a few entries, including this one (and many others) by our very own Brett Shoelson. I chose to highlight Yvan's X Collection, because he provides a suite of functions for performing various simple tasks within
Excel, using ActiveX. He uses intuitive function names, such as XConnect, XOpenBook, and XAddSheet. Each of these functions take in various handles (Excel application, workbook, or worksheet) as inputs, so it's easy to understand
how the functions are related to each other. In fact, if you open the files in the MATLAB Editor, you can place the cursor
at one of the handle variables to see how they are used or created. The editor will highlight the same variables, so you can
easily pick out where they are used.

One suggestion I have for improvement is to convert this into a class, instead of a suite of functions. Naturally, these functions
operate on handles that represent the Excel application, workbooks, and worksheets. Each of them can be a class, and each
class can contain methods that specifically operate on the instances of the class. This will make the functions inherently
invisible to everyone except to the class instances. You wouldn't have to catch errors resulting from invalid handle inputs.
Of course, you may still need to catch other errors that may occur. The other reason for creating a class in this case would
be to hide (encapsulate) the actual handles. The hExcel returned from XConnect is a COM.Excel_Application.

hExcel = XConnect()

hExcel =
COM.Excel_Application

This handle exposes the user to various methods and properties provided by the Excel application.

This may be good or bad, but perhaps the intent is for users to use the functions provided by this suite and not the ones
that are natively available. By creating a class that maintains the handles internally as a property, it protects the users
from misusing the handle.

Nonetheless, this is a great suite of functions for connecting MATLAB with Excel. I'm not doing justice to Yvan's entry if
I don't mention the true motivation for him creating this. If you need to repeatedly read from or write to an Excel file,
xlsread or xlswrite would suffer a performance hit because of the repeated connect/disconnect that it's doing with Excel. With Yvan's XCollection,
you connect to Excel once, do the necessary read/write operations, and disconnect once.

Note that Yvan's entry does not provide the functionality the my training customer was asking for, but I pointed him to this
entry as a starting point.

Comments

Yvan provides detailed documentation of how to use the functions. Give this a try and let us know what you think here or leave a comment for Yvan.

]]>New Video Player in Blogs Area
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/kChMGW4reBU/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2015/05/14/new-video-player-in-blogs-area/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201504/2512/62009828001_4162802751001_new-player-thumbnail10b.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 01:46</span>
</div>
</a></div><p>We have made some changes to how videos play in the Blogs area, including a new video player. There are one or two features that I’d like to make you aware of including:
The player works on devices that do not have Flash such as iPad, and most other mobile devices
The player... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2015/05/14/new-video-player-in-blogs-area/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1617Thu, 14 May 2015 20:32:47 +0000We have made some changes to how videos play in the Blogs area, including a new video player. There are one or two features that I’d like to make you aware of including:

The player works on devices that do not have Flash such as iPad, and most other mobile devices

The player can go full screen, which is useful for me as I will often be capturing a larger area of the screen in my videos than Doug did

The video auto-plays if you click on the thumbnails on the blog index page

]]>Format: VideoHow do you create a mask from a variable-thickness open, freehand curve?
http://feedproxy.google.com/~r/mathworks/pick/~3/AxM9HOuqTpo/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/pick/files/normalLine4.png"/></div><p>
Brett's Pick this week is 2D Line Curvature and Normals, by Dirk-Jan Kroon.Recently, I began working on an app to facilitate interactive image segmentations. Ultimately, I want to create a suite of tools that allow me to use some combination of automated and manual manipulations to create flexibly a complicated... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/05/08/how-do-you-create-a-mask-from-a-variable-thickness-open-freehand-curve/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5965Fri, 08 May 2015 13:00:09 +0000

Recently, I began working on an app to facilitate interactive image segmentations. Ultimately, I want to create a suite of tools that allow me to use some combination of automated and manual manipulations to create flexibly a complicated mask of my images. As part of that project, I want the ability to create closed or open (!)imfreehand or impoly regions, and to use those regions to modify my mask.

Creating masks from closed regions is easy with the imroi tools of the Image Processing Toolbox; the "createMask" method trivializes the process. But creating masks from open regions can produce some unexpected, and unwanted, results.

Suppose, for example, that I wanted to create a mask from a manual tracing of the lower left petal of the "yellow lily" image that ships with the Image Processing Toolbox:

I read in and display the image, and trace the boundary of interest:

lily = imread('C:\Program Files\MATLAB\R2015a\toolbox\images\imdata\yellowlily.jpg');
ax(1) = subplot(1,2,1);
imshow(lily)
title('Lovely Lily')
% I did this manually, saved the positions to a Mat-file:% h = impoly(imgca,'closed',false);
load outlinedLily
h = impoly(imgca,xsys,'closed',false);
xsys = getPosition(h);
xs = xsys(:,1);
ys = xsys(:,2);
% (Note: you can do this with |imfreehand|, too, but I find that I can be% much more accurate with |impoly|.)
defaultMask = createMask(h);
hold on
plot(xs,ys,'r.','markersize',16);

If I simply use the createMask method, I get a mask of the filled region within the closed ROI:

But that's not what I wanted! I could use bwperim (for example) on that mask, but then I would have to remove the effects of the unwanted region closing. What I really wanted was a thin "shell" around my drawn region, from which I could create my mask.

Enter Dirk's LineNormals2D.

N = LineNormals2D(xsys); %Dirk-Jan!!!! This rocks!
thicknessMultiplier = 2;
posn = [xs-thicknessMultiplier*N(:,1) ys-thicknessMultiplier*N(:,2);
flipud(xs+thicknessMultiplier*N(:,1)) flipud(ys+thicknessMultiplier*N(:,2))];
h = impoly(ax(2),posn);
xsys = getPosition(h);
xs = xsys(:,1);
ys = xsys(:,2);
plot(ax(1),xs,ys,'g.','markersize',16);

Now...

imshow(createMask(h),'parent',ax(2))
title('That''s more like it!')

The line above, in which I used the output of LineNormals2D to create my "posn" vector, might seem a bit puzzling. In a nutshell, Dirk-Jan's function calculated the position of the normal at each point on my imroi. I subtract those values from, and then add them to, the points on my curve to calculate desired positions. (thicknessMultiplier just scales the amount of offset.) I use flipud to have the "inward points" flip around and start where the "outward points" left off:

t = 0:pi/64:3*pi;
xy = [t',sin(t)'];
N = LineNormals2D(xy);
plot(xy(:,1),xy(:,2),'b.')
plot(xy(:,1)-N(:,1),xy(:,2)-N(:,2),'r.')
plot(xy(:,1)+N(:,1),xy(:,2)+N(:,2),'g.')
legend({'Original points','XY-N','XY+N'})
xlabel('t')
ylabel('sin(t)')

Zooming in, we can see what LineNormals2D did more explicitly:

h = impoly(ax(2),posn);
setColor(h,[1 0 1])
ax.xlim
set(ax,'xlim',[165 230],'ylim',[1120 1200])

There are undoubtedly other ways of creating this mask, but this is what initially occurred to me. Swag, of course, goes to Dirk-Jan for his fine code, but also to anyone who shows me another clever way to implement this functionality!

]]>PicksGray-scale dilation equation
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/AybkU0DxP2c/
<p>A question came into tech support a month or so ago regarding the documentation for imdilate. The question concerned an apparent discrepancy in the equations for binary and grayscale dilation. Here's the formula given for binary dilation:... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2015/05/06/gray-scale-dilation-equation/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1321Wed, 06 May 2015 20:51:59 +0000

A question came into tech support a month or so ago regarding the documentation for imdilate. The question concerned an apparent discrepancy in the equations for binary and grayscale dilation. Here's the formula given for binary dilation:

$$A \oplus B = \{ z | (\hat{B})_z \cap A \neq 0 \}$$

where $\hat{B}$ is the reflection of the structuring element $B$. The reference page goes on to say, "in other words, it is the set of pixel locations $z$ where the reflected structuring element overlaps with foreground pixels of $A$ when translated to $z$."

The reference page then gives another equation for gray-scale dilation.

where $D_B$ is the domain of the structuring element and $A(x,y)$ is assumed to be $-\infty$ outside the domain of the image.

The user who contacted tech support wondered if there might be an error in the equation for gray-scale dilation because the equation doesn't show a reflection of the structuring element. The case was escalated to me for comment.

The gray-scale dilation equation above is correct. (Well, it's correct for about half the world. The other half uses a slightly different form.)

But there are several different but equivalent mathematical equations that can be used to define dilation. Each of these equations corresponds to a different but equivalent geometric interpretation. The equation above can be interpreted as follows:

To compute the output at $(x,y)$, flip (or reflect) $A$ through the origin and then slide the origin pixel over to $(x,y)$. Form the sums of the $A$ pixels with the structuring element heights underneath. Find the maximum of these sums and record the result as the output at $(x,y)$.

As it turns out, dilation is commutative. That suggests that there is a form of the equation, and a corresponding geometric interpretation, in which the structuring element is reflected instead of the image. We can take the first step in that direction via a substitution of variables. Let $q = x - x'$ and $r = y - y'$. Then:

This second equation has the geometric interpretation of leaving the image in placing, flipping (reflecting) and sliding the structuring element, performing sums of the corresponding image pixels and structuring element heights, and then taking the maximum of the sums.

So which equation should we use? Well, both are correct. Which form to use for a real implementation is completely up to the implementer. And these are not the only two equations and geometric interpretations that are valid.

I do think, though, that there is some merit in modifying the gray-scale dilation equation in our documentation to make it more consistent with the form used for binary dilation.

Dear reader, I am curious: do you use nonflat grayscale dilation or erosion in your work? As far as I can tell, there are not many practical applications for using nonflat structuring elements. If you have a use for it, please leave me a comment below.

]]>UncategorizedMATLAB Used to Map EarthQuakes from Satellite Data
http://feedproxy.google.com/~r/mathworks/loren/~3/XCycbVyqNK0/
<p>A friend just pointed out to me a really cool article: <a rel="nofollow" target="_blank" href="http://www.wired.com/2015/04/turns-satellites-work-great-mapping-earthquakes/">Turns Out Satellites Work Great for Mapping Earthquakes</a>. It's about mapping earthquakes using satellite data. This sounded intriguing because I know earth scientists use MATLAB after the fact to analyze seismic data, but I was less certain what exactly they might do with satellite data and earthquakes. The article, centered on work by Bill Barnett from U. of Iowa, is very interesting. Rather than wait for data to be processed well after events occur, Bill demonstrates what can be done with the data in a much shorter time period, allowing first responders to expedite their decision making regarding where, when, and how to respond to damaging events.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/05/06/matlab-used-to-map-earthquakes-from-satellite-data/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1171Wed, 06 May 2015 14:14:22 +0000

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

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

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

]]>NewsAlt Up Down Keyboard Shortcut in MATLAB Editor
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/P-H-EbTXX7g/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2015/05/06/alt-up-down-keyboard-shortcut-in-matlab-editor/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201412/3184/62009828001_3956925921001_alt-up-down-thumbnail5.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 02:39</span>
</div>
</a></div><p>When editing code in the MATLAB editor, I often find that I’m interested in the use of one variable throughout the code. The keyboard shortcut Alt plus the Up or Down arrow keys lets me search through all instances of a variable, so I can navigate quickly through all related... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2015/05/06/alt-up-down-keyboard-shortcut-in-matlab-editor/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1615Wed, 06 May 2015 11:41:05 +0000When editing code in the MATLAB editor, I often find that I’m interested in the use of one variable throughout the code. The keyboard shortcut Alt plus the Up or Down arrow keys lets me search through all instances of a variable, so I can navigate quickly through all related points in the script or function.

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

]]>Format: VideoSpotted: a Real MATLAB Fan
http://feedproxy.google.com/~r/mathworks/desktop/~3/hGRCbQtFIDk/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/community/files/matlab-tattoo.jpg"/></div><p>Apparently, it’s not temporary. Impressive!
... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/community/2015/05/05/spotted-a-real-matlab-fan/">read more >></a></p>http://blogs.mathworks.com/community/?p=3191Tue, 05 May 2015 21:12:31 +0000Apparently, it’s not temporary. Impressive!
]]>UncategorizedParallel Random Number Generators
http://feedproxy.google.com/~r/mathworks/moler/~3/A0LwkZhsTKo/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/cleve/random_blog2_01.png"/></div><p>This is the second of a multi-part series about the MATLAB random number generators. If you ask for <tt>help rng</tt>, you will get lots of information, including the fact that there are three modern generators.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2015/05/04/parallel-random-number-generators/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1198Mon, 04 May 2015 17:00:54 +0000

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

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

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

Begin by opening a pool of workers.

parpool;

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

Four players each play twenty thousand hands of Blackjack.

p = 4;
n = 20000;

Collect the results in this array.

B = zeros(n,p);

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

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

Plot the results.

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

Parallel streams

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

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

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

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

combRecursive

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

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

A single precision random real u is then produced by

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

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

multFibonacci

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

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

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

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

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

Which one?

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

rng('shuffle')

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

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

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

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

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

]]>Remove Default Excel Sheets
http://feedproxy.google.com/~r/mathworks/pick/~3/nV_mzx8HE8U/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/pick/Sean/mainrmsheets/after.png"/></div><p>
Sean's pick this week is RemoveSheet123 by Noam Greenboim.
Working with MATLAB And Excel
Recently, I've been doing data analysis for colleagues who are not MATLAB users. ... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/05/01/remove-default-excel-sheets/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5961Fri, 01 May 2015 13:00:39 +0000

Recently, I've been doing data analysis for colleagues who are not MATLAB users. However, most of them are familiar with
Excel. Thus, I do the heavy lifting in MATLAB, such as aggregating millions of records from a database, or complex set analysis, and export the results to Excel so they can see it, manipulate it, or pivot it out in Excel.

To facilitate this, my mode of operation has been to take the data and stick it into a table and then use writetable to export to Excel. In general, writetable "does the right thing". If you have heterogenous data, it figures out how to write it, the variable names already exist
so those are used as column headers, etc. I find this to be a simple workflow for most cases where I don't usually need the
additional control xlswrite offers.

There are two things about this approach that have to be taken care of though. First, if you write to a named sheet, which
I always do, you'll get a warning for adding a new sheet to the Excel file.

Fortunately, this can be turned off pretty easily.

warning('off','MATLAB:xlswrite:AddSheet');

I actually turn off this warning in my startup.m file so I never need to see this message except when blogging about it!

The other thing that happens with this approach is that the default Excel sheets are not used, remaining empty and taking
up space or time when clicking on them to see what they contain. And this is where Noam's file comes in: It removes those
default sheets, plain and simple.

% Download Some Data from Yahoo! (Requires Datafeed Toolbox)
y = yahoo;
% Fetch
ClosePrice = fetch(y,'SAM','Close','04/01/15','04/30/15');
% Convert to table and datetime with descriptive variable names
TClosePrice = table(datetime(ClosePrice(:,1),'ConvertFrom','datenum'),ClosePrice(:,2));
TClosePrice.Properties.VariableNames = {'Date','SamAdamsClosingPrice'};
% Write it
writetable(TClosePrice,'StockData.xlsx','Sheet','SamAdams')

Before:

Note: This only works on Windows because it relies on the Excel API.

% Remove sheets
RemoveSheet123('StockData.xlsx');

sheet #1 - deleted.

After:

Thanks Noam for saving me the effort of researching the API for Excel.

I did make one small change which was to comment out the catch block after the delete (line 59) so if it it fails, it does so silently. The reason for this is that I have a different
Excel template that gives me just one sheet by default. I never have to delete all three, at least not on this machine.
This change takes care of the fact that when the error is thrown, the document is not saved or exited. I'd recommend using
an onCleanup object to prevent this state from being possible.

Comments

How do you integrate MATLAB and Excel? Any tips, tricks, or challenges, you'd like to share for some MathWorks swag?

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

]]>Undo Parameter Changes!
http://feedproxy.google.com/~r/SethOnSimulink/~3/HBe_8LFJuAQ/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/seth/2015Q2/undo.gif"/></div><p>This week I want to highlight an improvement to the Undo feature in Simulink.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/04/28/undo-parameter-changes/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4484Wed, 29 Apr 2015 01:12:25 +0000This week I want to highlight an improvement to the Undo feature in Simulink.

Undo Parameter Changes

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

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

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

Now it's your turn

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

]]>Under New Management
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/GpQZGAl9R_M/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2015/04/24/under-new-management/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201504/3049/62009828001_4194782607001_under-new-management-thumbnail10-210x118.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 02:04</span>
</div>
</a></div><p>As I will be taking over Doug’s blog, I thought I’d tell you a bit about myself and what I hope to cover in future posts.
Feel free to make requests for topics in the comments area at any time.
... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2015/04/24/under-new-management/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1612Fri, 24 Apr 2015 19:52:40 +0000As I will be taking over Doug’s blog, I thought I’d tell you a bit about myself and what I hope to cover in future posts.

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

]]>Format: VideoParsing JSON files
http://feedproxy.google.com/~r/mathworks/pick/~3/QlliKehTZgw/
<p>
Jiro's pick this week is JSONlab by Qianqian Fang.
Recently, I've been working quite a bit with JSON files. I needed a way to read the files and parse out some content. With structured files like these, it's not terribly difficult
... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/pick/2015/04/24/parsing-json-files/">read more >></a></p>http://blogs.mathworks.com/pick/?p=5953Fri, 24 Apr 2015 13:00:20 +0000

Recently, I've been working quite a bit with JSON files. I needed a way to read the files and parse out some content. With structured files like these, it's not terribly difficult
to find tags and extract out information. But it does take some time. I started to write a program to do this.

What was I thinking?! Why re-invent the wheel when, most likely, someone has already done this before. It turns out that there
are already a few File Exchange entries for dealing with JSON formats. Almost any of them would have helped me, but Qianqian's JSONlab caught my eyes. His entry
directly handled JSON files, as opposed to JSON strings.

Let's see it in action. Suppose that there is a JSON file, "example.json", whose content is this.

streetAddress: '3 Apple Hill Dr'
city: 'Natick'
state: 'MA'
postalCode: '01760'

disp(data.phoneNumber{1})

type: 'home'
number: '123 456 7890'

disp(data.phoneNumber{2})

type: 'cell'
number: '098 765 4321'

In addition to decoding JSON files into MATLAB structures, the package lets you convert MATLAB structures into JSON format
strings, as well as deal with binary JSON formats.

Thanks Qianqian! You saved me some time with my work!

Comments

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

]]>The Netflix Prize and Production Machine Learning Systems: An Insider Look
http://feedproxy.google.com/~r/mathworks/loren/~3/mA72stV8Su4/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/netflix_02.png"/></div><p>Do you watch movies on Netflix? Binge-watch TV series? Do you use their movie recommendations? Today's guest blogger, Toshi Takeuchi, shares an interesting blog post he saw about how Netflix uses machine learning for movie recommendations.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/04/22/the-netflix-prize-and-production-machine-learning-systems-an-insider-look/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1159Wed, 22 Apr 2015 13:10:27 +0000

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

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

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

How Recommender Systems Work

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

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

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

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

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

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

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

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

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

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

ngh_corr =
0.9661 0.9439 0.8528 0 -0.4402 -0.8315

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

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

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

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

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

RMSE =
0.5567

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

What Netflix did with the winning solutions

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

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

The number of ratings in the competition dataset was 100 million, but the actual production system had over 5 billion

The competition dataset was static, but the number of ratings in the production system keeps growing (4 million ratings per day when the blog post was written)

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

Additional accuracy gains that we measured did not justify the engineering effort needed to bring them into a production environment

The Netflix business model changed from DVD rental to streaming and that changed the way data is collected and recommendations are delivered.

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

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

So, Was It Worth $1M?

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

Lessons Learned: New Metrics

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

"75% of what people watch is from some sort of recommendation"

"continuously optimizing the member experience and have measured significant gains in member satisfaction"

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

Lessons Learned: System Architecture

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

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

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

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

Offline to process data - pre-computes the time-consuming steps in a batch process.

Online to process requests - responds to user action instantaneously by taking advantage of the Offlne and Nearline outputs

Nearline to process events - bridges the two sub-systems by pre-computing frequently performed operations and caching the result in advance of active user action

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

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

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

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

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

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

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

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

As noted earlier, Collaborative Filtering and Market Basket Analysis are closely related applications of data mining and machine learning. Both aim to learn from the customer behaviors, but Market Basket Analysis aims for high frequency transactions while Collaborative Filtering enables personalized recommendations. Naturally, you cannot restock your store shevles for individual shoppers, but you can in ecommerce. You can also think of Latent Semantic Analysis as a related technique applied in text analytics domain.

Just as Market Basket Analysis found its way into web usage data mining, you can probably use Colaborative Filtering for the web data analysis. There may be applications in other fields as well. Any thoughts about your creative use of "Collaborative Filtering" with your data? Let us know here.

]]>Random Number Generators, Mersenne Twister
http://feedproxy.google.com/~r/mathworks/moler/~3/FYIIKDK5X6o/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/cleve/random_blog_01.png"/></div><p>This is the first of a multi-part series about the MATLAB random number generators.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2015/04/17/random-number-generator-mersenne-twister/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1192Fri, 17 Apr 2015 16:26:15 +0000

This is the first of a multi-part series about the MATLAB random number generators.

If you issue the following commands at any point in any recent version of MATLAB, you will always get this plot.

rng default
hist(randn(10000,1),100)

The rng command controls the random number generator that is used by the rand, randn, and randi functions. When called with the default parameter, rng resets the generator to the condition that it has when a fresh MATLAB is started. To see what that condition is, just call rng by itself.

rng

ans =
Type: 'twister'
Seed: 0
State: [625x1 uint32]

We see that the default generator is 'twister', the default seed is 0, and that the state is a length 625 vector of unsigned 32-bit integers.

If you ask for help rng, you will get lots of information, including the fact that there are three modern generators.

The remainder of today's post is about 'twister'. I will cover the others in future posts.

Mersenne Twister

Mersenne Twister is, by far, today's most popular pseudorandom number generator. It is used by every widely distributed mathematical software package. It has been available as an option in MATLAB since it was invented and has been the default for almost a decade.

Mersenne Twister was developed by professors Makoto Matsumoto and Takuji Nishimura of Hiroshima University almost twenty years ago. Here is their home page. The C source code is available here.

Mersenne primes

Mersenne primes are primes of the form 2^p - 1 where p itself is prime. They are named after a French friar who studied them in the early 17th century. We learn from Wikipedia that the largest known prime number is the Mersenne prime with p equal to 57,885,161. The Mersenne Twister has p equal to 19937. This is tiny as far as Mersenne primes go, but huge as far as random number generators are concerned.

Why the name?

Matsumoto explains how the name "Mersenne Twister" originated in the Mersenne Twister Home Page.

MT was firstly named "Primitive Twisted Generalized Feedback Shift
Register Sequence" by a historical reason.
Makoto: Prof. Knuth said in his letter "the name is mouthful."
Takuji: ........
a few days later
Makoto: Hi, Takkun, How about "Mersenne Twister?" Since it uses Mersenne
primes, and it shows that it has its ancestor Twisted GFSR.
Takuji: Well.
Makoto: It sounds like a jet coaster, so it sounds quite fast, easy to
remember and easy to pronounce. Moreover, although it is a secret, it
hides in its name the initials of the inventors.
Takuji: .......
Makoto: Come on, let's go with MT!
Takuji: ....well, affirmative.
Later, we got a letter from Prof. Knuth saying "it sounds a nice name." :-)

Algorithm

The integer portion of the Mersenne twister algorithm does not involve any arithmetic in the sense of addition, subtraction, multiplication or division. All the operations are shifts, and's, or's, and xor's.

All the elements of the state, except the last, are unsigned 32 bit random integers that form a cache which is carefully generated upon startup. This generation is triggered by a seed, a single integer that initiates the whole process.

The last element of the state is a pointer into the cache. Each request for a random integer causes an element to be withdrawn from the cache and the pointer incremented. The element is "tempered" with additional logical operations to improve the randomness. When the pointer reaches the end of the cache, the cache is refilled with another 623 elements.

The algorithm is analyzed by investigating the group theoretic properties of the permutation and tempering operations. The parameters have been chosen so that the period is the Mersenne prime 2^19937-1. This period is much longer than any other random number generator proposed before or since and is one of the reasons for MT's popularity.

By design, the results generated satisfy an equidistribution property in a 623-dimensional cube.

Doubles

Here is the function in the Mersenne Twister source code that converts a pair of random uint32s into a random double. You can see that it takes the top 27 bits of one int and the top 26 bits of the other, sticks them together, and multiplies by what MATLAB would call eps/2. This is the only place in the code where floating point arithmetic is involved.

double genrand_res53(void)
{
unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
return(a*67108864.0+b)*(1.0/9007199254740992.0);
}

The result would be zero in the highly unlikely event that both a and b are zero. If that happens, the MATLAB interface rejects the result and calls this function again. So the smallest double results when a equals zero and b equals 1. The largest double results when both a and b are all 1's. Consequently, the output from rand is in the closed interval

$$ 2^{-53} \leq x \leq 1-2^{-53} $$

Seeds, Streams and State

For more about random number seeds, streams, and state, see Peter Perkins, guest blogger in Loren's blog. And, of course, see the documentation.

Thanks

Thanks to Peter Perkins for the work he has done on our random number suite over the years, and for enlightening me.

References

M. Matsumoto and T. Nishimura. "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudorandom Number Generator." ACM Transactions on Modeling and Computer Simulation, 8(1):3-30. 1998. Available online at: <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf>.

]]>Creating Test Harnesses with Simulink Test!
http://feedproxy.google.com/~r/SethOnSimulink/~3/mdeMYoP2pxU/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/seth/2015Q1/harnessInModel.png"/></div><p>In R2015a, we introduced a new product called Simulink Test. This product offers many great features like a Test Sequence block, various ways to test results of a model against validated data, and a test manager interface.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/04/11/creating-test-harnesses-with-simulink-test/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4503Sat, 11 Apr 2015 12:56:50 +0000In R2015a, we introduced a new product called Simulink Test. This product offers many great features like a Test Sequence block, various ways to test results of a model against validated data, and a test manager interface.

Among all the feature of Simulink Test, the one that I am the most interested in is the Test Harness. I think this will make developing and debugging models more efficient. Let's see how this works.

What is a Test Harness?

Let's say I am developing or debugging a model with multiple components. To illustrate the harness concept, we will use an example model I like: sf_electrohydraulics. This example is great because it contains multiple components of various domains: electrical, hydraulics, mechanical, etc.

Before R2015a, when a model like this one was giving me unexpected results, what I ended up doing is pasting the component I suspected to be problematic in a new empty model and feed it known inputs to see if I get the output I expect.

With the Simulink Test Test Harness feature, this test or debug model can now be part of the original model, making it very easy to manage (no additional files), and to switch between the large model and the test harness.

Creating a Harness

To get started, right-click on a subsystem and select Create Test Harness:

A dialog will open where you can set the properties of the harness. When the harness is created, it can be created with standard blocks like Inport, Signal Builder, From Workspace, etc.

Then you select the harness objective.

What is the meaning of those objective?

Prototyping: If your original model does not compile, that's what you need. The harness will be created without knowledge of input/output signals properties like dimensions and data type. This should allow you to debug and figure out why the original model does not update.

Refinement/Debugging: In this case, the original model will be compiled and blocks will be inserted in the harness to ensure the input and output signals have the same properties, dimensions, data types, etc.

Verification: In the harness, the subsystem will be marked as read-only so it can be validated, but not modified. In addition, SIL and PIL mode verification options are enabled for the component under test. The harness is also rebuilt each time it is opened so that the compiled attributes enforced in the harness are re-computed from the main model and brought up-to-date.

Custom: Four independent checkboxes will be enabled for you to combine the properties of the three above ojectives the way you want.

When you click OK the test harness opens.

Navigating Between the Harnesses and the Model

When a subsystem has harnesses, you will notice a new icon on its bottom right corner. Clicking on it will list the harnesses.

To avoid confusions, you can open only one harness at a time and the original model cannot be edited when a harness is open. In the harness, you can modify the subsystem as you wish, and the changes will be propagated to the original model as soon as you close the harness.

Now it's your turn

Do you think the Test Harness feature will be useful for you? Let us know how you plan to use it by leaving a comment here.

]]>Can You Find Love through Text Analytics?
http://feedproxy.google.com/~r/mathworks/loren/~3/A81epMQ9xVU/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/findingLoveUpdate2_02.png"/></div><p><a rel="nofollow" target="_blank" href="https://www.youtube.com/watch?v=qtsNbxgPngA">Jimmy Fallon Blew a Chance to Date Nicole Kidman</a>, but do you know there is supposedly a way to fall in love with anyone? Today's guest blogger, Toshi Takeuchi, would like to talk about finding love with MATLAB.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/04/08/can-you-find-love-through-text-analytics/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1134Wed, 08 Apr 2015 14:03:05 +0000

Jimmy Fallon Blew a Chance to Date Nicole Kidman, but do you know there is supposedly a way to fall in love with anyone? Today's guest blogger, Toshi Takeuchi, would like to talk about finding love with MATLAB.

"Two heterosexual strangers sat face to face in a lab and answered a series of 36 increasingly personal questions. Then they stared silently into each other's eyes for four minutes. Six months later, they were married."

I wanted to see if someone could try it. Luckily, a friend of mine in Japan was keen to give it a try, but there was one minor issue: she couldn't find any male counterpart who was willing to join her in this experiment.

This is a big issue in Japan where the birthrate went negative. There is even a new word, Konkatsu, for the intensive effort required to get married. Before we can do this experiment, we need to solve this problem first. A lot of people turn to online dating for that, but that is not so easy, either. Do you need some evidence?

In an online dating world you need to comb through a mind-numbing volume of profiles just to get started. Then came the idea: why not use MATLAB to mine online profiles to find your love?

We need data to analyze. I don't have access to real online dating profiles, but luckily I found Online Dating Ipsum by Lauren Hallden that randomly generates fictitious ones. I used Latent Semantic Analysis (LSA) to cluster online profiles based on the words they contain. I cooked up a MATLAB class myLSA.m to implement Latent Semantic Analysis methods. Let's initialize it into an object called LSA, and load the dataset and print one of those.

Working at a coffee shop adventures tacos medical school. Feminism going
to the gym strong and confident Family Guy listening to music, my beard
Kurosawa discussing politics trying different restaurants I know I listed
more than 6 things. Snowboarding no drama outdoor activities discussing
politics pickles my friends tell me they don't get why I'm single.

Not bad for a random word salad, except that they are all male profiles. If you need female profiles, you need to find other sources.

Text Processing Pipeline

Before we can analyze text, we need to process it into an appropriate form. There is a fairly standard process for English text.

Tokenization: split text into word tokens using white space, etc.

Standardization: standardize word forms, i.e., all lowercase

Stopwords: remove common words, such as 'the, a, at, to'

Stemming: reduce words into their root forms by trimming their endings

Indexing: sort the words by document and count word frequencies

Document-Term Frequency Matrix: turn indexed frequency counts into a document x term matrix

The tokenizer method takes care of the first four steps - tokenization, normalization, stopwords and stemming. Check out the before and after.

tokenized = LSA.tokenizer(profiles.Profile);
before = profiles.Profile(1)
after = {strjoin(tokenized{1},' ')}

before =
'Working at a coffee shop adventures tacos medical school. Feminism goi...'
after =
'work coffe shop adventur taco medic school femin go gym strong confid ...'

Next, the indexer method creates word lists and word count vectors.

Then we create a document-term frequency matrix from these using docterm. The minimum frequency is set to 2 and that drops any words that only occur once through the entire collection of documents.

docterm = LSA.docterm(word_lists,word_counts,2);

TF-IDF Weighting

You could use the document-term frequency matrix directly, but raw word count is problematic - it gives too much weight to frequent words, and frequent words that appear in many documents are usually not so useful to understand the differences among those documents. We would like to see the weight to represent the relevancy of each word.

TF-IDF is a common method for frequency weighting. It is made up of TF, which stands for Term Frequency, and IDF, Inverse Document Frequency. TF scales based on the number of times a given term appears in a document, and IDF inversely scales based on how many document a given term appears in. The more frequently a word appears in documents, the less weight it gets. TF-IDF is just a product of those two metrics. Let's use tfidf to apply this weighting scheme. It also optionally returns TF.

tfidf = LSA.tfidf(docterm);

I went through each step of text processing, but we could instead run vectorize to turn a raw cell array of online dating profiles into a TF-IDF weighted matrix in one shot.

tfidf = LSA.vectorize(profiles.Profile,2);

Low-Rank Approximation

Once the data is transformed into a matrix, we can apply linear algebra techniques for further analysis. In LSA, you typically apply singular value decomposition (SVD) to find a low-rank approximation.

Let's first get the components of SVD. U is the SVD document matrix, V is the SVD term matrix, and S is the singular values.

[U,S,V] = svd(tfidf);

If you square S and divide it by sum of S squared, you get the percentage of variance explained. Let's plot the cumulative values.

explained = cumsum(S.^2/sum(S.^2));
figure
plot(1:size(S,1),explained)
xlim([1 30]);ylim([0 1]);
line([5 5],[0 explained(5)],'Color','r')
line([0 5],[explained(5) explained(5)],'Color','r')
title('Cumulative sum of S^2 divided by sum of S^2')
xlabel('Column')
ylabel('% variance explained')

You see that the first 5 columns explain 60% of variance. A rank-5 approximation will retain 60% of the information of the original matrix. The myLSA class also provides lowrank that performs SVD and returns a low rank approximation based on some criteria, such as number of columns or the percentage of variance explained.

[Uk,Sk,Vk] = LSA.lowrank(tfidf,0.6);

Visualize Online Dating Profiles

We can also use the first 2 columns to plot the SVD document matrix U and SVD term matrix V in 2D space. The blue dots represent online dating profiles and words around them are semantically associated to those profiles.

figure()
scatter(U(:,1), U(:,2),'filled')
title('Online Dating Profiles and Words')
xlabel('Dimension 1')
ylabel('Dimension 2')
xlim([-.3 -.03]); ylim([-.2 .45])
for i = [1,4,9,12,15,16,20,22,23,24,25,27,29,33,34,35,38,47,48,53,57,58,...
64,73,75,77,80,82,83,85,88,97,98,103,113,114,116,118,120,125,131,...
136,142,143,156,161,162,166,174,181,185,187,199,200,204,206,212,...
234,251]
text(V(i,1).*3, V(i,2).*3, LSA.vocab(i))
end
text(-0.25,0.4,'Wholesome/Sporty','FontSize', 12, 'Color', 'b')
text(-0.15,-0.15,'Bad Boy/Colorful','FontSize', 12, 'Color', 'b')

You can see there are two main clusters - what I would call the "Wholesome/Sporty" cluster and one called the "Bad Boy/Colorful" cluster, based on the words associated with them. This makes sense, because Lauren provides two options in her profile generator:

Typical inane jabber

With a side of crazy sauce

Can you guess which cluster belongs to which category?

Now you can cluster a whole bunch of profiles at once and quickly eliminate those that don't match your taste. You can also add your own profile to see which cluster you belong to, and, if that puts you in a wrong cluster of profiles, then you may want to update your profile.

Computing Similarity

Say you find a cluster of profiles you are interested in. Among the profiles you see there, which one is the closest to your taste? To answer this question, we need to find a way to define the similarity of two documents. If you use the Euclidean distance between vectors, the longer documents and shorter documents can have very different values even if they share many of the same words. Instead, we can use the angle between the vectors to determine the similarity. This is known as Vector Space Model. For ease of computation, cosine is used for similarity computation.

For practical implementation, you can just length normalize vectors by the L2 norm, and compute the dot product.

cosine = dot(A/norm(A),B/norm(B))

You can apply length normalization ahead of the similarity computation. We will use the rank-5 approximation of the SVD document matrix to compare online dating profiles using normalize.

doc_norm = LSA.normalize(U(:,1:5));

Now we can compute cosine similarities between profiles with score. Let's compare the first profile to the first five profiles.

LSA.score(doc_norm(1:5,:),doc_norm(1,:))

ans =
1
0.20974
0.55248
0.97436
0.72994

The first score is 1, which means it is a perfect match, and that's because we are comparing the first profile to itself. Other profiles got lower scores depending on how similar they are to the first profile.

Getting the Ranked Matches

It's probably useful if you can describe your ideal date and find the profiles that match your description ordered by similarity. It is a bit like a search engine.

To compare the new text string to the pre-computed matrix, we need to apply the same pre-processing steps that we have already seen. query can take care of the tedious details.

q = 'someone fun to hang out with, good sense of humor, likes sushi,';
q = [q 'watches Game of Thrones, sees foreign films, listens to music,'];
q = [q 'do outdoor activities or fitness'];
weighted_q = LSA.query(q);

Now we need to transform the query vector into the rank-5 document space. This is done by transforming M = U*S*V' into U = M'*V*S^-1 and substituting M' with the query vector and V and S with their low rank approximations.

The myLSA class also provides the reduce method to perform the same operation.

q_reduced = LSA.reduce(weighted_q);

Then we can length-normalize the query vector and compute the dot products with documents. Let's sort the cosine similarities in descending order, and check the top 3 results.

q_norm = LSA.normalize(q_reduced);
[scores,idx] = sort(LSA.score(doc_norm,q_norm),'descend');
disp('Top 3 Profiles')
for i = 1:3
profiles.Profile(idx(i))
end

Top 3 Profiles
ans =
'Someone who shares my sense of humor fitness my goofy smile Oxford com...'
ans =
'My cats I'm really good at my goofy smile mountain biking. Fixing up m...'
ans =
'My eyes just looking to have some fun if you think we have something i...'

Looks pretty reasonable to me!

In this example, we applied TF-IDF weighting to both the document-term frequency matrix as well as the query vector. However, you only need to apply IDF just once to the query to save computing resources. This approach is known as lnc.ltc in SMART notation system. We already processed our query in ltc format. Here is how you do lnc for your documents - you use just TF instead of TF-IDF:

Can my Japanese friends benefit from this technique? Yes, definitely. Once you have the document-term frequency matrix, the rest is exactly the same. The hardest part is tokenization, because there is no whitespace between words in Japanese text.

Fortunately, there are free tools to do just that - they are called Japanese Morphological Analyzers. One of the most popular analyzers is MeCab (this link goes to a Japanese page). A binary package is available for installation on Windows, but it is for 32-bit, and doesn't work with MATLAB in 64-bit. My Japanese colleague, Takuya Otani, compiled the source code to run it on MATLAB 64-bit on Windows.

MATLAB provides an interface to shared libraries like DLLs, and we can use loadlibrary to load them into memory and access functions from those shared libraries. Here is an example of how to call the shared library libmecab.dll that Takuya compiled.

You may not have any particular need for handling Japanese text, but it gives you a good example of how to load a DLL into MATLAB and call its functions. Please note some requirements in case you want to try:

Have a 64-bit Japanese Windows computer with 64-bit MATLAB

Follow Takuya's instructions to compile your own 64-bit DLL and place it in your current folder along with its header file.

loadlibrary('libmecab.dll', 'mecab.h');

When you run this command, you may get several warnings, but you can ignore them. If you would like to see if the library was loaded, use libfunctionsview function to view the functions available in the DLL.

libfunctionsview('libmecab')

To call a function in the DLL, use calllib. In the case of Mecab, you need to initialize Mecab and obtain its pointer first.

As an example, let's call one of the Mecab functions you can use to analyze Japanese text - mecab_sparse_tostr.

text = 'Some Japanese text';
result = calllib('libmecab', 'mecab_sparse_tostr', mecab, text);

When finished, clear the pointer and unload the DLL from the memory using unloadlibrary.

clearvars mecab
unloadlibrary('libmecab')

Call for Action

If you happen to be single and are willing to try the experiment described in the New York Times article, please report back here with your results. The New York Times now provides a free app to generate 36 magical questions!

]]>So Long, and Thanks for all the Fantastic Videos
http://feedproxy.google.com/~r/DougsMatlabVideoTutorials/~3/gyeGVatvOo8/
<div class="thumbnail thumbnail_asset asset_overlay video"><a rel="nofollow" target="_blank" href="http://blogs.mathworks.com/videos/2015/04/08/so-long-and-thanks-for-all-the-fantastic-videos/?dir=autoplay"><img src="https://bcsecure01-a.akamaihd.net/6/62009828001/201504/937/62009828001_4160682258001_StuartMcGarrity-62009828001-4147528847001-farewell-to-doug-thumbnail1b.jpg?pubId=62009828001"/>
<div class="overlay_container">
<span class="icon-video icon_color_null"> 01:12</span>
</div>
</a></div><p>Doug has left the MathWorks but watch this video to hear from who is taking over his blog.
... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/videos/2015/04/08/so-long-and-thanks-for-all-the-fantastic-videos/">read more >></a></p>http://blogs.mathworks.com/videos/?p=1595Wed, 08 Apr 2015 13:49:45 +0000Doug has left the MathWorks but watch this video to hear from who is taking over his blog.

]]>Format: VideoDisplaying a color gamut surface
http://feedproxy.google.com/~r/SteveOnImageProcessing/~3/kWDPjdPhvdQ/
<div class="overview-image"><img src="http://blogs.mathworks.com/steve/files/srgb_gamut_surface_08.png" class="img-responsive attachment-post-thumbnail wp-post-image" alt="srgb_gamut_surface_08"/></div><p>Today I'll show you one way to visualize the sRGB color gamut in L*a*b* space with assistance with a couple of new functions introduced last fall in the R2014b release. (I originally planned to post this a few months ago, but I got sidetracked writing about colormaps.)... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/steve/2015/04/03/displaying-a-color-gamut-surface/">read more >></a></p>http://blogs.mathworks.com/steve/?p=1311Fri, 03 Apr 2015 20:40:04 +0000

Today I'll show you one way to visualize the sRGB color gamut in L*a*b* space with assistance with a couple of new functions introduced last fall in the R2014b release. (I originally planned to post this a few months ago, but I got sidetracked writing about colormaps.)

The first new function is called boundary, and it is in MATLAB. Given a set of 2-D or 3-D points, boundary computes, well, the boundary.

Here's an example to illustrate.

x = gallery('uniformdata',30,1,1);
y = gallery('uniformdata',30,1,10);
plot(x,y,'.')
axis([-0.2 1.2 -0.2 1.2])
axis equal

Now compute and plot the boundary around the points.

k = boundary(x,y);
hold on
plot(x(k),y(k))
hold off

"But, Steve," some of you are saying, "that's not the only possible boundary around these points, right?"

Right. The function boundary has an optional shrink factor that you can specify. A shrink factor of 0 corresponds to the convex hull. A shrink factor of 1 gives a compact boundary that envelops all the points.

k0 = boundary(x,y,0);
k1 = boundary(x,y,1);
hold on
plot(x(k0),y(k0))
plot(x(k1),y(k1))
hold off
legend('Original points','Shrink factor: 0.5 (default)',...'Shrink factor: 0','Shrink factor: 1')

Here's a 3-D example using boundary. First, the points:

P = gallery('uniformdata',30,3,5);
plot3(P(:,1),P(:,2),P(:,3),'.','MarkerSize',10)
grid on

k = boundary(P);
hold on
trisurf(k,P(:,1),P(:,2),P(:,3),'Facecolor','red','FaceAlpha',0.1)
hold off

The second new function I wanted to mention is rgb2lab. This function is in the Image Processing Toolbox. The toolbox could convert from sRGB to L*a*b* before, but this function makes it a bit easier. (And, if you're interested, it supports not only sRGB but also Adobe RGB 1998).

Just for grins, let's reverse the a* and b* color coordinates for an image.

Now let's get to work on visualizing the sRGB gamut surface. The basic strategy is to make a grid of points in RGB space, transform them to L*a*b* space, and find the boundary. (We'll use the default shrink factor.)

[r,g,b] = meshgrid(linspace(0,1,50));
rgb = [r(:), g(:), b(:)];
lab = rgb2lab(rgb);
a = lab(:,2);
b = lab(:,3);
L = lab(:,1);
k = boundary(a,b,L);
trisurf(k,a,b,L,'FaceColor','interp',...'FaceVertexCData',rgb,'EdgeColor','none')
xlabel('a*')
ylabel('b*')
zlabel('L*')
axis([-110 110 -110 110 0 100])
view(-10,35)
axis equal
title('sRGB gamut surface in L*a*b* space')

Here are another couple of view angles.

view(75,20)

view(185,15)

That's it for this week. Have fun with color-space surfaces!

]]>UncategorizedTeacher’s Pet Students’ Robotics Challenge – Enter today for a chance to win 3-D printer
http://feedproxy.google.com/~r/SethOnSimulink/~3/LD50d7lL_YU/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/seth/2015Q2/piKit.png"/></div><p>Today my colleague Madhu Govindarajan is here to talk about the Teacher’s Pet Students’ Robotics Challenge... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/seth/2015/04/03/teachers-pet-students-robotics-challenge-enter-today-for-a-chance-to-win-3-d-printer/">read more >></a></p>http://blogs.mathworks.com/seth/?p=4547Fri, 03 Apr 2015 10:15:28 +0000Today my colleague Madhu Govindarajan is here to talk about the Teacher’s Pet Students’ Robotics Challenge

If your proposal is chosen, you will have the choice between two kits.

The Arduino kit, with a wheeled platform:

The Raspberry Pi kit, including a Pi Camera.

Now it's your turn

What kind of robot are you dreaming of? Something to carry your books? To take care of your cat while you are studying? A robot that watches you to ensure you do not fall asleep while studying? To cook for you while you are studying? The options are endless!

Get to the drawing board, and let us know if you have questions.

]]>The Winter of Our Vectorization
http://feedproxy.google.com/~r/mathworks/loren/~3/ZFqHCthG1TA/
<div class="overview-image"><img class="img-responsive" src="http://blogs.mathworks.com/images/loren/2015/winterOfOurVectorization_01.png"/></div><p>Today I'd like to introduce guest blogger Matt Tearle who works on our MATLAB Product Training materials here at MathWorks. Matt is on a mission to teach the world MATLAB, but this winter is testing his resolve. Annoyed that 22" of snow forced him to reschedule a training, today he... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/loren/2015/04/01/the-winter-of-our-vectorization/">read more >></a></p>http://blogs.mathworks.com/loren/?p=1162Wed, 01 Apr 2015 13:19:12 +0000

Today I'd like to introduce guest blogger Matt Tearle who works on our MATLAB Product Training materials here at MathWorks. Matt is on a mission to teach the world MATLAB, but this winter is testing his resolve. Annoyed that 22" of snow forced him to reschedule a training, today he shows just how bad this winter had been.

Going the Whole 9 Feet

It has been a brutal winter at MathWorks' headquarters in Natick, MA (a little west of Boston). I recently read a fascinating article by MIT doctoral student Ben Letham who created a convincing visualization of just how relentless the snowfall has been this year. It was so convincing that I had to reproduce it myself (in MATLAB, of course!). In the process I had some fun with vectorization.

It's common to look at statistics like the largest amount of snowfall in a single day (23.5" on Feb 17, 2003) or the highest total for the entire winter (108.5" for 2014/15). But 6" of snow every day for a week would also cause chaos, even though it would seem mild by either of those measures. Hence, Mr. Letham's idea is to calculate the largest total of snow that fell during a range of given time periods -- three days, a week, two weeks, ... For example, 40.5" of snow fell during the worst week-long period this winter, compared to just over 31" in the worst week of 1995/96 (which had the most snow of an entire winter until this year) and only(!) 28.5" in 2002/03 (23.5" of which fell in one day).

When you plot these "worsts" for a range of time periods, you get:

Sure enough, 2014/15 had the worst week ever, the worst two weeks ever, the worst month ever, the worst of just about any time period ever! And by some distance, too. The winter of 1977/78 is included because Boston locals talk with nostalgic reverence about the blizzard of 1978 that brought the area to a standstill. 1978 was pretty bad, but nothing compared to the snowfall we've had this winter.

Enter MATLAB

Intrigued by Mr. Letham's analysis, I wanted to play with the numbers myself. Like any good researcher, Mr. Letham documented his sources, which turned out to be NOAA's National Climatic Data Center. This meant that I could retrieve the data for myself and start playing (and so can you!). Before doing my own thing (which turned out to be fun, but a story for another day), I wanted to check that I could reproduce Mr. Letham's results. As it turned out, this helped me remember why I love MATLAB and gain appreciation for a function I had previously overlooked.

I downloaded the data from NOAA, used the Import Tool to import it into MATLAB, and did some basic preprocessing. As I often do, at this point I saved the cleaned-up data as a MAT-file, so I could now play with impugnity.

load snowfalls

Time for Some Logical Indexing

To compare different winters, I need to be able to extract the snowfall data for the given time periods. This is easily achieved with a combination of the date and time variable types that were introduced in R2014b and my all-time favorite MATLAB construction: logical indexing. I used the Import Tool to read the dates as datetime variables. Now I needed to separate 2014/15 from the rest of the data. Because the Northern Hemisphere winter crosses the new year, I decided to break years by the summer solstice (June 21):

Numerical comparisons work naturally with datetime variables, so pre2015 is a logical variable which I can use to extract the snowfall values. To extract the year of the legendary blizzard, I used the handy isbetween function:

So now the crux of the problem: how to calculate the largest total over an $n$-day period. The typical programming language approach would be to loop over the array, taking $n$ elements starting with the $k^{th}$ element. But thinking about this process made me realize that I was looking for something very similar to a moving average. I knew that you can calculate moving averages in MATLAB by using convolution, implemented with the conv function. In this case, I want the sum (rather than the sum divided by $n$, as I'd want for an average). Hence, I can calculate all the $n$-day totals in a vectorized way:

n = 7;
ndaytotals = conv(snow2015,ones(n,1));
max7day = max(ndaytotals)

max7day =
40.4724

Fun with arrayfun

Now I need to loop over $n$. I couldn't see any neat way to vectorize this, so maybe I should use a loop (remembering to preallocate):

nmax = 180;
ndaysmax = zeros(nmax,1);
for n = 1:nmax
ndaysmax(n) = max(conv(snow2015,ones(n,1)));
end

I could also use arrayfun to perform this same operation without explicitly writing a for-loop:

But is this really achieving anything? The code is very slightly more compact, but arguably more complicated (at least to someone who isn't a hard-core MATLAB-phile). I admit that, for this very reason, I tend to look down a bit on arrayfun as "how to cheat at vectorizing". But in this case I realized I had actually stumbled upon a great reason to use arrayfun...

Having applied this clever convolution operation to the 2015 data, I now need to do the same to the pre-2015 data and the 1978 data (and, presumably, any others that might strike me as interesting). Should I copy and paste code? No! This sounds like a job for a function. The "classical" approach would be to write a function that took an array of data and a maximum value of $n$ as inputs. It would return the maximum total for each $n$ from 1 to $n_{max}$. This wouldn't be a particularly difficult function to write, given that the code is basically done already (above). But do I really need to write a function file just for this? No: anonymous functions provide a way to create functions without writing a separate function file.

The one important caveat, though, is that anonymous functions must have a one-line definition (i.e., no loops!). So here's a compelling reason to use arrayfun: I can define my operation in a single line, enabling me to create an anonymous function, which can then be applied to any data set of interest:

Now I have a fairly short script that performs the whole analysis, leveraging the power of an anonymous function vectorized with arrayfun:

load snowfalls% Extract data for 2014/15 winter and for everything else until then
pre2015 = (t <= datetime(2014,6,21));
snowPre2015 = snow(pre2015);
snow2015 = snow(~pre2015);
% Extract data for the infamous winter of 1977/78
blizz = isbetween(t,datetime(1977,6,22),datetime(1978,6,21));
snow1978 = snow(blizz);
% Define a function that calculates the maximum total value in n% consecutive elements of the array x, for n from 1 to nmax
nmax = 180;
maxXinNdays = @(x) arrayfun(@(n) max(conv(x,ones(n,1))),1:nmax);
% Use the function to compare snowfall totals
figure
semilogx(maxXinNdays(snowPre2015))
hold on
semilogx(maxXinNdays(snow2015))
semilogx(maxXinNdays(snow1978))
xlim([1 200])
% Annotate the graph
xlabel('Number of days')
ylabel('Total snowfall (in)')
title('Max total snowfall in an n-day period')
legend('pre-2015 record','2015','1978','Location','southeast')

Spring into Action!

You can download the data and code and try this out for yourself. (The data also includes temperature readings.) What other statistics or metrics can you calculate with this technique of vectorizing anonymous functions with arrayfun? Was 2015 really that bad? Let us know what you discover about the winter of our discontent here.

]]>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 peek 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 peek 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?

]]>FunNew 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.

]]>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 (Kelly Luetkemeyer) 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... <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 (Kelly Luetkemeyer) 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

]]>Welcome 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.

]]>Triple Precision Accumlated Inner Product
http://feedproxy.google.com/~r/mathworks/moler/~3/yfzyp7_446M/
<p>Single and double precision are combined to facilitate a triple precision accumulated inner product.... <a rel="nofollow" class="read-more" target="_blank" href="http://blogs.mathworks.com/cleve/2015/03/02/triple-precision-accumlated-inner-product/">read more >></a></p>http://blogs.mathworks.com/cleve/?p=1172Mon, 02 Mar 2015 17:00:18 +0000

Single and double precision are combined to facilitate a triple precision accumulated inner product.

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.

]]>UncategorizedLet’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

]]>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>

]]>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 transformsIntroducing 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.

]]>ColormapRobot 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.

]]>UncategorizedRobot 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

]]>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.

]]>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.

]]>UncategorizedIntroducing 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!

]]>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.

]]>UncategorizedAcquire 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: