# Loren on the Art of MATLAB

## Annual Meeting of the American Geophysical Union 2013

MathWorks will have a booth at the annual AGU (American Geophysical Union) meeting this year, where you can expect to see me, Sarah Wait Zaranek, and Shyam Yadati. We will all be happy to show you some of the newest advancements in MATLAB and other MathWorks products. In addition, I will be speaking on Thursday from 1-1:30 PM at the Exhibit Hall Product Theater. My talk is entitled What's new in MATLAB for Earth & Ocean Sciences.

### Contents

#### We're Back!

After several years of hiatus, MathWorks is back at AGU! Please join us either at my talk or at the MathWorks booth were we can update you on some of the newest capabilities in MATLAB relevant to Earth & Ocean sciences.

#### Concepts in the Talk

You can learn how to accelerate:

• Importing and accessing data from excel files, images and data acquisition hardware
• Modeling, analysis and visualization with built-in analysis algorithms, mapping capabilities
• Automating and publishing your work by using MATLAB's built-in publishing capabilities, MATLAB Apps, Community website and Code Generation
• Technical computing by taking advantage of additional cores and GPUs to run algorithms and analyses in parallel

#### Booth Demo Topics

Here's a list of topics we are preparing in advance so we can show them in the booth. Needless to say, our conversations need not be constrained by these topics!

• Working with Complex and Messy Data

Import and handle “real-world” data with MATLAB

• Parallel Computing

Accelerate your computations using multicore desktops, GPU and clusters with MATLAB

• Streaming Data

Process streaming signals and large data with System objects

• New Features in Image Processing and Mapping

Find circles, automatically register images, and create interactive web maps

• Mathematical Modeling

Create models from first principles, using machine learning, or by fitting equations

Build robust maintainable functions, perform unit testing, and use object-oriented programming with MATLAB

• Sharing Results and Encouraging Reproducibility

Automatically create reports, package applications, and generate C code from MATLAB code

#### Will We See You There?

If you plan to come by, especially to talk about any special topics, let me know here.

Get the MATLAB code

Published with MATLAB® R2013b

## In-memory Big Data Analysis with PCT and MDCS

Ken Atwell in the MATLAB product management group is guest blogging about using distributed arrays to perform data analytics on in-memory "big data".

Big data applications can take many forms. Unlike unstructured text processing, common in many business and consumer applications, big data in science and engineering fields can often take a more "regular" form, either a lengthy data set of observations or time series data, or a monolithic, multi-dimensional matrix capturing some problem space.

64-bit MATLAB is fundamentally well-suited to processing data of these types, as long as the data fits within the memory of your computer. When data size outgrows the capability of a single computer, MATLAB Distributed Computing Server (MDCS) can allow MATLAB applications to simply scale and leverage the aggregate memory and computational capabilities of a cluster of computers.

This blog post explores the use of MDCS to analyze a tabular data set that does not fit into the memory of the typical desktop computer. Highlights include:

• Importing a large data set sourced from multiple files with parallel file I/O, and post-processing that data set
• Performing parallel computation across the aggregate data set
• Visualizing the statistical characteristics of the aggregate data set
• Optimizing file I/O performance with MATLAB MAT-Files
• Rebalancing data and workload between workers

You will need Parallel Computing Toolbox (PCT) to access the distributed array, a data type for working with data storage across a cluster. Its usage is virtually identical to that of a normal MATLAB matrix, supporting easy and rapid adoption. The in-memory nature of the distributed array facilitates experimentation and the rapid iteration workflows that MATLAB users have come to expect.

Finally, note that this example has been designed to run in about 20 GB of memory. That is hardly "big data", but the problem size has been purposefully reduced to be small enough to fit into a mid- to high-end workstation, should you not have access to cluster to follow along.

### Contents

This example uses a publicly available data set of over 100 million airline flight records from over a 20-year time span from 1987 to 2008. Download years individually:

As is typical in many big data applications, data is spread over a collection of medium-to-large files, rather than a single very large file. This is in fact advantageous, allowing the data to be more easily distributed across a cluster running MDCS. During this development phase, use a modest pool (8 workers) and only the first eight files in the data set -- the full data set will be seen later in this post.

#### Identify Files to Process

close all
clear all
dataFolder = '/Volumes/HD2/airline/data';
allFiles = dir(fullfile(dataFolder, '*.csv'));
% Remove "dot" files (Mac)
allFiles = allFiles(~strncmp({allFiles(:).name}, '.', 1));
files = allFiles(1:8);  % For initial development


#### Examine the Data

The data set is comma-separated text. To keep things simple in this particular example, only read in a few of the columns in the data set. Importing from text is an expensive operation -- reducing the importation to only the needed columns can speed up a MATLAB application considerably. First look at left-most part of the raw data:

if ~ispc
system(['head -20 ' fullfile(dataFolder, files(8).name) ' | cut -c -60']);
end

Year,Month,DayofMonth,DayOfWeek,DepTime,CRSDepTime,ArrTime,C
1994,1,7,5,858,900,954,1003,US,227,NA,56,63,NA,-9,-2,CLT,ORF
1994,1,8,6,859,900,952,1003,US,227,NA,53,63,NA,-11,-1,CLT,OR
1994,1,10,1,935,900,1023,1003,US,227,NA,48,63,NA,20,35,CLT,O
1994,1,11,2,903,900,1131,1003,US,227,NA,148,63,NA,88,3,CLT,O
1994,1,12,3,933,900,1024,1003,US,227,NA,51,63,NA,21,33,CLT,O
1994,1,13,4,NA,900,NA,1003,US,227,NA,NA,63,NA,NA,NA,CLT,ORF,
1994,1,14,5,903,900,1005,1003,US,227,NA,62,63,NA,2,3,CLT,ORF
1994,1,15,6,859,900,1004,1003,US,227,NA,65,63,NA,1,-1,CLT,OR
1994,1,17,1,859,900,955,1003,US,227,NA,56,63,NA,-8,-1,CLT,OR
1994,1,18,2,904,900,959,1003,US,227,NA,55,63,NA,-4,4,CLT,ORF
1994,1,19,3,858,900,947,1003,US,227,NA,49,63,NA,-16,-2,CLT,O
1994,1,20,4,923,900,1015,1003,US,227,NA,52,63,NA,12,23,CLT,O
1994,1,21,5,NA,900,NA,1003,US,227,NA,NA,63,NA,NA,NA,CLT,ORF,
1994,1,22,6,859,900,1001,1003,US,227,NA,62,63,NA,-2,-1,CLT,O
1994,1,24,1,901,900,1006,1003,US,227,NA,65,63,NA,3,1,CLT,ORF
1994,1,25,2,859,900,952,1003,US,227,NA,53,63,NA,-11,-1,CLT,O
1994,1,27,4,910,900,1008,1003,US,227,NA,58,63,NA,5,10,CLT,OR
1994,1,28,5,858,900,1020,1003,US,227,NA,82,63,NA,17,-2,CLT,O
1994,1,29,6,859,900,953,1003,US,227,NA,54,63,NA,-10,-1,CLT,O


Use the textscan function to import the data, as it gives the flexibility needed to import particular columns from a file of mixed text of data; it also can contend with "NA" as an indicator of a missing value. Extract columns 1, 2, 4, and 5 (Year, Month, DayOfWeek, DepTime) from one file on the local machine.

tic
srcFile = fullfile(dataFolder, files(8).name);
f = fopen(srcFile);
C = textscan(f, '%f %f %*f %f %f %*[^\n]', 'Delimiter', ',', ...
fclose(f);
[Year, Month, DayOfWeek, DepTime] = C{:};
sprintf('%g seconds to load "%s".\n', toc, srcFile)

ans =



#### Open a Pool of Local Workers

The previous step probably took a bit of time to import, maybe a full minute. Luckily, this operation will run in parallel later. If necessary, open an 8-worker pool now:

try
parpool('local', 8);
catch
end

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


#### Load One File to All Workers

Move the textscan code to the cluster, and have all of the nodes in the cluster load this file. This is done with a single program/multiple data (SPMD) block. An spmd block of MATLAB code is executed on all eight of the workers in parallel:

tic
spmd
srcFile = fullfile(dataFolder, files(8).name);
f = fopen(srcFile);
C = textscan(f, '%f %f %*f %f %f %*[^\n]', 'Delimiter', ',', ...
[Year, Month, DayOfWeek, DepTime] = C{:};
fclose(f);
end
sprintf('%g seconds to load "%s" on all workers.\n', toc, srcFile{1})

ans =

65.2914 seconds to load "/Volumes/HD2/airline/data/1994.csv" on all workers.



#### Load a Different File on Each Worker

This previous step takes about the same length of time as the previous load, a near linear speed-up. But, observe that the same file was loaded eight times. To allow for different behaviors on different workers, each worker is assigned a unique index, called labindex. Use labindex to make each worker behave differently within an spmd block. In this case, by simply replacing the hard-coded index files(8) with files(labindex), the spmd block loads different files on different workers (recall that files is a vector of all filenames in the data set) -- this spmd block is otherwise identical to the previous one.

tic
spmd
srcFile = fullfile(dataFolder, files(labindex).name);
f = fopen(srcFile);
C = textscan(f, '%f %f %*f %f %f %*[^\n]', 'Delimiter', ',', ...
[Year, Month, DayOfWeek, DepTime] = C{:};
fclose(f);
end
sprintf('%g seconds to load different files on each worker.\n', toc)

ans =

61.4875 seconds to load different files on each worker.



#### Performing Non-cooperative Computation per Worker

Each worker now has its own set of data. In some applications, this may be enough, and each data set can be worked on independently from each other. This is sometimes referred to as an "embarrassingly parallel" problem. In the case of this data set, each file is data for an entire year, so calculations by year are trivial to compute. The value of Year should be identical within each worker. Confirm this in an spmd block, and then bring the results of that computation back locally for examination. Note the use of cell array-style indexing immediately after the spmd block; this is the means to retrieve all values a variable has on each worker. Since there are eight workers, eight values are returned -- this is referred to as a Composite in MATLAB.

spmd
myYear = Year(1);
end

[myYear{:}]

ans =

1     1     1     1     1     1     1     1

ans =

Columns 1 through 6

1987        1988        1989        1990        1991        1992

Columns 7 through 8

1993        1994



#### Visualize Departure Times

Create a histogram of departure hour over eight years to get a first real look at the data. Note that one year (1987) has significantly less magnitude than the other years -- this is because it is only a partial year's worth of data. A course-grain partitioning of your data using files and folder hierarchy can be a straightforward way of segmenting your data set into meaningful (and manageable) subsets.

spmd
H = hist(DepTime, 24);
end

plot(cell2mat(H(:))', '-o')
title('Histogram of departure times by year')
xlabel('Hour of the day')
ylabel('Number of flights')
legend(arrayfun(@(x) sprintf('%d', x), [myYear{:}], 'UniformOutput', false), ...
'Location', 'NorthWest')


#### Perform Cooperative Computation per Worker

But suppose analysis must be performed across all observations in the data set. For this, distributed arrays can be used. A distributed array allows MATLAB to work with an array spread across many workers as-if it were a single, "normal" MATLAB array. distributed arrays are supported by a large set of MATLAB functions, including linear algebra.

To create a distributed array from the data already in memory, first determine how big each worker's piece of the array is, and its total size. Calculate the individual and total lengths, and then create four distributed arrays, one for each column in the original data set, using a so-called codistributor (called a co-distributor because this is from the point of view of each worker; each needs to work cooperatively with other workers).

tic
spmd
lengths = size(Year,1);
end

lengths = cell2mat(lengths(:));
lengths = lengths';
totalLength = sum(lengths);

spmd
codistr = codistributor1d(1, lengths, [totalLength, 1]);
Year = codistributed.build(Year, codistr);
Month = codistributed.build(Month, codistr);
DayOfWeek = codistributed.build(DayOfWeek, codistr);
DepTime = codistributed.build(DepTime, codistr);
end

sprintf('%g seconds to create distributed arrays from Composites.\n', toc)

ans =

1.86178 seconds to create distributed arrays from Composites.



#### Visualize Year Across the Entire Data Set

There are now four variables (Year, Month, DayOfWeek, and DepTime) that contain the entire 8-year data set; these variables can be used much like normal MATLAB variables. This example is using a 1-D array, but multidimensional distributed arrays are also supported should you need to distribute "tiles" of data across a cluster. If you ever want to gather up the distributed array into one normal MATLAB array, you can do so with the gather function. Be careful when doing this, and ensure that your local computer has enough memory to hold the aggregate content of the entire distributed array. To create a histogram of departure year:

firstYear = gather(min(Year));
lastYear = gather(max(Year));
figure
hist(Year, firstYear:lastYear, 'EdgeColor', 'w')
title('Histogram of year')
xlabel('Year')
ylabel('Number of flights')


Note the use of gather here. Any operation on a distributed array produces another distributed array. Even though min and max return small arrays, they must nevertheless be gather'ed back to your local MATLAB before using them in the second argument to the hist function. You can see more starkly now that 1987 has far fewer observations than the subsequent years.

#### Visualize Departure Time Across the Entire Data Set

Next produce a histogram of the departure time. There is no need to gather anything this time because the second argument to hist (24, as in hours in a day) is a constant value and not something that needs to be computed from the distributed data.

figure;
hist(DepTime,24);
title('Histogram of departure time')
xlabel('Hour of Day')
ylabel('Number of flights')


#### Use 24 Histogram Bins

Readers will notice that the histogram bins run up to 2500, not 25. This is because the original data set expressed time as an integer encoding HHMM. To consider only the hour:

hist(floor(DepTime/100),24)
title('Histogram of departure time (Hour only)')
xlabel('Hour of Day')
ylabel('Number of flights')


#### Re-express the Departure Time

In fact, re-express HHMM as a number between 0 and 24, with the value after the decimal point expressing each minute as 1/60th. Update the value of the distributed array and then again plot the histogram with DepTime's improved representation.

clear C;
DepTime = floor(DepTime/100)+mod(DepTime,60)/60;
hist(DepTime,24)
title('Histogram of departure time (Fractional time)')
xlabel('Hour of Day')
ylabel('Number of flights')


Note the use of clear to rid the workspace of C, a temporary created earlier in the code when the file was imported. MATLAB generally avoids keeping unnecessary copies of data around, and keeping the temporary had no material effect on the program up till now. However, since the value of DepTime will change, MATLAB will be in a position where it needs to keep similar-but-distinct copies of the arrays. This is to be avoided for obvious reasons, so clear these older "snapshots" before modifying DepTime.

#### Removing Missing Data Points

Up to this point, the data has not been checked for the presence of missing data, represented by not-a-number (NaN) in the departure time. Look for NaNs, and then remove those bad rows from all of the distributed arrays.

tic
goodRows = ~isnan(DepTime);
fprintf('There are %d rows to remove from the data.\n', ...
totalLength-gather(sum(goodRows)));
Year = Year(goodRows);
Month = Month(goodRows);
DayOfWeek = DayOfWeek(goodRows);
DepTime = DepTime(goodRows);

sprintf('%g seconds to trim bad rows.\n', toc)

There are 419397 rows to remove from the data.

ans =

23.9631 seconds to trim bad rows.



#### Work With the Full Data Set

It is time to work with all of the files in the data set. Because the data files are text, this application will spend most of its runtime parsing text and converting to numeric. This may be acceptable for a one-time analysis, but if you anticipate the need to load the data set more than once, it is likely worth the effort to save that numeric representation for future use. This is done in a simple parfor loop, creating a MATLAB MAT-File for each file in the original data set. This loop takes a few minutes to run the first time. However, the code is also smart enough to not recreate a MAT-File when it already exists, so this loop finishes quickly when later rerun.

Note: If you have a cluster with more than eight cores available to you, now would be a good time to close the eight-worker pool and open a larger pool. If you do not have access to a cluster, that is okay too as long as your local computer has around 16 GB of RAM or more.

clear all
dataFolder = '/Volumes/HD2/airline/data';
files = dir(fullfile(dataFolder, '*.csv'));
files = files(~strncmp({files(:).name}, '.', 1));  % Remove "dot" files (Mac)

tic
parfor i=1:numel(files)
% See http://www.mathworks.com/support/solutions/en/data/1-D8103H for the trick below
parSave = @(fname, Year, Month, DayOfWeek, DepTime) ...
save(fname, 'Year', 'Month', 'DayOfWeek', 'DepTime');

srcFile = fullfile(dataFolder, files(i).name);
[path, name, ~] = fileparts(srcFile);
cacheFile = fullfile(path, [name '.mat']);
if isempty(dir(cacheFile))
f = fopen(srcFile);
C = textscan(f, '%f %f %*f %f %f %*[^\n]', 'Delimiter', ',', ...
fclose(f);
[Year, Month, DayOfWeek, DepTime] = C{:};
parSave(cacheFile, Year, Month, DayOfWeek, DepTime);
end
end

sprintf('%g seconds to import text files to save to MAT-Files.\n', toc)

ans =

220.984 seconds to import text files to save to MAT-Files.



#### Load the Entire Data Set

Now it is time to load the entire data set, which was just saved into MAT-Files. Since the number of workers may be fewer than the number of files, each worker will be responsible for a collection of files. Files are assigned to workers using their labindex and modulo arithmetic. Each worker reads all of its files, and vertcat's their contents into a single array. While importing the data, this is a convienent time to perform the processing introduced above -- filtering NaN rows, converting the HHMM representation of departure time, and noting the length of each worker's piece of the data set. Next, create the distributed arrays. Though this section of code is the longest of this post, it is largely a copy-and-paste of code already discussed.

tic
spmd
myIndicies = labindex:numlabs:numel(files);

Year = cell(size(myIndicies));
Month = cell(size(myIndicies));
DayOfWeek = cell(size(myIndicies));
DepTime = cell(size(myIndicies));

for i = 1:numel(myIndicies)
srcFile = fullfile(dataFolder, files(myIndicies(i)).name);
[path, name, ~] = fileparts(srcFile);
cacheFile = fullfile(path, [name '.mat']);
Year{i} = tmpT.Year;
Month{i} = tmpT.Month;
DayOfWeek{i} = tmpT.DayOfWeek;
DepTime{i} = tmpT.DepTime;
end
end

clear tmpT  % As above, avoid deep copies

spmd
Year = vertcat(Year{:});
Month = vertcat(Month{:});
DayOfWeek = vertcat(DayOfWeek{:});
DepTime = vertcat(DepTime{:});

goodRows = ~isnan(DepTime);
Year = Year(goodRows);
Month = Month(goodRows);
DayOfWeek = DayOfWeek(goodRows);
DepTime = DepTime(goodRows);

DepTime = floor(DepTime/100)+mod(DepTime,60)/60;
lengths = size(Year,1);
end

lengths = cell2mat(lengths(:));
lengths = lengths';
totalLength = sum(lengths);

spmd
codistr = codistributor1d(1, lengths, [totalLength, 1]);
Year = codistributed.build(Year, codistr);
Month = codistributed.build(Month, codistr);
DayOfWeek = codistributed.build(DayOfWeek, codistr);
DepTime = codistributed.build(DepTime, codistr);
end

sprintf('%g seconds to import from MAT-Files and create distributed arrays.\n', toc)

ans =

12.121 seconds to import from MAT-Files and create distributed arrays.



#### Rerun Visualizations

You can now perform some of the visualizations over the entire data set. DayOfWeek is visualized this time.

figure
hist(DayOfWeek, 1:7)
title('Histogram of day of the week')
xlabel('Day of the week')
ylabel('Number of flights')
set(gca, 'XTickLabel', {'Mon', 'Tues', 'Wed', 'Thurs', 'Fri', 'Sat', 'Sun'});


#### Rebalance Data Among the Pool of Workers

Load balancing is the final topic of discussion. The distributed array is partitioned at the file level, with zero or more files per worker. This is probably fine if your input files are roughly the same size and the number of workers is smaller than the number of files. However, uneven file size may create an imbalance of work, or, even worse, more workers than files will result in completely idle workers! If you are running more than 22 workers, you will see this effect when working with this data set.

It therefore can be helpful to rebalance the distribution of a distributed array. This will require some transfer of data between workers, but it may be a worthwhile investment to accelerate subsequent computations. Calculate the "ideal" (balanced) distribution, create a new codistributor with that ideal distribution, and then redistribute using that ideal.

tic
spmd
beforeSizes = size(getLocalPart(Year), 1);
end

newBreaks=round(linspace(0, totalLength, numel(lengths)+1));
newLengths=newBreaks(2:end)-newBreaks(1:end-1);

spmd
newCodistr = codistributor1d(1, newLengths, [totalLength, 1]);
Year = redistribute(Year, newCodistr);
Month = redistribute(Month, newCodistr);
DayOfWeek = redistribute(DayOfWeek, newCodistr);
DepTime = redistribute(DepTime, newCodistr);
end

spmd
afterSizes = size(getLocalPart(Year), 1);
end

sprintf('%g seconds to redistribute the arrays.\n', toc)

figure
bar([[beforeSizes{:}]', [afterSizes{:}]'])
title('Data per worker (Before and After)')
xlabel('Worker Number')
ylabel('Amount of data')
legend({'Before', 'After'}, 'Location', 'NorthWest')

ans =

14.1069 seconds to redistribute the arrays.



#### In Closing

Hopefully this paper and code has given you a sense of the kind of interactive workflow that is possible with big data when using a distributed array. Parallel file I/O, distributed and cooperative operations, visualization, and load balancing are introduced and explored. Though the examples used here were simple, these basic principles extend to other data sets similarly distributed among a collection of files.

Significantly, this post only scratched the surface of the actual numeric capabilities of distributed arrays -- higher order operations such as single-value decomposition (SVD) and other linear algebra functions, FFTs, and so on are also available. For further reading about distributed arrays, try the following:

How do you manage your large data analysis? Let us know here.

Get the MATLAB code

Published with MATLAB® R2013b

## Function is as functiontests

Today we have a post from guest blogger Andy Campbell to highlight a function-based approach for writing tests in MATLAB. We will be demonstrating how you can write simple functions that aid you in designing a physical system. This approach allows you to iterate quickly and improve your designs while protecting you from violating requirements met by earlier solutions.

With that, let's jump right in!

### Contents

#### Problem Statement

Suppose you are working on a design problem that is governed by a second order differential equation. The beauty of mathematics is that this actually has a wide range of applicability across many different areas. For our example, let's consider the design of a simple mass-spring-damper. Such a system may be used to model the suspension system on a car or truck. This is often represented using the following simplified model:

Where:

• $m$ = mass
• $k$ = spring constant
• $c$ = damping coefficient
• $x$ = displacement from equilibrium

Assuming there are no external forces, the equation of motion for this system is the solution to the second order differential equation $m\ddot{x} + c\dot{x} + kx = 0$. Following a similar approach to an earlier blog post, we can prepare this equation to be solved by one of MATLAB's ode solvers and write the following function used to evaluate our design:

type simulateSystem

function [x, t] = simulateSystem(design)

% Design variables
c = design.c;
k = design.k;

% Constant variables
z0 = [-0.1; 0];  % Initial Position and Velocity
m = 1500;        % Mass

odefun = @(t,z) [0 1; -k/m -c/m]*z;
[t, z] = ode45(odefun, [0, 3], z0);

% The first column is the position (displacement from equilibrium)
x = z(:, 1);


In this function, $z$ is a vector representing the state of the system. It contains $x$ and its derivative $\dot{x}$. Note that the initial conditions are that the mass has a position of -0.1 and its velocity is 0. This state may model the case where you are driving on the highway and hit a pothole. The initial position models the mass at the bottom of the pothole and the zero initial velocity models the suspension compressed, stopped, and about to recoil. For the sake of this example, let's say that the mass is constant at 1500 kg, which is reasonable for the load carried by a single wheel of a large light duty truck.

The design variables are the spring constant, $k$, and the damping coefficient, $c$. Initially they both have an arbitrary value of 1.

type springDamperDesign0

function design = springDamperDesign0
% Design variables
design.c = 1;          % Damping Coefficient (initial design)
design.k = 1;          % Spring Constant (initial design)


#### Enter functiontests

If I were tasked to design a system like this I would undoubtedly look for some requirements of the design and ensure my design meets them all. However, a problem can surface if I implement a design that fits one requirement but then adjust the design to fit a second requirement that ends up breaking the first. For one or two requirements this is fairly easy to manage but it quickly becomes intractable with many requirements or if new requirements are added at a later time. Wouldn't it be great if we had some simple mechanism to capture the first requirement to ensure it is never violated? We do! Let's write a test!

An easy way to start writing your own test function is to use the functiontests function. The functiontests function can be called at the top of your test function file, and it allows you to create and group many related tests in the same file.

type testTemplate

function tests = testTemplate

tests = functiontests(localfunctions);

end


Let's dissect what is happening here. By convention, the name of this file starts or ends with the word "test" (case insensitive) to help identify it as a function file that contains tests. All it needs to do is call the localfunctions function as input into the functiontests function, and then assign the result to the output of the main function.

#### What is functiontests?

functiontests is a MATLAB function (new in R2013b) that accepts a cell array of function handles to local functions within a MATLAB file. It then determines which of these local functions are tests based on naming convention and creates a Test element for each of them. All local functions that start or end with the word "test" (case insensitive) are considered tests.

#### What is localfunctions?

The localfunctions function is another MATLAB function that is new in R2013b. This function takes no input arguments and simply returns a cell array of handles to all of the local functions within the calling file. When testing with functiontests, using localfunctions is highly recommended so that when you add additional test functions to the file their function handles do not need to be manually added to a cell array at the top of the file. Instead, localfunctions automatically discovers these new tests and adds them to the cell array of function handles passed to functiontests.

#### The First Test

With this basic framework in place we can write our first test. The first design requirement might be that the system returns to equilibrium within a certain amount of time, for example, 2 seconds. Let's write a test to see how we are doing against this requirement.

type design0Test

function tests = design0Test

tests = functiontests(localfunctions);

end

function testSettlingTime(testCase)
% Test that the system settles to within 0.001 of zero under 2 seconds.

[position, time] = simulateSystem(springDamperDesign0);

positionAfterSettling = position(time > 2);

% For this example, verify the first value after the settling time.
verifyEqual(testCase, positionAfterSettling(1), 0, 'AbsTol', 0.001);
end


This test:

1. Simulates the system using the specified design variables.
2. Captures the position after the required settling time.
3. Verifies that the position reaches 0 within +/- 0.001 by the required time. Typically for a real test you would want to verify that all positions after the settling time are sufficiently close to zero, not just the point at the desired time.

Note that this test function uses, and requires, a single testCase argument. This argument contains important test context and enables the use of a rich library of qualification functions such as verifyEqual.

What happens when we call this function?

test = design0Test;
disp(test)

  Test with properties:

Name: 'design0Test/testSettlingTime'
SharedTestFixtures: []


You can see a single test was created from this function. To run this test, simply call the run function with this test as an argument.

run(test);

Running design0Test

================================================================================
Verification failed in design0Test/testSettlingTime.

---------------------
Framework Diagnostic:
---------------------
verifyEqual failed.
--> NumericComparator failed.
--> The values are not equal using "isequaln".
--> AbsoluteTolerance failed.
--> The value was not within absolute tolerance.

Tolerance Definition:
abs(expected - actual) <= tolerance
Tolerance Value:
1.000000000000000e-03

Actual Value:
-0.099863405108097
Expected Value:
0

------------------
Stack Information:
------------------
In H:\Documents\LOREN\MyJob\Art of MATLAB\Andy Campbell\informal test interface\design0Test.m (testSettlingTime) at 15
In C:\Program Files\MATLAB\R2013b\toolbox\matlab\testframework\+matlab\+unittest\FunctionTestCase.m (FunctionTestCase.test) at 90
================================================================================
.
Done design0Test
__________

Failure Summary:

Name                          Failed  Incomplete  Reason(s)
===========================================================================
design0Test/testSettlingTime    X                 Failed by verification.



As you can see, the test failed and printed some information from verifyEqual to help diagnose the failure. Let's see how close we are graphically. Note that, for this example, in order to show the dynamics we need to show a much larger time span than our system should require. The simulateSystemOverLargeTime contains the same content as simulateSystem defined above but it has a different time span.

[x, t] = simulateSystemOverLargeTime(springDamperDesign0);
plot(t, x);
title('Undesigned Response');
xlabel('Time');
ylabel('Position')


This plot demonstrates how this design doesn't come close reaching equilibrium under 2 seconds (it is not even close after 1000 seconds). However, this is to be expected since we haven't designed anything yet.

#### The First Design

Assume we come up with the following design:

type springDamperDesign1

function design = springDamperDesign1
design.k = 5e6; % Spring Constant
design.c = 1e4; % Damping Coefficient


Does this design meet our criteria?

run(design1Test);

Running design1Test
.
Done design1Test
__________



Yes it does, we've passed! Let's take a look at our response:

[x, t] = simulateSystem(springDamperDesign1);
plot(t, x)
title('Design #1')
xlabel('Time')
ylabel('Position')


Great, we reach equilibrium rapidly enough. However, unfortunately it doesn't look like a comfortable ride in your brand new truck. It looks like all the vibration from hitting a pothole will shake the poor passengers' eyes right out of their sockets. In addition, the system has overshot the equilibrium point by 80% of the pothole depth (0.08 maximum position with a pothole depth of 0.1). We can do better. Let's write a test to prove it.

#### The Second Test

The missing requirement here could be a restriction on how far the position can overshoot the equilibrium point. We can now add a second test to our arsenal:

type design1SecondTest

function tests = design1SecondTest

tests = functiontests(localfunctions);

end

function testSettlingTime(testCase)
% Test that the system settles to within 0.001 of zero under 2 seconds.

[position, time] = simulateSystem(springDamperDesign1);

positionAfterSettling = position(time > 2);

% For this example, verify the first value after the settling time.
verifyEqual(testCase, positionAfterSettling(1), 0, 'AbsTol', 0.001);
end

function testOvershoot(testCase)
% Test to ensure that overshoot is less than 0.01

position = simulateSystem(springDamperDesign1);
overshoot = max(position);

verifyLessThan(testCase, overshoot, 0.01);
end


This overshoot test uses verifyLessThan to ensure that the position never moves too far above the equilibrium point. We would expect this test to fail with our current design. First, we again create the Test array. Note that the second test was picked up automatically and the Test array now contains two elements instead of one.

tests = design1SecondTest;
disp(tests)
run(tests);

  1x2 Test array with properties:

Name
SharedTestFixtures

Running design1SecondTest
.
================================================================================
Verification failed in design1SecondTest/testOvershoot.

---------------------
Framework Diagnostic:
---------------------
verifyLessThan failed.
--> The value must be less than the maximum value.

Actual Value:
0.082943282378938
Maximum Value (Exclusive):
0.010000000000000

------------------
Stack Information:
------------------
In H:\Documents\LOREN\MyJob\Art of MATLAB\Andy Campbell\informal test interface\design1SecondTest.m (testOvershoot) at 24
In C:\Program Files\MATLAB\R2013b\toolbox\matlab\testframework\+matlab\+unittest\FunctionTestCase.m (FunctionTestCase.test) at 90
================================================================================
.
Done design1SecondTest
__________

Failure Summary:

Name                             Failed  Incomplete  Reason(s)
==============================================================================
design1SecondTest/testOvershoot    X                 Failed by verification.



The overshoot test fails, and this is not a bad thing! Not only does it inform us that we are missing this requirement, it also serves as a sanity check to confirm that our tests are indeed doing what we think they are doing.

#### The Second Design

Based on the response of the system, it seems our current design is underdamped, so why don't we just add some damping to prevent our overshoot and run the tests again.

type springDamperDesign2

function design = springDamperDesign2

design.k = 5e6;   % Spring Constant
design.c = 2.5e6; % Increase the Damping Coefficient from 1e4


Let's see how this does

[x, t] = simulateSystem(springDamperDesign2);
plot(t, x)
title('Design #2')
xlabel('Time')
ylabel('Position')


Looking at the response of this system it seems like a nice smooth ride that gently moves the truck back to equilibrium after the disturbance. Let's go with it! First let's just make sure it passes our tests.

run(design2Test);

Running design2Test

================================================================================
Verification failed in design2Test/testSettlingTime.

---------------------
Framework Diagnostic:
---------------------
verifyEqual failed.
--> NumericComparator failed.
--> The values are not equal using "isequaln".
--> AbsoluteTolerance failed.
--> The value was not within absolute tolerance.

Tolerance Definition:
abs(expected - actual) <= tolerance
Tolerance Value:
1.000000000000000e-03

Actual Value:
-0.001823826310161
Expected Value:
0

------------------
Stack Information:
------------------
In H:\Documents\LOREN\MyJob\Art of MATLAB\Andy Campbell\informal test interface\design2Test.m (testSettlingTime) at 15
In C:\Program Files\MATLAB\R2013b\toolbox\matlab\testframework\+matlab\+unittest\FunctionTestCase.m (FunctionTestCase.test) at 90
================================================================================
..
Done design2Test
__________

Failure Summary:

Name                          Failed  Incomplete  Reason(s)
===========================================================================
design2Test/testSettlingTime    X                 Failed by verification.



Look at that! This design passes the overshoot test, but we have regressed and, once again, failed the settling time test. This is where we see the value in encoding these requirements through testing. There was some point after we secured the settling time requirement that we changed our focus to the overshoot requirement. It ended up being way too easy to evolve our design when tuning the overshoot and lose sight of the settling time. Focusing on overshoot, the latest design seems so reasonable, and since it is close to meeting the settling time requirement it was hard to notice through manual inspection. Obviously, close is not good enough! These tests are now protecting us from violating previously met requirements. These tests can be used to lockdown all past requirements, even when they are no longer in focus. Adding more tests can only make these requirements stricter or more complete.

This is especially important as the number of requirements for your system grows. You can become more protected by a larger and more comprehensive set of tests. You might imagine a case where you adjust the design in one way and 1 test out of 100 fails, but then you adjust the design in a different way and 75 tests out of 100 fail. When MATLAB is used to execute all of the tests in a single call you are able to gain insight into such cases, suggesting that the former design is likely closer to meeting your full set of requirements than the latter.

Since these tests use the standard MATLAB test framework they can be shared with colleagues and/or customers who can easily run them and verify that they capture the requirements sufficiently. The tests act as an executable specification of the requirements, as well as the validation that these requirements are met.

#### The Final Design

For completeness, let's go ahead and make this test pass and then we can debrief. With a second order system we can design for critical damping and run the tests again.

type springDamperDesign3.m
result = run(design3Test);
[x, t] = simulateSystem(springDamperDesign3);
plot(t, x)
title('Design #3')
xlabel('Time')
ylabel('Position')

function design = springDamperDesign3
m = 1500; % Need to know the mass to determine critical damping

design.k = 5e6;                   % Spring Constant
design.c = 2*m*sqrt(design.k/m);  % Damping Coefficient to be critically damped

Running design3Test
..
Done design3Test
__________



The tests pass and the response looks as expected for critical damping. Now we can move forward improving our design with the knowledge that we are protected against failing to meet our first two requirements.

#### Conclusion

Through a simplified example, we have demonstrated the power in using tests to validate engineering and scientific problems. All it took was writing a MATLAB function containing a couple of structured local functions. What's next? Well, for this design there might be any number of next steps, such as:

1. Testing to ensure the spring constant is low enough to have a smoother, less jerky response.
2. Testing to limit the cost of the springs and/or dampers (perhaps stiffer springs cost more money).
3. Testing to ensure the natural frequency of the system lies within a certain range.
4. Improving the model. In the real world an accurate model is likely much more complicated than a canonical second order system.
5. Testing the design against a payload. For example, if the mass needs to operate with cargo with additional mass within the range of 0 to 500 kg.
6. Testing the shape and characteristics of the frequency response.

The tests we wrote today are testing only the results of a given design. Remember the underlying system does not need to be a mass-spring-damper at all. The design and/or model can completely change, but these tests still are available to ensure that the response of the system meets the requirements.

The initial investment of writing a test up front yields lasting dividends and provides a safety net for analysis and design. It is interesting that these activities are what we, as scientists and engineers, already do anyway. We always test our designs, we just don't always capture the process in the form of an executable test for future execution. This approach gets easier and easier to do the more you get the hang of it. Furthermore, any time invested in such tests up front quickly pays off in productivity gains by allowing safe, rapid design iterations.

This example only begins to scratch the surface of MATLAB testing capabilities. For more information, take a look at the documentation for the MATLAB Unit Test Framework, which includes this approach using test functions. Let us know what you think here!

Get the MATLAB code

Published with MATLAB® R2013b

## Timing Code

and Revisiting Fibonacci Festival from 2006

I originally wrote about Fibonacci numbers in 2006. Today, along with Sarah Wait Zaranek, we'd like to discuss some of those same ideas, but focus only on execution time. In R2013b, MATLAB includes a new timing function, timeit, one that I mentioned in some earlier posts when it was only available on the MATLAB File Exchange.

### Contents

#### Using timeit

To use timeit, you simply pass in a function handle that requires no inputs. If your function requires input values, you can pass them in by creating an anonymous function, where the values are captured from the current workspace. You can learn more about function handles in the MATLAB documentation, and find examples and discussions on this blog.

timeit runs the code several times to ensure that any first-run timing effects do not affect the results - and returns the median time. If you wish to time your function in the situation where it returns more than 0 or 1 output values, you can specify that as one of the inputs to timeit.

#### Four Algorithms

In the earlier Fibonacci post, I introduce and times four algorithms. Today we will use the same algorithms, coded as separate programs, in conjunction with timeit. Since we have four algorithms, we will initialize a vector for collecting the timing information.

fibTimes = zeros(1,4);


#### Straightforward for-loop

Let's calculate 102 Fibonacci values (since the first 2 are known already, [1, 1]).

nf = 102;
type fibloop

function fibs = fibloop(nf)
fibs(1) = 1;
fibs(2) = 1;
for n = 3:nf
fibs(n) = fibs(n-1)+fibs(n-2);
end


To use timeit, we simply give it a function handle to the function we wish to time.

fibTimes(1) = timeit(@()fibloop(nf));


#### for-loop with Preallocation

You've all heard it before. One reason the former calculation took so much time is at least alleged to be because of growing the array fibf each time through the loop. So now let's use the same loop, but preallocate the output array so it doesn't need to grow during the calculation.

type fibloopPrealloc

function fibs = fibloopPrealloc(nf)
fibs = ones(1,nf);
for n = 3:nf
fibs(n) = fibs(n-1) + fibs(n-2);
end



Now time it.

fibTimes(2) = timeit(@()fibloopPrealloc(nf));


#### Recursive Function

We can use recursion to calculate successive terms, as I showed in that earlier post. Here's the code.

type fibRecursive

function fibs = fibRecursive(n)
if n == 1,
fibs = 1;   % First element is 1.
return;
elseif n == 2
fibs = [1 1];  % First two elements are 1.
else
% Call fibrec with previous result and compute next one from it.
fm1 = fibRecursive(n-1);
fibs = [fm1 fm1(end)+fm1(end-1)];
end


And now let's use it to calculate the first chunk of Fibonacci numbers.

fibTimes(3) = timeit(@()fibRecursive(nf));


#### Use filter to Calculate the Sequence

I also explained in the earlier Fibonacci post how to use filter to calculate the sequence as well. Here's the code.

type fibFilter
fibTimes(4) = timeit(@()fibFilter(nf));

function fibf = fibFilter(n)
x = [1 zeros(1,n-1)];
a = [1 -1 -1];
b = 1;
fibf = filter(b, a, x);


#### Compare Computation Times

It's time to see how quickly the various methods performed. We have multiplied the times by 1000 to make reading the numbers easier.

fprintf('%s\n\n','Compare Times (in milliseconds)')
meths = ['loop           ',...
'prealloc       ',...
'recursive      ',...
'filter         '];
fprintf('%s\n\n',meths)
fprintf('%2.9f    ', fibTimes*1000), fprintf('\n')

Compare Times (in milliseconds)

loop           prealloc       recursive      filter

0.040766819    0.010003257    0.517234655    0.010549081


#### Do You Do Lots of Timing?

Do you do lots of timing? How do you do it? Profiler? Will timeit help? Let us know here.

Get the MATLAB code

Published with MATLAB® R2013b

## Introduction to the New MATLAB Data Types in R2013b

Today I’d like to introduce a fairly frequent guest blogger Sarah Wait Zaranek who works for the MATLAB Marketing team here at The MathWorks. She and I will be writing about the new capabilities for MATLAB in R2013b. In particular, there are two new data types in MATLAB in R2013b – table and categorical arrays.

### Contents

#### What are Tables and Categorical Arrays?

Table is a new data type suitable for holding heterogenous data and metadata. Specifically, tables are useful for mixed-type tabular data that are often stored as columns in a text file or in a spreadsheet. Tables consist of rows and column-oriented variables. Categorical arrays are useful for holding categorical data - which have values from a finite list of discrete categories.

One of the best ways to learn more about tables and categorical arrays is to see them in action. So, in this post, we will use tables and categoricals to examine some airplane flight delay data. The flight data is freely available from the Bureau of Transportation Statistics (BTS). You can download it yourself here. The weather data is from the National Climatic Data Center (NCDC) and is available here.

#### Importing Data into a Table

You can import your data into a table interactively using the Import Tool or you can do it programmatically, using readtable.

FlightData = readtable('Jan2010Flights.csv');

whos FlightData

  Name                Size              Bytes  Class    Attributes

FlightData      17816x7             9007742  table



The entire contents of the file are now contained in a single variable – a table. Here you are reading in your data from a csv file. readtable also supports reading from .txt,.dat text files and Excel spreadsheet files. Tables can also be created directly from variables in your workspace.

#### Looking at Variable Names (Column Names)

You can see all the variable names (column names in our table) by looking at the VariableNames properties.

FlightData.Properties.VariableNames

ans =
Columns 1 through 5
'FL_DATE'    'CARRIER'    'ORIGIN'    'DEST'    'CRS_DEP_TIME'
Columns 6 through 7
'DEP_TIME'    'DEP_DELAY'


This particular table does not contain any row names, but for a table with row names you can access the row names using the RowNames property.

#### Accessing the Data in Your Table

There are multiple ways to access the data in your table. You can use dot indexing to access or modify a single table variable, similar to how you use fieldnames in structures.

For example, using dot indexing you can plot a histogram of the departure delays (in minutes).

hist(FlightData.DEP_DELAY)
title('Histogram of Flight Delays in Minutes')


You can also display the first 5 departure delays.

FlightData.DEP_DELAY(1:5)

ans =
49
-7
-5
-8
-10


You can also extract data from one or more variables in the table using curly braces. Within the curly braces you can use numeric indexing or variable and row names. For example, you can extract the actual departure times and scheduled depature times for the first 5 flights.

SomeTimes = FlightData{1:5,{'DEP_TIME','CRS_DEP_TIME'}};
disp(SomeTimes)

        1149        1100
1053        1100
1055        1100
1052        1100
1050        1100


This is similar to indexing with cell arrays. However, unlike with cells, this concatenates the specified variables into a single array. Therefore, the data types of all the specified variables need to be compatible for concatenation.

#### Converting Data to Categorical Arrays

You can convert some of the variables in your table using categorical. Categorical arrays are more memory efficent than holding cell arrays of strings when you have repeated data. Categorical arrays store only one copy of each category name, reducing the amount of memory required to store the array. You can use whos to see the amount of memory you save by converting the data to a categorical array.

whos FlightData

  Name                Size              Bytes  Class    Attributes

FlightData      17816x7             9007742  table


FlightData.ORIGIN = categorical(FlightData.ORIGIN);
FlightData.DEST = categorical(FlightData.DEST);
FlightData.CARRIER = categorical(FlightData.CARRIER);

whos FlightData

  Name                Size              Bytes  Class    Attributes

FlightData      17816x7             2857264  table



Using categories, you can find all the distinct categories in your array. By default, categorical arrays do not define a definite order. If your data contains categories with a definite order, you can set the 'Ordinal' flag to true when creating your categorical array. The default the order will be alphabetical, but you can prescribe you own order instead.

categories(FlightData.CARRIER)

ans =
'9E'
'AA'
'AS'
'B6'
'CO'
'DL'
'F9'
'FL'
'MQ'
'OH'
'UA'
'US'
'WN'
'XE'
'YV'


Categorical arrays are also faster and more convenient than cell arrays of strings for indexing and searching. By converting to categorical arrays, you can then mathematically compare sets of strings just like you would do with numeric values. You can use this functionality to create a new table containing only the flights that left from Boston.

#### Creating a New Table

You can create a new table from a section of an existing table using parentheses with numerical indexing, variable names, or row names. Since the flight origin is now a categorical array, you can use logical indexing to find all flights that left from Boston.

idxBoston = FlightData.ORIGIN == 'BOS' ;
BostonFlights = FlightData(idxBoston,:);

height(FlightData)
height(BostonFlights)

ans =
17816
ans =
8904


You can also modify your table by adding and removing variables and rows. All variables in a table must have the same number of rows, but they can be of different widths.

Let's add a new variable (DATE) to represent the serial date number for the various flight dates.

BostonFlights.DATE = datenum(BostonFlights.FL_DATE);


The origin now has all BOS values and you are not going to use destination information right now, so those variables can be removed. HOUR can also be calculated, as well as a LATE variable which indicates if the flight was 15 minutes late or more.

BostonFlights.ORIGIN = [];
BostonFlights.DEST = [];
BostonFlights.FL_DATE = [];

BostonFlights.HOUR = floor(BostonFlights.CRS_DEP_TIME./100);
BostonFlights{:,'LATE'} = BostonFlights.DEP_DELAY > 15;


#### Removing Missing Data

Tables have supported functions for finding and standardizing missing data. In this case, you can find any missing data using ismissing and remove it. You can use height, which gives you the number of table rows, to see how many flights were removed from the table. We exploit logical indexing to get only the flights that have no missing data.

height(BostonFlights)

ans =
8904

TF = any(ismissing(BostonFlights),2);
BostonFlights = BostonFlights(~TF,:);
height(BostonFlights)

ans =
8640


#### Summarizing a Table

You can then see descriptive statistics for each variable in this new table by using summary.

summary(BostonFlights)

Variables:
CARRIER: 8640x1 categorical
Values:
9E      85
AA     913
AS      55
B6    1726
CO     321
DL    1215
F9      26
FL     536
MQ     751
OH     341
UA     643
US    1494
WN     384
XE      93
YV      57
CRS_DEP_TIME: 8640x1 double
Values:
min        500
median    1215
max       2359
DEP_TIME: 8640x1 double
Values:
min          2
median    1224
max       2400
DEP_DELAY: 8640x1 double
Values:
min       -25
median     -3
max       419
DATE: 8640x1 double
Values:
min       7.3414e+05
median    7.3415e+05
max       7.3417e+05
HOUR: 8640x1 double
Values:
min        5
median    12
max       23
LATE: 8640x1 logical
Values:
true     1471
false    7169


#### Sorting Data

There are additional functions to sort tables, apply functions to table variables, and merge tables together. For example, you can sort your BostonFlights by departure delay.

BostonFlights = sortrows(BostonFlights,'DEP_DELAY','descend');
BostonFlights(1:10,:)

ans =
CARRIER    CRS_DEP_TIME    DEP_TIME    DEP_DELAY       DATE       HOUR
_______    ____________    ________    _________    __________    ____
DL         1830             129        419          7.3416e+05    18
DL         1850             111        381          7.3414e+05    18
CO         1755              10        375          7.3414e+05    17
AA         1710            2323        373          7.3416e+05    17
DL          630            1240        370          7.3416e+05     6
FL         1400            1951        351          7.3414e+05    14
FL         1741            2330        349          7.3414e+05    17
UA         1906              22        316          7.3414e+05    19
AA          840            1355        315          7.3416e+05     8
AA          905            1420        315          7.3414e+05     9

LATE
_____
true
true
true
true
true
true
true
true
true
true


#### Applying Functions to Table Variables

You can apply functions to work with table variables, with varfun.

varfun has optional additional calling inputs such as 'InputVariables' and 'GroupingVariables'. 'InputVariables' lets you specific which variables you want to operate on instead of operating on all the variables in your table. 'GroupingVariables' let you define groups of rows on which to operate. varfun would then apply your function to each group of rows within each of the variables of your table, rather than to each entire variable.

You can use varfun to calculate the mean delay for all flights and the fraction of late flights for a given hour on a given day. The default output of varfun is a table.

ByHour = varfun(@mean, BostonFlights, ...
'InputVariables', {'DEP_DELAY', 'LATE'},...
'GroupingVariables',{'DATE','HOUR'});

disp(ByHour(1:5,:))

                   DATE       HOUR    GroupCount    mean_DEP_DELAY    mean_LATE
__________    ____    __________    ______________    _________
734139_5    7.3414e+05    5        5               3.8                  0
734139_6    7.3414e+05    6       19            10.579            0.31579
734139_7    7.3414e+05    7       18            6.7778            0.11111
734139_8    7.3414e+05    8       21            8.8571            0.33333
734139_9    7.3414e+05    9       17            5.0588            0.23529


#### Joining (Merging) Tables

Weather might have an important role in determining if a flight is delayed. For a given hour, you might want to know both the delayed flight information and the weather at the airport. So, you can start by reading in another table containing weather data for Boston Logan Airport. Then, you can merge that table with the existing ByHour table.

Since there are a lot of variables in this file, you can specify the input format when using readtable. This allows you to use * to skip variables that you aren't interested in loading into the table. For more information about specifying formating strings, see here in the documentation. Since this data uses 'M' to represent missing data, you can use 'TreatAsEmpty' to replace any instances of 'M' with the standard missing value indicator (NaN for numeric values).

FormatStr = ['%*s%s%f' repmat('%*s',1,9) '%f' repmat('%*s',1,7) '%f',...
repmat('%*s',1,3),'%f', repmat('%*s',1,18)];

'Format',FormatStr,'TreatAsEmpty','M');

WeatherData.Properties.VariableNames

ans =
'Date'    'Time'    'DryBulbCelsius'    'DewPointCelsius'    'WindSpeed'


WeatherData contains the date, time, dew point and dry bulb temperature in Celsius, and wind speed. Let's convert DATE to a serial date number and round to the hour for the time measurement.

WeatherData.DATE = datenum(WeatherData.Date,'yyyymmdd');
WeatherData.Date = [];
WeatherData.HOUR = floor(WeatherData.Time/100);


Since there are multiple weather measurements per hour, you can average the data by hour using varfun.

ByHourWeather = varfun(@mean, WeatherData, ...
'InputVariables', {'DryBulbCelsius','DewPointCelsius','WindSpeed'},...
'GroupingVariables',{'DATE','HOUR'});


Now, you can merge the two tables using join which matches rows using key variables (columns) common to both tables. join keeps all the variables from the first input table and appends the corresponding variables from the second input table. The table that join creates will use the key variable values as the row names. In this case, that means the row names will represent the date and hour data.

AllData = join(ByHour,ByHourWeather,'Keys',{'DATE','HOUR'});


#### Plotting Data From the Final Table

Let's now plot this final data set to get an idea of the effect of weather on the flight delays.

AllData.TDIFF =  ...
abs(AllData.mean_DewPointCelsius - AllData.mean_DryBulbCelsius);

scatter(AllData.TDIFF, AllData.mean_DEP_DELAY,...
[], AllData.mean_DEP_DELAY,'filled');
xlabel('abs(DewPoint-Temperature)')
ylabel('Average Departure Delay')

scatter3( AllData.HOUR,AllData.TDIFF,AllData.mean_DEP_DELAY,...
[],AllData.mean_DEP_DELAY,'filled');
xlabel('Hour of Flight')
ylabel('abs(DewPoint-Temperature)')
zlabel('Average Departure Delay')


Qualitatively it looks having a temperature near the dew point (greater change of precipitation) effects the departure time. There are other factors at work as well, but it is nice to know that our intuition (flights later in the day and when it is cold and snowing might have a greater chance to be delayed) seems to work with the data.

Can you see yourself using tables and categorical arrays? Let us know what you think or if you have any questions by leaving a comment here.

Get the MATLAB code

Published with MATLAB® R2013b

## What Kind of MATLAB File is This?

We just had an interesting thread at MathWorks prompted by a customer request to programmatically distinguish between script and function MATLAB files. The File Browser already shows these separately. I will show some ideas we had and what we agreed was a good way to attack this.

### Contents

#### Function Vs. Script: The First Idea

The first idea was to create a function, perhaps called isfunction, and might be implemented by opening and reading the file, and looking for a function signature. Of course, you'd need to be sure that the function signature was not embedded in a comment. The tedium of doing this struck me and led me to ...

#### Function Vs. Script: The Next Idea

The next thing I thought of was to see what the function nargin did with various files on my path. I encourage you to read the documentation for this function. It's not very long.

Let's first see what happens for a "typical" function.

maxInputs = nargin('fftshift')

maxInputs =
2


We see that nargin finds that fftshift is called with up to 2 inputs.

Now let's try a function that can have any number of inputs.

maxInputs = nargin('ndgrid')

maxInputs =
-1


Here, the negative value indicates that you can input a variable number of inputs, and the first formal input is varargin.

If I were to try a script, e.g., this one, here's what I would type

    nargin('whatIsThis')

And here's the output from the command window.

  Error using nargin
whatIsThis is a script.
Error in whatIsThis (line 31)
nargin('whatIsThis')

We get an error message for this case. If we want to write a program to determine the file type, we can't have the program error out. So, in order to robustly check to see if the file is a function, we need to plan for scripts. For that, I recommend using the try/catch construct.

try
maxInputs = nargin('whatIsThis')
catch exception
if strcmp(exception.identifier, 'MATLAB:nargin:isScript')
disp('This file is a script.')
else
% We are only looking for scripts and functions so anything else
% will be reported as an error.
disp(exception.message)
end
end

This file is a script.

try
maxInputs = nargin('blogTopics.doc')
catch exception
if strcmp(exception.identifier, 'MATLAB:nargin:isScript')
disp('This file is a script.')
else
% We are only looking for scripts and functions so anything else
% will be reported as an error.
disp(exception.message)
end
end

Not a valid MATLAB file.


So now we can identify scripts without causing an error.

There's another kind of code file in MATLAB, pertaining to classes (the updated style, since R2008a, where classes are defined with (http://www.mathworks.com/help/matlab/ref/classdef.html classdef>), and we haven't done anything special to see if the file in question is for a class. I will demonstrate with the class dataset from Statistics Toolbox.

maxInputs = nargin('dataset')

maxInputs =
-1


The fact that dataset can take a variable number of inputs doesn't tell us that dataset is a class file. To probe for that information, we can use the function exist. Looking at the help, we see that if the file in question is a class file, then exist returns the value 8. To ask specifically if it's a class, we use this code.

classy = exist('dataset','class')

classy =
8


To be sure the files we probed originally are not class files, we need to check them out as well.

classy = exist('fftshift','class')

classy =
0

classy = exist('ndgrid','class')

classy =
0

classy = exist('whatIsThis','class')

classy =
0


#### Do You Need to Distinguish File Types in Your Work?

Do you have a similar need, where you need to characterize different file types, not based on their extension (e.g., .m)? I'd love to hear about cases where you need similar functionality. Let me know here.

Get the MATLAB code

Published with MATLAB® R2013a

## Zero Evolution

Today's post will take us on an historical tour of the function zeros, and pertains, in various ways, to the related functions ones, eye, false, true, rand, randn, and complex.

You might think "How can a function like zeros need to evolve at all?" - but you will see some of the MATLAB changes that lead to the evolution.

### Contents

#### Original Behavior of zeros

In Cleve's original Fortran MATLAB, there was no zeros function. However, the function ones had the following help entry:

  ONES  All ones. ONES(N) is an N by N matrix of ones.  ONES(M,N)
is an M by N matrix of ones.  ONES(A) is the same size as
A and all ones.

When I first started using MATLAB in 1987(version 3.05a), the behavior for the function zeros matched this ones design, with the obvious difference being the matrix had values of 0 instead of 1 for all of its entries.

Notice what this design means. I always get an array the same size as A unless A is a scalar. And then we get something with size A by A. This means that if we want to ever get a scalar value for zeros, we must constantly program defensively and ask if A is a scalar, and then create zeros(1), else create zeros(A). Constantly. This was a source of both confusion and bugs in code where people didn't recognize the discontinuity in behavior when they were programming.

#### Ambiguity Led to Change around MATLAB Version 4

To alleviate this nuisance, in MATLAB Version 4, we no longer used the syntax zeros(A) but instead expected the inputs to be non-negative integers. MATLAB still only had 2-dimensional double arrays (with an attribute describing when to interpret the values as characters). With this design, there was no more need for such defensive programming.

#### Extension to N Dimensions in MATLAB Version 5

Once we introduced higher dimensional arrays in MATLAB (version 5), we changed the input syntax for zeros to allow 2 scalar non-negative integer inputs OR a vector of such values for N-Dimensional arrays. While it was nice to be able to create higher dimensional arrays this way, it was tricky to create them for non-double datatypes, another new feature in that release of MATLAB.

The most common way to generate an array of zeros of another type was to first create the correctly sized double array and then convert it to, let's say uint8. This was inefficient, first creating an array using 8 times more memory than necessary, and then converting it to the uint8 array, with its smaller datatype.

Of course there were ways around this too. You could create a scalar zero of the right datatype, e.g., and then set that to be the final element in a newly created array, like this for your 4-D array:

     B(m,n,p,q) = uint8(zeros(1));

Effective, but perhaps a bit cryptic.

#### Ability to Create zeros of Particular Class in MATLAB R14 (Version 7)

In MATLAB 7, we introduced another new syntax for zeros and companions, allowing you to specify the datatype of the output, provided that it was one of the built-in datatypes. Here's a snippet of the relevant information.

   X = zeros(..., classname) classname restricted to built in types:
'double', 'single', 'int8', 'uint8', 'int16', 'uint16', 'int32',
'uint32', 'int64', or 'uint64'

This was definitely more memory efficient than creating an array of doubles and then coverting them. And less mysterious than setting the end element of a new array to a zero of the right type.

There was a remaining limitation however. It was still cumbersome to create an array of zeros of a user-defined datatype, including classes that shipped as part of our products.

#### Latest Syntax Updates in MATLAB R2013a (Version 8.1)

In addition to creating zero-arrays in all the ways possible since MATLAB Version 4, you can now create an array of zeros to be like another array, saving you the hassle of finding out that array's size and type and passing that information along. And it also allows you to create arrays of any datatype, built-in or user-defined. How do you do this? Here's one of the possible syntax choices.

    X = zeros(sz1,...,szN,'like',p)

#### Have You Been Able to Capitalize on These Changes?

Have you been able to take advantage of some of these changes over time? Did you know about them? I'd love to hear your thoughts here.

Get the MATLAB code

Published with MATLAB® R2013a

## Deconstructing Destructors

I'm pleased to introduce guest blogger Jennifer Black, manager of the MATLAB object system team. Today Jennifer will be sharing some thoughts on what happens when objects are destroyed in MATLAB and how you can control aspects of this process.

### Contents

#### Why Define a Destructor

When an instance of a class is cleared, MATLAB automatically cleans up the values of the properties owned by that object. For example, if I have a Signal class with a property that stores a matrix of double values, that matrix is destroyed when a Signal is destroyed. Sometimes, however, additional actions might be needed to handle external resources.

Let’s consider the following class:

type FileReader1

classdef FileReader1 < handle
properties
FileID
end
end



In MATLAB, a file identifier is an integer returned by the fopen function. It is used by file I/O functions such as fread and fwrite to access the file. If fopen cannot open the file, it returns -1.

If I have an open FileReader named reader in my workspace, clearing reader will cause MATLAB to clear the value stored in my FileID property. Now the number is gone, but the file is still open. Clearing FileID without first closing the file means the file is left open even though it is no longer being used. If this happens enough times I might run out of system file handles.

Let's look at defining a destructor to close the file handle. A destructor is a method of a class responsible for cleaning up resources owned by objects. In MATLAB, the handle superclass is used for all kinds of objects that have a unique identity independent of their current state. Unlike numbers, matrices, etc., handle objects represent unique things that have a beginning and an end and may change internal state along the way. Any subclass of handle may define a special method named delete, often referred to as a destructor. In the case of my FileReader class, I can correct the problem of leaking file handles by implementing my own destructor:

type FileReader2

classdef FileReader2 < handle
properties
FileID
end
methods
end
end
end



#### Defining a Destructor

A MATLAB destructor takes a single input argument - the object being destroyed - and returns no outputs. The input object is always scalar even if an array of MATLAB objects goes out of scope all at once.

In a MATLAB class you can define a method named delete that is not a destructor. A delete method with a different signature is not a destructor, nor is a static delete method. For example, defining the method to take more than one input argument or to return any output arguments means the method is not treated as an object destructor.

#### Calling Destructors

An object destructor can be called either explicitly by calling the delete method, or implicitly by MATLAB, for example when a variable is cleared or goes out of scope. Let's consider the difference between these two actions. To do so I will add a constructor and a readData method to my class:

type FileReader3

classdef FileReader3 < handle
properties(SetAccess = protected)
FileID = -1;
FileName
end

methods
if ischar(fname)
end
end

error('No file name has been specified for this FileReader.  No data will be read.');
else
colorData = colorData';
end
end

disp(s);
end
end
end
end



First, let's consider the case where I have one variable holding a FileReader:

myReader = FileReader3('colorData.txt');


Explicitly calling delete on myReader invokes the destructor and then destroys the FileReader. The myReader variable remains in my workspace, but its handle is no longer valid:

delete(myReader);

Closing colorData.txt

isvalid(myReader)

ans =
0


An implicit call to a destructor will occur if I clear myReader from my workspace:

myReader = FileReader3('colorData.txt');

Closing colorData.txt


The destructor will also be implicitly called if my variable goes out of scope, for example because the end of a function has been reached. To illustrate, let's add a helper function in which I create a FileReader:

type readDataFromFile

function colorData = readDataFromFile(filename)
end


When I call my new helper function, myReader is created in the function and will go out of scope when the function ends. This causes MATLAB to clear the variable, which results in an implicit call to the delete method:

colordata = readDataFromFile('colorData.txt');

Closing colorData.txt


Now, let's contrast what we have just discussed about having a single handle with the case where I have multiple handles to the same FileReader. This happens, for example, when I create a handle object and then assign that handle to another workspace variable:

reader1 = FileReader3('colorData.txt');


With these two handles to the same FileReader now in my workspace, I will explicitly call delete, and then use the isvalid method to see what happened to the FileReader object:

delete(reader1);

Closing colorData.txt

isvalid(reader1)

ans =
0

isvalid(reader2)

ans =
0


As expected, the reader1 handle is no longer valid, but note that reader2 is also no longer valid. Why did this happen? Because when I explicitly call delete, the destructor is called and the object is destroyed no matter how many handles there are referencing that object.

Now let's see what happens when a destructor is called implicitly. Once again let's create two handles to the same FileReader:

reader1 = FileReader3('colorData.txt');


But this time, rather than call the destructor explicitly, I will just clear one of the handles. As we saw earlier, this can result in an implicit call to the delete method:

clear reader1;

ans =
1


Why is reader2 still valid, and why was my destructor not called? Because MATLAB will only implicitly call a destructor when the last reference to an object is cleared. As reader2 still exists in my workspace, the underlying FileReader is not destroyed.

Most often a destructor is defined as a public method, but it can also be declared as a private or protected method. Making it private will prevent code outside the class from explicitly calling the destructor. Similarly, a protected destructor can only be explictly called from methods of the same class or from subclass methods. MATLAB will always be able to implicitly call a destructor, even one declared private or protected. You might choose to restrict access to your destructor in a situation where you want to prevent the accidental deletion of an object, such as when you have a singleton object.

#### Handles Contained within Other Structures

What happens when a handle is stored as a field of a struct or a cell of a cell array, or as a property of another object? When that top-level container variable is destroyed, the handle object will also be destroyed and its delete method implicitly called if and only if no other references to the object exist. Let's look at an example with structs:

s.myReader = FileReader3('colorData.txt');
clear s

Closing colorData.txt


Note that the FileReader stored in the myReader field of s has been destroyed. However, as we saw previously, another handle to the same FileReader will prevent the destruction of the object:

reader4 = FileReader3('colorData.txt');
clear s

ans =
1


Even after clearing s, we see that reader4 is still valid. The second handle prevented the implicit destruction of the object.

#### Classes in Hierarchies

So far we have been looking at examples of stand-alone classes, but what about classes that are part of a hierarchy? In MATLAB, every class has the ability to control its own destruction behavior. In a hierarchy of classes, every class can define its own destructor. As a result, a destructor cannot be defined as a Sealed method.

When an object is destroyed, MATLAB will call the destructor of the class of the object first, if it has been defined. Next, the destructor of each superclass is called, starting with the most direct superclass.

A destructor can reference the properties of the class itself, including properties inherited from superclasses, but it generally should not reference properties of a subclass. Why? Because when a superclass destructor is executing, the destructor of its subclass has already executed, and could have invalidated the values in its properties.

#### How Do You Use Destructors?

Now that we've discussed the basics of working with destructor methods in MATLAB, I'd like to hear of your experiences. Are you already using destructors? If you have any interesting applications or questions, I'd be very happy to hear about them here.

Get the MATLAB code

Published with MATLAB® R2013a

## Using memmapfile to Navigate through “Big Data” Binary Files

This week, Ken Atwell from MATLAB product management weighs in with using a memmapfile as a way to navigate through binary files of "big data".

memmapfile (for "memory-mapped file") is used to access binary files without needing to resort to low-level file I/O functions like fread. It includes an ability to declare the structure of your binary data, freely mixing data types and sizes. Originally targeted at easing the reading of lists of records, memmapfile also has application in big data. Today's post will examine column-wise access of big binary files, and how to navigate through metadata that sometimes is at the beginning of binary files.

### Contents

#### Experiment Parameters

To get started, create a potentially large 2D matrix that is stored on disk. numRows and numColumns can be changed to experiment with different sizes. To keep things simple and snappy here, the matrix is under a gigabyte in size. This is hardly "big data", and you can adjust the parameters here to create a larger problem. Do note that, of course, the disk space required to run this code will grow with the matrix size you create.

scratchFolder = tempdir;
numRows = 1e5;
numColumns = 1e3;


#### Create Test File

Create the scratch file. This can take from a moment to many minutes to run, depending on the sizes declared above. Because data of type double is being created, the file will consume 8*numRows*numColumns bytes of free disk space.

The value of [r,c] in the matrix is set to be c*1,000,000+r. This will make it easy to glance at our output and recognize that we are getting the values that are expected.

filename = ['mmf' int2str(numRows) 'x' int2str(numColumns) '.dat'];
filename = fullfile(scratchFolder, filename);
f = fopen(filename, 'w');
for colNum = 1:numColumns
column = (1:numRows)' + colNum*1000000;
fwrite(f,column,'double');
end
fclose(f);


#### memmapfile for Entire Data Set

To create a memory-mapped file, we call memmapfile with these two arguments:

1. The filename containing the data
2. The 'Format' of the data, which is a cell array with three components: a. The data type (double in this example), b. the size of the data (a matrix of size numRows by numColumns in this example), and c. a name to assign to this data (m for "matrix" in this example)

This is basic usage of memmapfile, and it encapsulates the entire data set in a single access. When working with "big data", you will want to avoid singular accesses like this. If the size of the data is large enough, your computer may become unresponsive (" thrash ") as it busily creates swap space in an effort to read in the entire matrix. The if statement is here to prevent you from doing this accidentally. If you are experimenting with data sizes larger than the physical memory available in your computer, you will want to skip this step.

% Prevent a memory-busting matrix from being created.
if numRows*numColumns*8 > 1e9
error('Size possibly too big; are you sure you want to do this?')
end

mm = memmapfile(filename, 'Format', {'double', [numRows numColumns], 'm'});
m = mm.Data.m;  %#ok<NASGU>


Regardless, clear m to free up whatever memory was used.

clear('m');


#### memmapfile with Columnwise Access

Here is a smarter way to access the big data a column at a time. Instead of creating a single variable that is numRows * numColumns large, we create a numRows * 1 vector, which is repeated numColumns times (note this code is now using the optional 'Repeat' argument to memmapfile). This subtle difference allows the big matrix to be read in one column at a time, presumably staying within available memory. The variable is named mj to indicate the 'j''th column of data.

mm = memmapfile(filename, 'Format', {'double', [numRows 1], 'mj'}, ...
'Repeat', numColumns);


The code spot-checks the 17th column.

if ~isequal(mm.Data(17).mj, (1:numRows)' + 17*1000000)
error('The data was not read back in correctly!');
end


memmapfile allows for creative uses of 'Repeat' if your application need it. For example, rather than a vector of an entire column, you can read in blocks of half a column:

memmapfile(filename, 'Format', {'double', [numRows/2 1], 'mj'}, 'Repeat', numColumns*2);


or blocks containing multiple columns:

memmapfile(filename, 'Format', {'double', [numRows*10 1], 'mj'}, 'Repeat', numColumns/10);


Of course, first ensure that your data's size is evenly divisible by these multiples, or you will create a memmapfile that does not accurately reflect the actual file that underlies it.

A note about memory-mapped files and virtual memory: If your application loops over many columns of memory-mapped data, you may find that memory usage as reported by the Windows Task Manager or the OS X Activity Monitor will begin to climb. This can be a little misleading. While memmapfile will consume sections of your computer's virtual memory space (only of practical consequence if you are still using a 32-bit version of MATLAB), physical memory (RAM) will not be used. The assignment of m above has the potential to fail only because that operation is pulling the contents of the entire memmapfile into a workspace variable, and workspace variables (including ans) reside in RAM. A comprehensive discussion of virtual memory is beyond the scope of this blog; the Wikipedia article on virtual memory is a starting point if you want to learn more.

#### Data File with XML Header

The above code assumes that the matrix appears at the very beginning of the data file. However, a number of data files begin with some form of metadata, followed by the "payload", the data itself.

For this blog, a file with some metadata followed by the "real" data will be created. The metadata is expressed using XML-style formatting. This particular format was created for this post, but it is representative of actual metadata. Typically, the metadata indicates an offset into the file where the actual data begins, which is expressed here in the headerLength attribute in the first line of the header. What follows next is a var to declare the name, type, and size of the variable contained in the file. This file will contain only one variable, but conceptually the file could contain multiple variables.

strNumC = int2str(numColumns);
strNumR = int2str(numRows);

'  <var name="mj" type="double" size="' strNumR ',' strNumC '"/>' char(10) ...
'</datFile>' char(10) ...
];


<datFile headerLength=00000095>
<var name="mj" type="double" size="100000,1000"/>
</datFile>


filename = ['mmf' int2str(numRows) 'x' int2str(numColumns) '_header.dat'];
filename = fullfile(scratchFolder, filename);
f = fopen(filename, 'w');
for colNum = 1:numColumns
column = (1:numRows)' + colNum*1000000;
fwrite(f, column, 'double');
end
fclose(f);


The header will now be read back in and parsed. While xlmread could be used to get a DOM node to traverse the XML data structure, regular expressions can often be used as a quick and dirty way to scrape information from XML. If you are unfamiliar with regular expressions, it is sufficient for this example just to understand that:

• (\d+) extracts a string of digits
• (\w+) extracts a word (an alphanumeric string)
• \s+ skips over whitespace

The first line of the file is read to determine the length of the header (extracted by a regular expression), and then the full header is read using this information. Finally, a second, more complex regular expression is used to extract the name, type, and size information for the variable contained in the binary data "blob" that follows the header.

f = fopen(filename, 'r');
firstLine = fgetl(f);
fclose(f);

firstLine %#ok<NOPTS>

firstLine =

% Get the length and convert the string to a double

headerLength =
95

f = fopen(filename, 'r');
fclose(f);

% Scan the metadata for type, size, and name
'tokens');


#### Create the Memory-mapped File

Lastly, create a memmapfile for the variable . The cell array returned by regexp is transformed into a new cell array that matches the expected input arguments to the memmapfile function.

% Reorganize the data from XML into the form expected by memmapfile
mmfFormater = {...
'Format', ...
{vars{1}{2}, ...
[str2double(vars{1}{3}), 1], ...
vars{1}{1}} ...
'Repeat', str2double(vars{1}{4})};

mm = memmapfile(filename, 'Offset', headerLength, mmfFormater{:});
mj = mm.Data(17).mj;  % Check the 17th column
if ~isequal(mj, (1:numRows)' + 17*1000000)
error('The matrix ''mj'' was not read in correctly!');
end


#### Conclusion

I hope this blog will be useful to those readers struggling to import big blocks of binary data into MATLAB. Though not covered in this post, memmapfile can also be used to load row-major data, and 2D "tiles" of data.

When you are done experimenting, remember to delete the scratch files you have been creating.

Have you used memmapfile or some other technique to incrementally read from large binary files? Share your tips here!

Get the MATLAB code

Published with MATLAB® R2013a

## Running Monte Carlo Simulations on Multiple GPUs

Today I'd like to introduce James Lebak. James is a developer who works on GPU support in the Parallel Computing Toolbox.

### Contents

#### Basic Option Pricing on the GPU

One common use of multi-GPU systems is to perform Monte Carlo simulations. The use of more GPUs can increase the number of samples that can be simulated, resulting in a more accurate simulation.

Let's begin by looking at the basic option pricing simulation from the Parallel Computing Toolbox example "Exotic Option Pricing on a GPU using a Monte-Carlo Method". In this example, we run many concurrent simulations of the evolution of the stock price. The mean and distribution of the simulation outputs gives us a sense of the ultimate value of the stock.

The function simulateStockPrice that describes the stock price is a discretization of a stochastic differential equation. The equation assumes that prices evolve according to a log-normal distribution related to the risk-free interest rate, the dividend yield (if any), and the volatility in the market.

type simulateStockPrice;

function finalStockPrice = simulateStockPrice(stockPrice, rate,...
dividend, volatility,...
Tfinal, dT)
% Discrete simulation of stock price
t = 0;
while t < Tfinal
t = t + dT;
dr = (rate - dividend - volatility*volatility/2)*dT;
perturbation = volatility*sqrt(dT)*randn();
stockPrice = stockPrice*exp(dr + perturbation);
end
finalStockPrice = stockPrice;
end



Achieving good performance on a GPU requires executing many operations in parallel. The GPU has thousands of individual computational units, and you get best performance from them when you execute hundreds of thousands or millions of parallel operations. This is similar to the advice we often give to vectorize functions in MATLAB for better performance. The function simulateStockPrice produces a single, scalar output, so if we run it as it is, it won't achieve good performance.

We can execute the simulation many times on the GPU using ARRAYFUN. The arrayfun function on the GPU takes a function handle and a gpuArray input, and executes the function on each element of the input gpuArray. This is an alternative to vectorizing the function, which would require changing more code. It returns one output for each element in the input. This is exactly what we want to do in a Monte Carlo simulation: we want to produce many independent outputs.

The runSimulationOnOneGPU function uses arrayfun to perform many simulations in parallel and returns the mean price.

type runSimulationOnOneGPU.m;

function mfp = runSimulationOnOneGPU(Nsamples)
% Run a single stock simulation on the GPU and return the
% mean file price on the CPU.

stockPrice   = 100;   % Stock price starts at $100. dividend = 0.01; % 1% annual dividend yield. riskFreeRate = 0.005; % 0.5 percent. timeToExpiry = 2; % Lifetime of the option in years. sampleRate = 1/250; % Assume 250 working days per year. volatility = 0.20; % 20% volatility. % Create the input data. Any scalar inputs are expanded to the size of the % array inputs. In this case, the starting stock prices is a vector whose % length is the number of simulations to perform on the GPU at once, and % all of the other inputs are scalars. startPrices = stockPrice*gpuArray.ones(Nsamples, 1); % Run all Nsamples simulations in parallel with one call to arrayfun. finalPrices = arrayfun( @simulateStockPrice, ... startPrices, riskFreeRate, dividend, volatility, ... timeToExpiry, sampleRate ); mfp = gather(mean(finalPrices)); end  We execute this function with a relatively small number of samples and display the computed mean. nSamples = 1e6; meanFinalPrice = runSimulationOnOneGPU(nSamples); disp(['Calculated mean final price of ', num2str(meanFinalPrice)]);  Calculated mean final price of 98.96  #### Using Multiple GPUs on One Machine Assume that we want to run this simulation on a desktop machine with two GPUs. This will allow us to have more samples in the same amount of time, which by the law of large numbers should give us a better approximation of the mean final stock price. With Parallel Computing Toolbox, MATLAB can perform calculations on multiple GPUs. In order to do so, we need to open a MATLAB pool with one worker for each GPU. One MATLAB worker is needed to communicate with each GPU. The workers in the pool will do all the calculations, so the client doesn't need to use a gpuDevice. Deselect the device on the client, so that all the memory on all the devices is completely available to the workers. gpuDevice([]);  Determine how many GPUs we have on this machine and open a local MATLAB pool of that size. This gives us one worker for each GPU on this machine. nGPUs = gpuDeviceCount(); matlabpool('local', nGPUs);  Starting matlabpool using the 'local' profile ... connected to 2 workers.  The runSimulationOnManyGPUs function contains a parfor loop that executes nIter times. Each iteration of the parfor loop executes nSamples independent Monte Carlo simulations on one GPU and returns the mean final price over all those simulations. The output of the overall simulation is the mean of the individual means, because each parfor iteration executes the same number of independent simulations. When the loop finishes we have executed nIter*nSamples simulations. type runSimulationOnManyGPUs.m;  function [tout, mfp] = runSimulationOnManyGPUs(nSamples, nIter) % Executes a total of nSamples*nIter simulations, distributing the % simulations over the workers in a matlabpool in groups of nSamples each. tic; parfor ix = 1:nIter meanFinalPrice(ix) = runSimulationOnOneGPU(nSamples); end mfp = mean(meanFinalPrice); tout = toc;  To run nSamples iterations on each GPU, we call runSimulationOnManyGPUs with nIter equal to the number of GPUs. [tout, meanFinalPrice] = runSimulationOnManyGPUs(nSamples, nGPUs); disp(['Performing simulations on the GPU took ', num2str(tout), ' s']); disp(['Calculated mean final price of ', num2str(meanFinalPrice)]);  Performing simulations on the GPU took 8.2226 s Calculated mean final price of 98.9855  It is important that the results are returned as regular MATLAB arrays rather than as gpuArrays. If the result is returned as a gpuArray then the client will share a GPU with one of the workers, which would unnecessarily use more memory on the device and transfer data over the PCI bus. #### Multi-GPU Execution Details The recommended setup for MATLAB when using multiple GPUs is, as we discussed above, to open one worker for each GPU. Let's dig in a little more to understand the details of how the workers interact with GPUs. When workers share a single machine with multiple GPUs, MATLAB automatically assigns each worker in a parallel pool to use a different GPU by default. You can see this by using the spmd command to examine the index of the device used by each worker. spmd gd = gpuDevice; idx = gd.Index; disp(['Using GPU ',num2str(idx)]); end  Lab 1: Using GPU 1 Lab 2: Using GPU 2  NVIDIA GPUs can be operated in one of four compute modes: default, exclusive thread, exclusive process, or prohibited. The compute mode is shown in the 'ComputeMode' field of the structure returned by gpuDevice. If a GPU is in prohibited mode, no worker will be assigned to use that GPU. If a GPU is in exclusive process or exclusive thread mode, only one MATLAB worker will attempt to access that GPU. It is possible for multiple workers to share the same GPU, if the GPU is in 'Default' compute mode. When this is done, the GPU driver serializes accesses to the device. You should be aware that there will be a performance penalty because of the serialization. There will also be less memory available on the GPU for each worker, because two or more workers will be sharing the same GPU memory space. In general it will be hard to achieve good performance when multiple workers share the same GPU. You can customize the way that MATLAB assigns workers to GPUs to suit your own needs by overriding the MATLAB function selectGPU. See the help for selectGPU for more details. MATLAB initializes each worker to use a different random number stream on the GPU by default. Doing so ensures that the values obtained by each worker are different. In this example, we have for simplicity elected to use the default random number generation stream on the GPU, the well-known MRG32K3A stream. While MRG32K3A is a highly regarded algorithm with good support for parallelism, there are other streams that you can select that may provide better performance. The documentation page Using GPUArray describes the options available for controlling random number generation on the GPU, and lists the different streams that you can select. #### Using Multiple GPUs in a Cluster With all of these details in place, it's time to consider how to put multiple GPUs to use when hundreds of millions or even billions of simulations are desired. Such a situation might arise, for example, when we want different values of key parameters such as volatility or the dividend in each simulation. Let's assume we want to run 800 million simulations, and that we have a cluster with 16 NVIDIA C2050 GPUs available. The problem is that a single GPU has a limited amount of memory. A single NVIDIA C2050 compute card has about 3 GB of memory and can compute on the order of$5\times10^7$points in runSimulationOnOneGPU without running out of memory. We can call runSimulationOnOneGPU 16 times to complete all the desired simulations. On a single machine, doing so takes around 10 minutes. We measured the time that it takes to run all$8\times 10^8\$ simulations with up to 16 GPUs. As you can see, we got nearly-perfect linear scaling in this experiment.

How many GPUs do you use in your application, and how do they interact? Let me know in the comments.

Get the MATLAB code

Published with MATLAB® R2013a

Loren Shure works on design of the MATLAB language at MathWorks. She writes here about once a week on MATLAB programming and related topics.

These postings are the author's and don't necessarily represent the opinions of MathWorks.