Tucker McClure is an Application Engineer with The MathWorks. He spends his time helping our customers accelerate their work with the right tools and problem-solving techniques. Today, he'll be discussing how "functional programming" can help create brief and powerful MATLAB code.
Note that the recursive sequences we created in the last part could also have been implemented with for loops. For instance, here's factorial of n:
factorial = 1;
for k = 1:n
factorial = k * factorial;
end
Many times, recursive functions can be written iteratively in loops. However, we can't use for or while in an anonymous function, so instead of asking how we can unwrap recursive functions into iterative loops, let's ask the reverse: how can we implement loops with recursive functions?
Loops via Recursion
To loop properly, one must know:
What to do each iteration
If the process should continue to the next iteration
What's available when the loop begins
Allowing the "what to do" to be a function (fcn) of some state (x), the "if it should continue" to be another function (cont) of the state, and "what's available when the loop begins" to be the initial state (x0), we can write a loop function. This is a big step, so bear with me for some explanation!
On each step, the loop function will call the cont function, passing in all elements of the state, x, as in cont(x{:}). If that returns false (meaning we shouldn't continue), the current state, x, is returned. Otherwise, if we should continue, it calls fcn with all elements of the current state, as in fcn(x{:}), and passes the output from that to the next iteration. Letting this single iteration be denoted as f, we can build the anonymous function loop using our recur function.
For this trivial example, the state is simply the iteration count. We'll increase the count every iteration until the count >= n and return the final count. All this does therefore is count from 0 to the input n. Not very interesting, but it demonstrates the loop.
count = @(n) loop({0}, ... % Initialize state, k, to 0
@(k) k < n, ... % While k < n
@(k) {k + 1}); % k = k + 1 (returned as cell array)
arrayfun(count, 1:10)
ans =
[1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
I suppose that worked, but why are we using cell arrays to store the state, such as {0} and {k+1}? There are two reasons. First, if x is a cell array, then when we dump all elements of x into fcn, they become multiple arguments! That is, fcn(x{:}) is the same as fcn(x{1}, x{2}, ...). So instead of our function taking a big cell array for an input, it can take named arguments, which we'll use below. Second, we do this because it allows a function to return multiple elements that will be used by the next iteration, so if a function needed to return y and z, which would be arguments to the next iteration, it can simply return one cell array, {y, z}. It makes it easy to use. Here's a factorial example demonstrating this. The state is two different things: the iteration count, k, and factorial of the previous number, x. Note that both values of the state, k and x, are inputs to all of the functions. Note here how we're using @(k, x) for our functions. By allowing x to be a cell array, each element of the array becomes an argument such as k or x!
factorial = @(n) loop({1, 1}, ... % Start with k = 1 and x = 1
@(k, x) k <= n, ... % While k <= n
@(k, x) {k + 1, ... % k = k + 1;
k * x}); % x = k * x;
Call it:
factorial(5)
ans =
[6] [120]
Wait, we wanted 120 (the fifth number of the factorial sequence), so what's that 6 doing there?
A Better Loop
Remember how we return the full state? That's not very useful for this factorial example, as we get both k and the number we want as outputs in that cell array. Because the whole state isn't generally useful, let's add a cleanup function to our loop. We'll execute this when the loop is done (right after ~cont(...) returns false). Our cleanup function will take the full state and return only the important parts.
factorial = @(n) loop({1, 1}, ... % Start with k = 1 and x = 1
@(k,x) k <= n, ... % While k <= n
@(k,x) {k + 1, ... % k = k + 1;
k * x}, ... % x = k * x;
@(k,x) x); % End, returning x
The result:
factorial(5)
ans =
120
First seven numbers of factorial:
arrayfun(factorial, 1:7)
ans =
Columns 1 through 6
1 2 6 24 120 720
Column 7
5040
That's better.
I'll be the first to admit that the loop is a bit longer and much more rigid than a normal MATLAB loop. On the other hand, it can be used in anonymous functions, and its syntax has a certain cleanliness to it in that it doesn't modify any variables that live outside the loop; it has its own scope. This is one nice feature of loop being a function that takes code (functions) as arguments.
Doing More in a Loop
Let's say we want to do something else in the loop, but don't want its output passed to the next iteration, like printing something out. Remember the do_three_things example from last time? We executed numerous statements by putting them in a cell array and used curly to access the output we cared about. We can do that here, in a loop. For example, let's write out a function to print n digits of the factorial sequence. We'll use an array to store two things. The first will be the number that fprintf returns, which we don't care about. The second will be the cell array we want to return, k and x. We'll access that cell array with curly, as in curly({..., {k, x}}, 2), which just returns {k, x}.
say_it = @(k, x) fprintf('Factorial(%d): %d\n', k, x);
print_factorial = @(n) loop({1, 1}, ... % Start with k=1, x=1
@(k,x) k <= n, ... % While k <= n
@(k,x) curly({say_it(k,k*x),... % Print, discard
{k + 1, ... % k = k + 1;
k * x}}, ... % x = k * x;
2), ... % Return {k+1,k*x}.
@(k,x) x); % End, returning x
Now we're executing multiple things and only returning what we want while inside a loop built built on recursion and anonymous conditionals! We've come a long way since Part 1.
As a practical note, recall that because these loops use recursion, there's a limit to the number of times they can loop (MATLAB has a recursion limit, which is a setting in Preferences). Also, a recursive implementation of a loop isn't the most efficient. For this reason, it's best to implement loop itself in a file that can then be used in the same way. If it's in a file, it can also be kept on the MATLAB path so that it can be used anywhere.
function x = loop(x, cont, f, cleanup)
while cont(x{:})
x = f(x{:});
end
if nargin == 4
x = cleanup(x{:});
end
end
Final Example
This brings us to our final example. Below, we'll simulate a simple harmonic oscillator over time, using a structure to store dissimilar states, including a complete time history of the oscillator. This might simulate, for example, the sway of a lamp that's hanging from the ceiling after an earthquake.
% First, calculate a state transition matrix that represents a harmonic% oscillator with damping. Multiplying this by |x| produces |x| at a% slightly later time. The math here isn't important to the example.
Phi = expm(0.5*[0 1; -1 -0.2]);
% Now create the loop.
x = loop({[1; 0], 1}, ... % Initial state, x = [1; 0]
@(x,k) k <= 100, ... % While k <= 100
@(x,k) {[x, Phi * x(:, end)], ... % Update x
k + 1}, ... % Update k
@(x,k) x); % End, return x% Create a plot function.
plot_it = @(n, x, y, t) {subplot(2, 1, n), ... % Select subplot.
plot(x(n, :)), ... % Plot the data.
iif(nargin==4, @() title(t), ... % If there's a
true, []), ... % title, add it.
ylabel(y), ... % Label y
xlabel('Time (s)')}; % and x axes.% Plot the result.
plot_it(1, x, 'Position (m)', 'Harmonic Oscillator');
plot_it(2, x, 'Velocity (m/s)');
Summary
That's it for loops via recursion!
Let's look back at what we did over these three parts. First, we started with a simple map utility function to demonstrate the function-of-functions idea. Then we created our ubiquitous inline if, which further enabled recursion (a conditional is necessary to make recursion stop!). We also showed using multiple statements by storing their outputs in a cell array. Finally, we created a loop construct on top of our recursion functions.
At this point, we've done more than just scratch the surface of functional programming. We've used MATLAB's interesting constructs, such as function handles, cell arrays, and varargin to implement a functional programming framework, allowing a new syntax within MATLAB, where code can be arguments to flow control functions. Here's a roundup of what we created.
These have also been programmed as "normal" MATLAB functions so that they can be kept on the path and used whenever they're needed. These can be found under "Functional Programming Constructs" in File Exchange, here.
Thanks for reading. I hope this has both enabled a new level of detail in anonymous functions in MATLAB and helped demonstrate the wide range of possibilities available within the MATLAB language.
Do you have other functional programming patterns you use in your code? For instance, a do-while loop is just like our loop above except that it always runs at least one iteration. Any ideas how to program this or other interesting constructs in anonymous functions? Please let us know here!
Tucker McClure is an Application Engineer with The MathWorks. He spends his time helping our customers accelerate their work with the right tools and problem-solving techniques. Today, he'll be discussing how "functional programming" can help create brief and powerful MATLAB code.
Last time, we said that functional programming was marked by storing functions as variables (function handles) and working with functions that act on other functions. We put these ideas together to implement our own version of a map function for handling multiple inputs and outputs from multiple functions simultaneously, and we created iif, an "inline if", to allow the use of conditional statements inside of anonymous functions. So how might we work with recursive functions -- functions of themselves? We'll see how a functional programming style allows us to implement recursive functionality inside anonymous functions, and this will pave the way for the final part, in which we'll implement loops, without ever using for or while (which we can't use in anonymous functions).
Before we get started, let's implement iif again; we're going to need it frequently.
Recall that a recursive function is a function that calls itself. It therefore needs some way to refer to itself. When we write an anonymous function, it isn't "named" (hence, "anonymous"), so it can't call itself by name. How can we get around this?
Let's start with a Fibonacci sequence example. Recall that the nth number of the Fibonacci sequence is the sum of the previous two numbers, starting with 1 and 1, yielding 1, 1, 2, 3, 5, 8, 13, 21, etc. This is easy to implement recursively.
fib = @(n) iif(n <= 2, 1, ... % First two numbers
true, @() fib(n-1) + fib(n-2)); % All later numbers
But hey, that can't work! We haven't defined fib yet, so how could this anonymous function call it? In fact, the anonymous function will never "know" we're referring to it as fib, so this won't work at all. Therefore, instead of trying to call fib directly, let's provide another input: the handle of a function to call, f.
fib = @(f, n) iif(n <= 2, 1, ... % First two numbers
true, @() f(f, n-1) + f(f, n-2)); % All later numbers
Getting closer. Now, if we pass fibtofib along with the number we want, it will call fib, passing in fib as the first argument, recursively until we get our answer.
fib(fib, 6)
ans =
8
Ok, that's right. The sixth number of the sequence is 8. On the other hand, the syntax we've created is terrible. We have to provide the function to itself? I'd rather not. Instead, let's just write a new function that hands fib to fib along with the input n.
fib2 = @(n) fib(fib, n);
fib2(4)
fib2(5)
fib2(6)
ans =
3
ans =
5
ans =
8
That's a lot closer to what we want, but there's one more step. Let's write a function called recur to hand a function handle to itself, along with any other arguments. This makes recursion less cumbersome.
recur = @(f, varargin) f(f, varargin{:});
That was simple, so now let's re-write fib. The first argument to recur is the function, which we'll define inline. The second is n. That's all there is to it. It now reads as "Recursively call a function that, if k <= 2, returns one, and otherwise returns the recursive function of k-1 plus that of k-2, starting with the user's input n." (If it doesn't read quite this clearly at first, that's ok. It takes some getting used to. Comment liberally if necessary!)
That was a big step for anonymous functions, using both recursion and an inline conditional together with ease. Like map and iif, recur, looks strange at first, but once it's been seen, it's hard to forget how it works (just make one of the inputs a function handle and pass it to itself). And recursion doesn't have to stop at interesting mathematical sequences of numbers. For instance, in the next part, we'll use this to implement loops in, but first, we'll need a some helper functions and a good way to execute multiple statements in an anonymous function.
Helpers
These little functions are useful in many circumstances, and we're going to need curly frequently.
They allow us to write x(3, 4) as paren(x, 3, 4) and similarly for curly braces. That is, now we can think of parentheses and curly braces as functions! At first this might not seem useful. However, imagine writing a function to return the width and height of the screen. The data we need is available from this call:
get(0, 'ScreenSize')
ans =
1 1 1920 1200
However, we don't need those preceeding ones. We could save the output to a variable, say x, and then access x(3:4), but if we need this in an anonymous function, we can't save to a variable. How do we access just elements 3 and 4? There are numerous ways, but paren and curly are similar to constructs found in other languages and are easy to use, so we'll use those here.
Now we can write our screen_size function to return just the data we want.
With curly in place, let's examine something a little different. Consider the following:
do_three_things = @() {fprintf('This is the first thing.\n'), ...
fprintf('This is the second thing.\n'), ...
max(eig(magic(3)))};
do_three_things()
This is the first thing.
This is the second thing.
ans =
[25] [26] [15]
We've executed three statements on a single line. All of the outputs are stored in the cell array, so we have three elements in the cell array. The first two outputs are actually garbage as far as we're concerned (they're just the outputs from fprintf, which is the number of bytes written, which we don't care about at all). The last output is from max(eig(magic(3))); That is, the biggest eigenvalue of magic(3) is exactly 15. Let's say we just wanted that final value, the eigenvalue. It's the third element of the cell array, so we can grab it with curly.
do_three_things = @() curly({fprintf('This is the first thing.\n'), ...
fprintf('This is the second thing.\n'), ...
max(eig(magic(3)))}, 3);
do_three_things()
This is the first thing.
This is the second thing.
ans =
15
For a more complex example, let's say we want to write a function to:
Create a small figure in the middle of the screen
Plot some random points
Return the handles of the figure and the plot
Then by storing all of the outputs in a cell array and using curly to access the outputs we care about, we can make a multi-line function with multiple outputs, all in a simple anonymous function.
(As a quick aside, note that if a statement doesn't return anything, we can't put it in a cell array, and so we can't use it this way. There are ways around this, discussed here.)
To Be Continued
Today, we've come a long way, from a simple condition through recursion and executing multiple statements. Here's a roundup of the functions so far.
These can also be found here, implemented as regular MATLAB functions that can be kept on the path.
Next time, we'll look at loops. Until then, have you worked with functions such as paren or curly? How else are people implementing these or similar operations? Let us know here.
Tucker McClure is an Application Engineer with The MathWorks. He spends his time helping our customers accelerate their work with the right tools and problem-solving techniques. Today, he'll be discussing how "functional programming" can help create brief and powerful MATLAB code.
I use a lot of anonymous functions. They're nice and compact and almost invisible in their simplicity. Plus, if I can write an anonymous function to do something, I don't need to save the function in a file, and that can save me from file clutter on larger projects and from having to send someone a dozen files instead of sending one clean script. However, it seems at first glance like anonymous functions must necessarily be simple. No if... else, while, for, or any other keywords can be used. So how could we possibly write sophisticated programs in anonymous functions? We'll see, and it will involve some ideas from functional programming.
The goal of this introduction is to demonstrate how a few of these techniques can change the way we work in MATLAB, allowing greater brevity and simultaneously increasing readability. There are three parts. In this first part, we'll present creating functions of functions and treating functions as variables (in MATLAB, that means function handles), and from there, we'll move on to implementing conditional statements (like if... else) in anonymous functions. In the next part, we'll add recursion and executing multiple statements inside an anonymous function. In the final part, we'll develop a loop function. But first, if "function handle" or "anonymous function" is new to you, go check out Loren’s great introductions to those ideas in her previous posts!
Minimum and Maximum Example
Let's say we want to write a function to find the minimum and maximum of a set of numbers and store the results in an array. Here's a first pass:
Our min_and_max function takes in an array that we'll call x, finds the minimum and maximum, and stores the two results in an output array. Clear? Good. But now let's make it more difficult. The min and max functions both return two outputs if desired (both the minimum or maximum and the index at which they occur in the input array). Our simple min_and_max function can't get those secondary outputs! How can we access them? Consider this odd-looking line.
Well, that clearly worked. The minimum, 1, occurs at index 3. The maximum, 6, occurs at index 4, but what is this line actually doing? First, recall how cellfun behaves. The first argument is a function handle. The second argument is a cell array of whatever. Each element of the cell array is given as an argument to the provided function handle. Most of the time, that cell array is full of data, and each piece of data is passed to the function. However, we could just as easily put function handles in the cell array. Then the first function (@(f) f(...)) acts on all the other functions. So, first @min is passed in for f and the outputs from min([3 4 1 6 2]) are stored. Then, @max is passed in for f, and its outputs are stored.
Ok, now that we're working with functions of functions, let's remove that hard-coded [3 4 1 6 2] and write a new min_and_max function by simply adding a @(x) out front and changing [3 4 1 6 2] to x.
That might have looked a little funny, but it's pretty easy to think about, right? Now let's make it look a little nicer too.
Map
Above, we're mapping each function to our input x. More generally, we might write a "map" function to map a series of functions to the input values. We'll make val a cell array so we can also send multiple inputs to multiple functions all at once. This is like what we had before, but rearranged a bit.
map = @(val, fcns) cellfun(@(f) f(val{:}), fcns);
Look how simple this makes min_and_max (below), while still accessing both outputs. Not only is it shorter to write than any other versions so far, it's easier to read, with hardly anything but a single occurrence of each variable or function name. "Map x to the min and max functions". No problem.
What if outputs are different sizes? We'll write mapc (as in MAP with Cell array outputs) to handle this; all it needs is an extra argument to cellfun to say that our output isn't uniform in size.
Send pi to multiple functions that return differently-sized arrays. The first output is a scalar, the second is a scalar, and the third is a string.
mapc({pi}, {@(x) 2 * x, ... % Multiply by 2
@cos, ... % Find cosine
@(x) sprintf('x is %.5f...', x)}) % Return a string
ans =
[6.2832] [-1] 'x is 3.14159...'
That takes care of map, which we can now use anywhere to send a set of inputs to numerous functions and collect their multiple outputs with brief and easy-to-read code.
By the way, writing these functions that operate on other functions is part of the "functional programming" style, and we're just scratching the surface. Let's go a little deeper and see how we can write a function to choose which function to apply from a list of functions.
Inline Conditionals
Sometimes an anonymous function might need a condition, like if...else. However, normal MATLAB syntax doesn't allow program flow statements like these in anonymous functions. Hope it not lost. We can implement an "inline if" in a single line:
Alright, that looks decidedly strange, so before we discuss how it works, take a look at how easy it is to use:
[out1, out2, ...] = iif( if this, then run this, ...
else if this, then run this, ...
...
else, then run this );
All the "if this" conditions should evaluate to true or false. The "then run this" action next to the first true condition is executed. None of the other actions are executed! We could use this to make, for example, a safe normalization function to do the following:
If not all values of x are finite, throw an error.
Else if all values of x are equal to 0, return zeros.
Else, return x/norm(x).
This is implemented below. Note the @() out in front of the actions. This means, "don't do this action, but refer to this action". That is, we're passing pieces of code to the iif function as arguments. In this way, we aren't actually doing all three things; we'll only call the action for the single case we need.
try normalize([0 inf 2]), catch err, disp(err.message); end
Must be finite!
Test the all-zeros condition.
normalize([0 0 0])
ans =
0 0 0
Easy to use, right? We've implemented if... else behavior without needing an actual if or else anywhere! So now it's time to see how this thing works.
First, the iif function takes any number of arguments, thanks to varargin. These arguments will be condition 1 (true or false), action 1 (a function), condition 2, action 2, etc. First, the iif function selects all of the conditions (that's the odd numbered items in varargin) via [varargin{1:2:end}]. For our safe norm, this returns:
[~all(isfinite(x)), all(x == 0), true]
Next, it finds the index of the first true value in those conditions with find(..., 1, 'first'). E.g., if ~all(isfinite(x)) was false, but all(x == 0) was true, the index would be 2.
The actions to perform are the even-numbered items of varargin, so we just multiply that index by 2 to get the index of the action to perform. Finally, we execute the action by appending () on the end, as in
varargin{...}()
Did you catch what was happening there? We're passing little pieces of code as inputs to the iif function. Functions as arguments. See why this is called "functional" programming? I'll admit it looks weird at first, but once you've seen it, the pattern is hard to forget.
Ok, that's it for today. Here are the functions we developed today.
These can also be found here, implemented as regular MATLAB functions that can be kept on the path.
To Be Continued
In the next installment, we'll build on these to enable recursion and to make anonymous functions that execute multiple statements. In the mean time, there are multiple ways to accomplish an inline if. For instance, I've seen a very brief ternary operator in our forums in the past. I'm curious what other implementations people might suggest for this type of behavior, so please post suggestions here.
I just read the latest post on Steve Eddins image processing blog. I recommend you read it too. It discusses how to count objects in an image without introducing bias. It would be easy to count all the blobs you see, including partial ones cut by the border. But that would be an overestimate. So how do we get a better estimate?
The first technique Steve discusses is one I am familiar with already. The basic idea is to count all objects that don't touch any edges of the image, and only count objects if they touch two adjacent edges, say the north and west ones, equivalent to counting objects via the "lower right corner". That will, statistically speaking, give you an unbiased estimate of the count.
Weighted Count of Interior Objects
The second method instead does a weighted count of only the objects that touch no edges. The weighting accounts for the relative likelihood for an object of that size to encounter an edge from a particular view. See Steve's code and explanation to learn more of the details.
Today I welcome back guest blogger Ben Tordoff who previously wrote here on how to generate a fractal on a GPU. He is going to continue this GPU theme below, looking at how to measure the performance of a GPU.
Whether you are thinking about buying yourself a new beefy GPU or have just splashed-out on one, you may well be asking yourself how fast it is. In this article I will describe and attempt to measure some of the key performance characteristics of a GPU. This should give you some insight into the relative merits of using the GPU over the CPU and also some idea of how different GPUs compare to each other.
There are a vast array of benchmarks to choose from, so I have narrowed this down to three tests.:
How quickly can we send data to the GPU or read it back again?
How fast can the GPU kernel read and write data?
How fast can the GPU do computations?
After measuring each of these, I can compare my GPU with other GPUs.
How Timing Is Measured
In the following sections, each test is repeated many times to allow for other activities going on on my PC and the first-call overheads. I keep the minimum of the results, because external factors can only ever slow down execution.
To get accurate timing figures I use wait(gpu) to ensure the GPU has finished working before stopping the timer. You should not do this in normal code. For best performance you want to let the GPU carry on working whilst the CPU gets on with other things. MATLAB automatically takes care of any synchronisation that is required.
I have put the code into a function so that variables are scoped. This can make a big difference in terms of memory performance since MATLAB is better able to re-use arrays.
function gpu_benchmarking
gpu = gpuDevice();
fprintf('I have a %s GPU.\n', gpu.Name)
I have a Tesla C2075 GPU.
Test Host/GPU Bandwidth
The first test tries to measure how quickly data can be sent-to and read-from the GPU. Since the GPU is plugged into the PCI bus, this largely depends on how good your PCI bus is and how many other things are using it. However, there are also some overheads that are included in the measurements, particularly the function call overhead and the array allocation time. Since these are present in any "real world" use of the GPU it is reasonable to include these.
In the following tests, data is allocated/sent to the GPU using the gpuArray function and allocated/returned to host memory using gather. The arrays are created using uint8 so that each element is a single byte.
Note that PCI express v2, as used in this test, has a theoretical bandwidth of 0.5GB/s per lane. For the 16-lane slots (PCIe2 x16) used by NVIDIA's Tesla cards this gives a theoretical 8GB/s.
Peak send speed is 5.70217 GB/s
Peak gather speed is 3.99077 GB/s
On the plot, you can see where the peak was achieved in each case (circled). At small sizes, the bandwidth of the PCI bus is irrelevant since the overheads dominate. At larger sizes the PCI bus is the limiting factor and the curves flatten out. Since the PC and all of the GPUs I have use the same PCI v2, there is little merit in comparing different GPUs. PCI v3 hardware is starting to appear though, so maybe this will become more interesting in future.
hold off
semilogx(sizes, sendBandwidth, 'b.-', sizes, gatherBandwidth, 'r.-')
hold on
semilogx(sizes(maxSendIdx), maxSendBandwidth, 'bo-', 'MarkerSize', 10);
semilogx(sizes(maxGatherIdx), maxGatherBandwidth, 'ro-', 'MarkerSize', 10);
grid on
title('Data Transfer Bandwidth')
xlabel('Array size (bytes)')
ylabel('Transfer speed (GB/s)')
legend('Send','Gather','Location','NorthWest')
Test Memory-Intensive Operations
Many operations you might want to perform do very little computation with each element of an array and are therefore dominated by the time taken to fetch the data from memory or write it back. Functions such as ONES, ZEROS, NAN, TRUE only write their output, whereas functions like TRANSPOSE, TRIL/TRIU both read and write but do no computation. Even simple operators like PLUS, MINUS, MTIMES do so little computation per element that they are bound only by the memory access speed.
I can use a simple PLUS operation to measure how fast my machine can read and write memory. This involves reading each double precision number (i.e., 8 bytes per element of the input), adding one and then writing it out again (i.e., another 8 bytes per element).
sizeOfDouble = 8;
readWritesPerElement = 2;
memoryTimesGPU = inf(size(sizes));
for ii=1:numel(sizes)
numElements = sizes(ii)/sizeOfDouble;
data = gpuArray.zeros(numElements, 1, 'double');
for rr=1:repeats
timer = tic();
for jj=1:100
data = data + 1;
end
wait(gpu);
memoryTimesGPU(ii) = min(memoryTimesGPU(ii), toc(timer)/100);
endend
memoryBandwidth = readWritesPerElement*(sizes./memoryTimesGPU)/1e9;
[maxBWGPU, maxBWIdxGPU] = max(memoryBandwidth);
fprintf('Peak read/write speed on the GPU is %g GB/s\n',maxBWGPU)
Peak read/write speed on the GPU is 110.993 GB/s
To know whether this is fast or not, I compare it with the same code running on the CPU. Note, however, that the CPU has several levels of caching and some oddities like "read before write" that can make the results look a little odd. For my PC the theoretical bandwidth of main memory is 32GB/s, so anything above this is likely to be due to efficient caching.
memoryTimesHost = inf(size(sizes));
for ii=1:numel(sizes)
numElements = sizes(ii)/sizeOfDouble;
for rr=1:repeats
hostData = zeros(numElements,1);
timer = tic();
for jj=1:100
hostData = hostData + 1;
end
memoryTimesHost(ii) = min(memoryTimesHost(ii), toc(timer)/100);
endend
memoryBandwidthHost = 2*(sizes./memoryTimesHost)/1e9;
[maxBWHost, maxBWHostIdx] = max(memoryBandwidthHost);
fprintf('Peak write speed on the host is %g GB/s\n',maxBWHost)
% Plot CPU and GPU results.
hold off
semilogx(sizes, memoryBandwidth, 'b.-', ...
sizes, memoryBandwidthHost, 'k.-')
hold on
semilogx(sizes(maxBWIdxGPU), maxBWGPU, 'bo-', 'MarkerSize', 10);
semilogx(sizes(maxBWHostIdx), maxBWHost, 'ko-', 'MarkerSize', 10);
grid on
title('Read/Write Bandwidth')
xlabel('Array size (bytes)')
ylabel('Speed (GB/s)')
legend('Read+Write (GPU)','Read+Write (host)','Location','NorthWest')
Peak write speed on the host is 44.6868 GB/s
It is clear that GPUs can read and write their memory much faster than they can get data from the host. Therefore, when writing code you must minimize the number of host-GPU or GPU-host transfers. You must transfer the data to the GPU, then do as much with it as possible whilst on the GPU, and only bring it back to the host when you absolutely need to. Even better, create the data on the GPU to start with if you can.
Test Computation-Intensive Calculations
For operations where computation dominates, the memory speed is much less important. In this case you are probably more interested in how fast the computations are performed. A good test of computational performance is a matrix-matrix multiply. For multiplying two NxN matrices, the total number of floating-point calculations is
$FLOPS(N) = 2N^3 - N^2$
As above, I time this operation on both the host PC and the GPU to see their relative processing power:
sizes = power(2, 12:2:24);
N = sqrt(sizes);
mmTimesHost = inf(size(sizes));
mmTimesGPU = inf(size(sizes));
for ii=1:numel(sizes)
A = rand( N(ii), N(ii) );
B = rand( N(ii), N(ii) );
% First do it on the hostfor rr=1:repeats
timer = tic();
C = A*B; %#ok<NASGU>
mmTimesHost(ii) = min( mmTimesHost(ii), toc(timer));
end% Now on the GPU
A = gpuArray(A);
B = gpuArray(B);
for rr=1:repeats
timer = tic();
C = A*B; %#ok<NASGU>
wait(gpu);
mmTimesGPU(ii) = min( mmTimesGPU(ii), toc(timer));
endend
mmGFlopsHost = (2*N.^3 - N.^2)./mmTimesHost/1e9;
[maxGFlopsHost,maxGFlopsHostIdx] = max(mmGFlopsHost);
mmGFlopsGPU = (2*N.^3 - N.^2)./mmTimesGPU/1e9;
[maxGFlopsGPU,maxGFlopsGPUIdx] = max(mmGFlopsGPU);
fprintf('Peak calculation rate: %1.1f GFLOPS (host), %1.1f GFLOPS (GPU)\n', ...
maxGFlopsHost, maxGFlopsGPU)
hold off
semilogx(sizes, mmGFlopsGPU, 'b.-', sizes, mmGFlopsHost, 'k.-')
hold on
semilogx(sizes(maxGFlopsGPUIdx), maxGFlopsGPU, 'bo-', 'MarkerSize', 10);
semilogx(sizes(maxGFlopsHostIdx), maxGFlopsHost, 'ko-', 'MarkerSize', 10);
grid on
title('Matrix-multiply calculation rate')
xlabel('Matrix size (edge length)')
ylabel('Calculation Rate (GFLOPS)')
legend('GPU','Host','Location','NorthWest')
Comparing GPUs
After measuring both the memory bandwidth and calculation performance, I can now compare my GPU to others. Previously I ran these tests on a couple of different GPUs and stored the results in a data-file.
These tests reveal a few things about how GPUs behave:
Transfers from host memory to GPU memory and back are relatively slow, <6GB/s in my case.
A good GPU can read/write its memory much faster than the host PC can read/write its memory.
Given large enough data, GPUs can perform calculations much faster than the host PC, more than four times faster in my case.
Noticable in each test is that you need quite large arrays to fully saturate your GPU, whether limited by memory or by computation. You get the most from your GPU when working with millions of elements at once.
Today I would like to introduce guest blogger Jeremy Greenwald who works in the Development group here at MathWorks. Jeremy works on the Code Analyzer and will be discussing when preallocating MATLAB arrays is useful and when it should be avoided.
There are numerous resources that discuss preallocation, such as sections of our documentation and articles discussing improvements to MATLAB allocation strategies. While we will quickly review the topic of preallocation here, readers unfamiliar with this topic are encouraged to read some of the provided links.
Imagine we write the following small function to fetch our data from some external source. The function returns the variable data after assigning to it, one element at a time.
function data = fillData
for idx = 1:100
data(idx) = fetchData();
endend
MATLAB will reallocate memory numerous times while executing this loop. After reallocating memory, MATLAB has to copy the old values to the new memory location. This memory allocation and copying of values can be very expensive in terms of computation time. It also has the effect of increasing peak memory usage, since the old and new copy must both exist for a period of time.
In this example we know that the final size of the variable data is 1-by-100, so we can easily fix the issue by preallocating the variable with the zeros function. In this version of the function, there will only be a single memory allocation and the values of data never have to be copied from one location to another.
function data = fillDataWithPreallocation
data = zeros(1,100);
for idx = 1:100
data(idx) = fetchData();
endend
While this may not be an important optimization for small data sizes (such as 1-by-100), it can be a significant improvement if the size of the data is large. For example, in an image processing application, the data may consist of thousands of high resolution images, each image using hundreds of megabytes of memory. With such applications, correct usage of preallocation can lead to a significant improvement in execution time.
The Code Analyzer and the MATLAB Editor
The MATLAB Editor uses a feature called the Code Analyzer to detect certain programming patterns that may not be optimal. The Code Analyzer offers suggestions on how to rewrite these patterns. It then communicates with the Editor to underline such code. If you copy-and-paste the first function above into the MATLAB Editor, the variable data appears underlined in orange. Hovering over the variable with the cursor causes a tooltip to appear with the following message.
The variable 'data' appears to change size on every loop iteration.
Consider preallocating for speed.
The tooltip also contains a button labeled Details. Clicking on that button causes the tooltip box to expand and contain a fuller explanation of the message. Finally, inside the fuller explanation is a link to the section of the MATLAB documentation already mentioned in this post. MATLAB tries to offer a lot of guidance on when and how to preallocate. For the first function shown above
There are other code patterns that can also cause the size of a variable to change in such a way that preallocation would help. The Code Analyzer can catch many of these common patterns. The function below contains several examples.
function data = fillLotsOfData
% all three different variables are growing inside the loop% and all three are underlined in the MATLAB Editor
data2 = [];
data3 = [];
for idx = 1:100
data1(idx) = fetchData();
data2(end+1) = fetchSomeOtherData();
data3 = [ data3 fetchYetMoreData() ];
end
data = { data1, data2, data3 };
end
A Common Misunderstanding
Users have been told so often to preallocate that we sometimes see code where variables are preallocated even when it is unnecessary. This not only complicates code, but can actually cause the very issues that preallocation is meant to alleviate, i.e., runtime performance and peak memory usage. The unnecessary preallocation often looks something like this.
function data = fillDataWithUnecessaryPreallocation
% note the Code Analyzer message% The value assigned to variable 'data' might be unused.
data = zeros(1,100);
data = fetchAllData();
end
The variable data is first preallocated with the zeros function. Then it is reassigned with the return value of fetchAllData. That second assignment would not have caused the issue preallocation is meant to avoid. The memory allocated by the call to zeros cannot be reused for the data that is returned from fetchAllData. Instead, it is thrown away once the call to fetchAllData successfully returns. This has the effect of requiring twice as much memory as needed, one chunk for the preallocated zeros and one chunk for the return value of fetchAllData.
Note that if you copy-and-paste the above code into the MATLAB Editor, the following Code Analyzer message appears.
The value assigned to variable 'data' might be unused.
This is an indication that the values (and hence the underlying memory) first assigned to data will never be used. The appearance of this message on a line of code that is preallocating a variable is a good sign that the preallocation is unneeded. Since the Code Analyzer can detect numerous patterns that would benefit from preallocation, if the Code Analyzer does not detect such a pattern and it detects an unused variable, together these indicate a high likelihood that the preallocation is not needed. While the Code Analyzer may occasionally miss code patterns that could benefit from preallocation, it can be relied on to catch the most common such patterns.
Conclusions
Preallocating is not free. Therefore you should not preallocate all large variables by default. Instead, you should rely on the Code Analyzer to detect code that might benefit from preallocation. If a preallocation line causes the unused message to appear, try removing that line and seeing if the variable changing size message appears. If this message does not appear, then the original line likely had the opposite effect you were hoping for.
Did you see the variable unused message? Have you been confused by this message? What could the Code Analyzer have done to make it more clear that there was an issue? Let us know here.
Recently on the MATLAB newsgroup, there was a thread asking how to split up a vector into pieces which were each monotonically increasing by the value 1. The post got several answers which I did not read first. Here's my thinking.
I want to break this into 4 pieces: the runs 2:4 and 12:15, plus the scalars 7 and 9. Here's how I thought about the problem.
First I want to find the runs, which are elements that differ by 1. So I calculate the differences, making sure I include the final value in the array (which is why I append a 0 below).
ad = [diff(a) == 1 0]
ad =
1 1 0 0 0 1 1 1 0
Next, I figure out how many runs there are, by seeing how many 0 values are represented in the differences. This tells me how many arrays I will split my original array into.
numcells = sum(ad==0)
out = cell(1,numcells);
numcells =
4
Next I find the ending indices of the chunks by looking where ad is 0. Then I start to create contents for each cell in the cell array out, starting with index 1 in the original array.
indends = find(ad == 0);
ind = 1;
for k = 1:numcells
out{k} = a(ind:indends(k));
ind = indends(k)+1;
end
Here's what out looks like now.
out
out{:}
out =
[1x3 double] [7] [9] [1x4 double]
ans =
2 3 4
ans =
7
ans =
9
ans =
12 13 14 15
File Exchange Function
On the File Exchange, you can get a more general purpose function called SplitVec on SplitVec by Bruno. It does this and a whole lot more.
Do You Need to Partition Data?
If you need to partition data, I'd like to understand why, and how you do so. Let me know here.
Today I would like to introduce guest blogger John Kotwicki, who works in the Training Services group here at MathWorks. John manages the Certification Program and will be introducing the new MATLAB Certification offering.
I am excited to be Loren’s guest this week to talk about the new MathWorks Certification Program. Released earlier this year, the program focuses on validating skills within the core set of MATLAB competencies. There are currently two credentials available:
MathWorks Certified MATLAB Associate
MathWorks Certified MATLAB Professional
These new credentials target two distinct levels of proficiency. At the associate-level, candidates passing the exam possess the expertise to navigate the MATLAB environment, import data, call analysis functions, produce plots, and write scripts. In essence, they are confident using MATLAB for most daily tasks. Building on this foundation, the professional-level adds the ability to write applications in MATLAB, including graphical user interfaces. The emphasis at this level is the ability to produce efficient MATLAB code with a full understanding of the data types and function types available in the MATLAB language.
A big congratulations to the first group of individuals in both North America and Europe who have passed the MathWorks Certified MATLAB Associate exam! This achievement marks an important milestone in the certification track.
Why Certification?
For many years now, members of the MATLAB community have been asking for a way to demonstrate their proficiency to their employers or to their clients. Similarly, employers have been looking for ways to develop and validate skills for groups who rely heavily on MATLAB to get their work done, and to remain competitive in the marketplace. MathWorks certification aims to address these needs by providing a clear pathway to mastering the technology. When we designed the program we did not simply want to test memorization of syntax or knowledge of a few slick tricks. Instead the exam focuses on measuring a comprehensive ability to use MATLAB for the tasks valued most by employers. Candidates passing the exam are able to use MATLAB efficiently and produce quality code.
How do I become MathWorks Certified?
We have strived to make the process for participating in the program as easy as possible. To begin, start by reviewing the objectives tested for the Associate-level. If you have been using MATLAB for a number of years, you may feel confident in your abilities and ready to take the exam. However, if you would like a refresher, or you are just starting on your journey to MATLAB mastery, you may want to consider starting with a training course. Training courses are available that cover all of the objectives tested. (By the way, these courses are not designed specifically for the exam. They are also for those who need to get up-and-running with the tools and may not wish to pursue certification at this time.)
After taking a course, you will need to spend time practicing within MATLAB. Course materials provide many interesting examples and exercises that you can use to sharpen your skills. After you feel comfortable with the tools, try the practice questions to assess your readiness.
Candidates passing the exam will receive a certificate recognizing their achievement. They also receive a corresponding logo, which they can place on resumes, business cards, and Web sites to identify themselves as having earned the credential. Everyone who is certified will be enrolled in a MathWorks database with a unique ID, where employers can verify the credential and date. Last, but certainly not least, you will earn bragging rights. When it comes to MATLAB, you will be the expert!
How is the Credential Protected?
As is probably evident by now, creating a certification program is a lot of work. Of utmost concern to any certification program is safeguarding the value of the credential. I thought it would be interesting to discuss some of what we do to protect the integrity of the credential that candidates work hard to earn.
To begin, every exam is proctored in person. Candidates receive sealed copies of exams along with instructions for completing the exam, which are sealed once again at the end of the session. We do not work with any test prep companies, and have the ability to switch exams in the event we believe content has become comprised in any way.
The design of the program also features a strict separation from MathWorks service activities. No one working for MathWorks is eligible to take the exam. Not even the instructors teaching training courses know the actual questions asked on the exam, helping to guarantee everyone has a fair experience regardless of their interaction with MathWorks staff.
These practices represent just a sampling of the security policies in place to ensure fairness across all participants. Securing exams and establishing a fair experience for all participants demands constant attention.
What’s Next?
Now that several folks have taken and passed the Associate-level exam, we will be turning our attention to the Professional-level and looking for ways to expand access. We want to hear from you. What would you like to see in the program? If you have taken one of the exams, what has been your experience?
PS: For those of you who are signed up for any of the exams coming up in November: Good luck!
Today I’d like to introduce guest blogger Sarah Wait Zaranek who works for the MATLAB Marketing team here at MathWorks. Sarah previously has written about using GPUs in MATLAB. Sarah will be discussing how she got started using regular expressions.
Over the past few years, I have had the honor of doing several guest blog posts for Loren. Usually, I am blogging about something that I know quite well - but this time it is different. I wanted to write about regular expressions, talk a little about how I am starting to use them, and show some examples that I created along the way. My background is in computational geophysics, so I am pretty comfortable with numbers, parallel computing and a whole bunch of other MATLAB stuff. But, I never had to really manipulate strings. In my minor working with strings, I found that functions like strfind were enough for me to get the proverbial job done.
Well, then I found Cody. If you haven't started using Cody - you might not understand how addictive it can be! All those cool little coding puzzles in MATLAB, I just couldn't stop. However, I found myself pretty consistently skipping over any challenge that had to do with manipulating strings. I figured this was a sign that I had a hole in my MATLAB skills, and I needed to start remedying it. So, I started on a quest to learn more about regular expressions. The more I learn and play with them, the more I am impressed with just how powerful they are.
If you are new to regular expressions, I hope this blog post will inspire you to start embracing them as well. If you are an experienced regular expression user, hopefully you will enjoy some of my examples and find my newly budding excitement about regular expression amusing. You might want to check out this guest post made by one of our developers, Jason Breslau, which discusses the differences between the Perl and MATLAB implementations of regular expressions. Also, please consider posting your favorite examples in the comments at end of the post.
The Basics
Regular expressions are a way to describe a pattern within text. With regular expressions, you can match or alter parts (substrings) of a text string that match the described pattern. Regular expressions are found in text editors and in a range of languages including Perl, Java, Ruby, and of course, MATLAB.
In this post, I am going to focus on the function regexp. There are several other regular expression related functions in MATLAB, so I encourage you to read more about them as well.
In MATLAB, the calling syntax for regexp that we will be using is:
string is the string or the cell array of strings that I want to search for the pattern. expr is the regular expression that specifies the pattern I want to match. outselect specifies the output I want from the function, including such options as the location of the start or end of the substring that matches the expression, and the text of the substring of the input string that matches the pattern. All the possible output options are explained in more detail in the documentation.
Enough background, let's look at three examples of places where I have been using regular expressions lately.
Example #1 - Splitting a String into Separate Words
This example is so deceptively easy - but is so totally useful, I just had to include it. Have you ever wanted to split a string into separate words? This does it for you in one easy step.
Let's go through the basic syntax. First, I need to define the expression to specify the pattern I want to match in the string. In this case, I want to find the spaces - so I choose '\s' which represents any white space character. It is what is known as a "character type", and it represents either a specific set of characters or a certain type of character. The documentation has a list of possible character types available to use with regexp.
Then I have to decide what I want as outputs from regexp. In this case, I pick split, which indicates that I want to split the string into parts determined by the substring that matches the expression (i.e., break up the string into substrings based on where there is a space).
mystring = 'My name is Sarah Wait Zaranek';
splitstring = regexp(mystring,'\s','split');
disp(splitstring)
'My' 'name' 'is' 'Sarah' 'Wait' 'Zaranek'
The initial string is now broken up into a cell array containing all the separate words in the sentence.
I could do a similar thing if I wanted to break up a string based on sentences. In this case, I want to split at a space immediately proceeded by a !, ., or ?. This is slightly more complicated, and I might want to use a lookaround operator.
First, I have to figure out how to let regexp know that I want to match from a set of possible characters; I do this by enclosing the possible matches in square brackets. [!.?] means match any one of the character's listed.
Second, I indicate that I want to match a space preceded by any of those characters. To do so I use a lookaround operator. Lookaround operators let you look around a current position in a string. For instance, to look ahead (to the right) of a position to test if a particular expression is found, you use the lookahead operator (?=expr). In this case, I use a lookbehind operator (?<=expr) which allows me to look behind a current position to test if an expression is found. In particular, I am looking for matches where I find a space and when I "look behind" the space (to the left of it) I find a !, . and ?. I, again, can use the split output option to split by the matched substring.
mystring = 'My name is Sarah. I love MATLAB! Do you?';
splitstring = regexp(mystring,'(?<=[!.?])\s','split');
disp(splitstring)
'My name is Sarah.' 'I love MATLAB!' 'Do you?'
Example #2 - Creating Short Labels for a Plot
In this example, I was working on a problem that involved data from several locations in California. I wanted to plot data from several locations on the same plot and label each set of data accordingly. However, since the city names were long and there was a lot of data in my plot, I wanted to create abbreviated city names to do my labeling.
My first attempt was to do the command listed below. First, I want to find the locations of the capital letters in the cell array of strings. To do this, I can use [A-Z] character range operator which allows me to specify any character within the range of capital A to capital Z (aka any capital letter). The default output of regexp with a single output variable is to give me the position of the start of the match string, and I use that here. I can then use these locations to create my abbreviated city names by taking the capital letter and one letter to the right of it to create my abbreviation.
regexp returns a 1 x 6 cell array, each element holding the location of the capital letters for the corresponding input strings.
To extract the capital letters and the letters next to them, I use cellfun to operate on each element of the output indices and input string. I use sort to sort the indices into monotonically increasing order. This method assumes I have no capital letters in a row.
When I learned more about regular expressions, I discovered a new and cleaner way to accomplish the same task. I can extend the pattern to be any capital letter followed by any character. A dot (.) is used to represent any single character.
Since this actually matches the substrings I am interested in extracting, I don't need to output the indices. Instead, I just indicate I want the matched substrings by specifying 'match' as my output option. I, then, concatenate the matched substrings from each city name into a single abbreviation for that city by using cellfun.
Example #3 - Finding Data Marked by Repeat Letters
In this case, I am working with text strings that look something like - 'CC=0/CT=1/TT=5375'. I wanted to extract out the numeric values that follow a repeated letter. This is actually genetic data, and in this particular string I want to extract the number of people associated with either the CC or TT genotype.
There are a couple ways to approach this. Since there are only a few letters that could be present (A,G,C,T), I could use | as an alternative match operator. Unlike using [] as above in Example 1, | allows me to combine multiple expressions as the possible alternatives to match. If you enclose these alternatives with the square brackets, [], it will take each character including the | as the list of possible characters to match, so be sure to enclose with parentheses, ().
This difference is explained nicely in Friedl's Mastering Regular Expressions. Note: he uses character class to refer to using []. He states, "Don't confuse the alternation with a character class...A character class can match exactly one character, and that's true no matter how long or short the specified list of acceptable characters might be. Alternation, on the other hand, can have arbitrarily long alternatives, each textually unrelated to the other".
I first create an expression which represents the possible double letters followed by an equals sign, (CC|TT|AA|GG)=. Then, I place this expression in a lookbehind operator because I want to find the numbers immediately proceeded by this pattern, (?<=(CC|TT|AA|GG)=).
I use \d+ to define the number I want to match. This is made up of the metacharacter \d which represents any numeric digit and the + which is a quantifier. Quantifiers are used to match consecutive occurrences of a pattern with a string. In this case, + means match the pattern one or more times.
When these pieces are put together, They make an expression that represents a number that is preceded by a double CC, TT, AA, or GG and an equals sign, (?<=(CC|TT|AA|GG)=)\d+.
Alternatively, I could use tokens to find where letters were repeated. Although in this case since I have a relatively short list of possible repeated letters, tokens might be a bit heavy-handed. However, tokens can help me find any repeated letter and are useful if I don't want to write out all possible double letter combinations.
Parentheses allow you to group multiple characters and designate matched expressions found as tokens. Tokens allow you to remember matched elements and allow you to match other parts of the string with these captured elements. Using \N, I can reference the Nth matched token in my expression.
In this case, I want to find the letters in the string and see if the next letter matches the previous letter found. I do this by using \w to find a letter and grouping it in parentheses to make it a token. Then, I reference that token by a \1 to indicate we want to find instances where the matched letter was repeated. I follow that by an = and a \d+ as before. I capture the numeric value as output by exporting just the tokens and by indicating I want \d+ to be a token by enclosing it in parentheses. By choosing tokens as the output option, I can just get the matched tokens and not the whole matched string. I can then use cellfun to extract the 2nd token (the numeric values) from the output.
I hope you enjoyed this post on regular expressions. This is only the tip of the iceberg, and there is much more that regular expressions can do. Check out the documentation for more examples, and have fun!
If you are currently using regular expressions, do you have any advice for those new to using regular expressions? If you are new to using regular expressions, do you have any questions on getting started? Let me know by leaving a comment for this post.
Answer: When the sum of all divisors of a positive integer, except the number itself, equals the number. Here's a reference with more information from Wikipedia.
You may think that all we need to do is find all the prime factors and then we can test the sum of those versus the number itself. Searching MATLAB documentation for prime factors yields the function factor.
help factor
FACTOR Prime factors.
FACTOR(N) returns a vector containing the prime factors of N.
This function uses the simple sieve approach. It may require large
memory allocation if the number given is too big. Technically it is
possible to improve this algorithm, allocating less memory for most
cases and resulting in a faster execution time. However, it will still
have problems in the worst case, so we choose to impose an upper bound
on the input number and error out for n > 2^32.
Class support for input N:
float: double, single
See also PRIMES, ISPRIME.
Overloaded methods:
sym/factor
Reference page in Help browser
doc factor
Let's see if the number 6 is a perfect number (having divisors 1, 2, and 3).
n = 6;
facs = factor(n);
perfectOrNot = sum(facs) == n
perfectOrNot =
0
Really? But 6 is a perfect number! What's going wrong? Let's first look at facs.
facs
facs =
2 3
The number 1 is not a prime factor, but is a divisor and needs to be included in the calculation. We can either add it into the divisors or subtract it from n.
divisors = [1 factor(n)]
perfectOrNot = sum(divisors) == n
divisors =
1 2 3
perfectOrNot =
1
Try Some More Numbers
Make a function to do our comparisons now.
perfectOrNot = @(n) sum([1 factor(n)]) == n;
Check out some other numbers. First my favorite, 17.
n = 17;
perfectOrNot(n)
ans =
0
Next, another known perfect number, 28.
n = 28;
perfectOrNot(n)
ans =
0
What's Wrong Now?
Let's trace through the steps again.
divisors = [1 factor(n)]
sum28 = sum(divisors)
divisors =
1 2 2 7
sum28 =
12
Actually, these are not all the divisors of 28. Plus the number 2 is duplicated. We need to add all combinations of products of the prime factors, deduplicate the products, and omit the product of all of them (which is the number itself). For 28, that means including 2*2 and 2*7.
moreDivs = [2*2 2*7]
moreDivs =
4 14
Combine prime factors with combinations, introduce the 1, and deduplicate.
We next need an algorithm to calculate all the unique divisors of n. If you don't want to make up the algorithm yourself, you can check the File Exchange for entries that list divisors and not polynomials. You will find at least one candidate function there.
Or you could turn to the Symbolic Math Toolbox which includes some library functions for number theory - a likely place for us to look for some help. On the MuPAD side of the toolbox, we can take advantage of the MuPAD function numlib::divisors(). Since there is currently not a MATLAB version of this function available, I'm using a function from the MuPAD engine. To do that, I supply the library and name of the function, and the inputs.
Recent Comments