{"id":1162,"date":"2015-04-01T08:19:12","date_gmt":"2015-04-01T13:19:12","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/?p=1162"},"modified":"2016-08-04T09:16:18","modified_gmt":"2016-08-04T14:16:18","slug":"the-winter-of-our-vectorization","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2015\/04\/01\/the-winter-of-our-vectorization\/","title":{"rendered":"The Winter of Our Vectorization"},"content":{"rendered":"<div class=\"content\"><!--introduction--><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#d34be643-051d-4777-8e23-0e0fa4fbc451\">Going the Whole 9 Feet<\/a><\/li><li><a href=\"#99044452-f5eb-418b-aea2-e866d2b0816c\">Enter MATLAB<\/a><\/li><li><a href=\"#6b5aad42-37e3-4d1d-a00c-87e93b3eb0e3\">Time for Some Logical Indexing<\/a><\/li><li><a href=\"#cd0cd71c-4737-4acd-bafe-7b8f1befc6b1\">Convoluted, Yet Simpler<\/a><\/li><li><a href=\"#61ecc103-058a-4fff-9e93-789e754adada\">Fun with <tt>arrayfun<\/tt><\/a><\/li><li><a href=\"#1d421e24-f8e8-450a-89c6-33a1931686af\">Putting it All Together<\/a><\/li><li><a href=\"#bbc77f33-8100-443c-a5c2-70d269bd865c\">Spring into Action!<\/a><\/li><\/ul><\/div><p>Today I'd like to introduce guest blogger Matt Tearle who works on our MATLAB Product Training materials here at MathWorks. Matt is on a mission to teach the world MATLAB, but this winter is testing his resolve. Annoyed that 22\" of snow forced him to reschedule a training, today he shows just how bad this winter had been.<\/p><h4>Going the Whole 9 Feet<a name=\"d34be643-051d-4777-8e23-0e0fa4fbc451\"><\/a><\/h4><p>It has been a brutal winter at <a href=\"https:\/\/goo.gl\/maps\/gKw0r\">MathWorks' headquarters in Natick, MA<\/a> (a little west of Boston). I recently read a fascinating article by MIT doctoral student Ben Letham who created a convincing visualization of just how relentless the snowfall has been this year. It was so convincing that I had to reproduce it myself (in MATLAB, of course!). In the process I had some fun with vectorization.<\/p><p>It's common to look at statistics like the largest amount of snowfall in a single day (23.5\" on Feb 17, 2003) or the highest total for the entire winter (108.5\" for 2014\/15). But 6\" of snow every day for a week would also cause chaos, even though it would seem mild by either of those measures. Hence, Mr. Letham's idea is to calculate the largest total of snow that fell during a range of given time periods -- three days, a week, two weeks, ... For example, 40.5\" of snow fell during the worst week-long period this winter, compared to just over 31\" in the worst week of 1995\/96 (which had the most snow of an entire winter until this year) and only(!) 28.5\" in 2002\/03 (23.5\" of which fell in one day).<\/p><p>When you plot these \"worsts\" for a range of time periods, you get:<\/p><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2015\/worstNDays.png\" alt=\"\"> <\/p><p>Sure enough, 2014\/15 had the worst week ever, the worst two weeks ever, the worst month ever, the worst of just about <b>any<\/b> time period ever! And by some distance, too. The winter of 1977\/78 is included because Boston locals talk with nostalgic reverence about the <a href=\"http:\/\/en.wikipedia.org\/wiki\/Northeastern_United_States_blizzard_of_1978\">blizzard of 1978<\/a> that brought the area to a standstill. 1978 was pretty bad, but nothing compared to the snowfall we've had this winter.<\/p><h4>Enter MATLAB<a name=\"99044452-f5eb-418b-aea2-e866d2b0816c\"><\/a><\/h4><p>Intrigued by Mr. Letham's analysis, I wanted to play with the numbers myself. Like any good researcher, Mr. Letham documented his sources, which turned out to be <a href=\"http:\/\/www.ncdc.noaa.gov\/\">NOAA's National Climatic Data Center<\/a>. This meant that I could <a title=\"http:\/\/www.ncdc.noaa.gov\/cdo-web\/search?datasetid=GHCND (link no longer works)\">retrieve the data for myself<\/a> and start playing (and so can you!). Before doing my own thing (which turned out to be fun, but a story for another day), I wanted to check that I could reproduce Mr. Letham's results. As it turned out, this helped me remember why I love MATLAB and gain appreciation for a function I had previously overlooked.<\/p><p>I downloaded the data from NOAA, used the <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/import_export\/import-data-interactively.html\">Import Tool<\/a> to import it into MATLAB, and did some basic preprocessing. As I often do, at this point I saved the cleaned-up data as a MAT-file, so I could now play with impugnity.<\/p><pre class=\"codeinput\">load <span class=\"string\">snowfalls<\/span>\r\n<\/pre><h4>Time for Some Logical Indexing<a name=\"6b5aad42-37e3-4d1d-a00c-87e93b3eb0e3\"><\/a><\/h4><p>To compare different winters, I need to be able to extract the snowfall data for the given time periods. This is easily achieved with a combination of the <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/matlab_prog\/represent-date-and-times-in-MATLAB.html\">date and time variable types that were introduced in R2014b<\/a> and my all-time favorite MATLAB construction: <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/math\/matrix-indexing.html#bq7egb6-1\">logical indexing<\/a>. I used the <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/import_export\/import-formatted-dates-and-times.html\">Import Tool to read the dates<\/a> as <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/datetime.html\"><tt>datetime<\/tt><\/a> variables. Now I needed to separate 2014\/15 from the rest of the data. Because the Northern Hemisphere winter crosses the new year, I decided to break years by the summer solstice (June 21):<\/p><pre class=\"codeinput\">pre2015 = (t &lt;= datetime(2014,6,21));\r\nsnowPre2015 = snow(pre2015);\r\nsnow2015 = snow(~pre2015);\r\n<\/pre><p>Numerical comparisons work naturally with <tt>datetime<\/tt> variables, so <tt>pre2015<\/tt> is a logical variable which I can use to extract the snowfall values. To extract the year of the legendary blizzard, I used the handy <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/isbetween.html\"><tt>isbetween<\/tt><\/a> function:<\/p><pre class=\"codeinput\">blizz = isbetween(t,datetime(1977,6,22),datetime(1978,6,21));\r\nsnow1978 = snow(blizz);\r\n<\/pre><h4>Convoluted, Yet Simpler<a name=\"cd0cd71c-4737-4acd-bafe-7b8f1befc6b1\"><\/a><\/h4><p>So now the crux of the problem: how to calculate the largest total over an $n$-day period. The typical programming language approach would be to loop over the array, taking $n$ elements starting with the $k^{th}$ element. But thinking about this process made me realize that I was looking for something very similar to a moving average. I knew that you can calculate moving averages in MATLAB by using convolution, implemented with the <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/conv.html\"><tt>conv<\/tt><\/a> function. In this case, I want the sum (rather than the sum divided by $n$, as I'd want for an average). Hence, I can calculate all the $n$-day totals in a vectorized way:<\/p><pre class=\"codeinput\">n = 7;\r\nndaytotals = conv(snow2015,ones(n,1));\r\nmax7day = max(ndaytotals)\r\n<\/pre><pre class=\"codeoutput\">max7day =\r\n   40.4724\r\n<\/pre><h4>Fun with <tt>arrayfun<\/tt><a name=\"61ecc103-058a-4fff-9e93-789e754adada\"><\/a><\/h4><p>Now I need to loop over $n$. I couldn't see any neat way to vectorize this, so maybe I should use a loop (remembering to preallocate):<\/p><pre class=\"codeinput\">nmax = 180;\r\nndaysmax = zeros(nmax,1);\r\n<span class=\"keyword\">for<\/span> n = 1:nmax\r\n    ndaysmax(n) = max(conv(snow2015,ones(n,1)));\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><p>I could also use <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/arrayfun.html\"><tt>arrayfun<\/tt><\/a> to perform this same operation without explicitly writing a <tt>for<\/tt>-loop:<\/p><pre class=\"codeinput\">nmax = 180;\r\ncalcndaymax = @(n) max(conv(snow2015,ones(n,1)));\r\nndaysmax = arrayfun(calcndaymax,1:nmax);\r\n<\/pre><p>But is this really achieving anything? The code is very slightly more compact, but arguably more complicated (at least to someone who isn't a hard-core MATLAB-phile). I admit that, for this very reason, I tend to look down a bit on <tt>arrayfun<\/tt> as \"how to cheat at vectorizing\". But in this case I realized I had actually stumbled upon a great reason to use <tt>arrayfun<\/tt>...<\/p><p>Having applied this clever convolution operation to the 2015 data, I now need to do the same to the pre-2015 data and the 1978 data (and, presumably, any others that might strike me as interesting). Should I copy and paste code? No! This sounds like a job for a function. The \"classical\" approach would be to write a function that took an array of data and a maximum value of $n$ as inputs. It would return the maximum total for each $n$ from 1 to $n_{max}$. This wouldn't be a particularly difficult function to write, given that the code is basically done already (above). But do I really need to write a function file just for this? No: <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/matlab_prog\/anonymous-functions.html\">anonymous functions<\/a> provide a way to create functions without writing a separate function file.<\/p><p>The one important caveat, though, is that anonymous functions must have a one-line definition (i.e., no loops!). So here's a compelling reason to use <tt>arrayfun<\/tt>: I can define my operation in a single line, enabling me to create an anonymous function, which can then be applied to any data set of interest:<\/p><pre class=\"codeinput\">nmax = 180;\r\nmaxXinNdays = @(x) arrayfun(@(n) max(conv(x,ones(n,1))),1:nmax);\r\n<\/pre><p>Now I can apply this function to the snowfall data:<\/p><pre class=\"codeinput\">ndaysmaxPre2015 = maxXinNdays(snowPre2015);\r\nndaysmax2015 = maxXinNdays(snow2015);\r\nndaysmax1978 = maxXinNdays(snow1978);\r\n<\/pre><h4>Putting it All Together<a name=\"1d421e24-f8e8-450a-89c6-33a1931686af\"><\/a><\/h4><p>Now I have a fairly short script that performs the whole analysis, leveraging the power of an anonymous function vectorized with <tt>arrayfun<\/tt>:<\/p><pre class=\"codeinput\">load <span class=\"string\">snowfalls<\/span>\r\n\r\n<span class=\"comment\">% Extract data for 2014\/15 winter and for everything else until then<\/span>\r\npre2015 = (t &lt;= datetime(2014,6,21));\r\nsnowPre2015 = snow(pre2015);\r\nsnow2015 = snow(~pre2015);\r\n<span class=\"comment\">% Extract data for the infamous winter of 1977\/78<\/span>\r\nblizz = isbetween(t,datetime(1977,6,22),datetime(1978,6,21));\r\nsnow1978 = snow(blizz);\r\n\r\n<span class=\"comment\">% Define a function that calculates the maximum total value in n<\/span>\r\n<span class=\"comment\">% consecutive elements of the array x, for n from 1 to nmax<\/span>\r\nnmax = 180;\r\nmaxXinNdays = @(x) arrayfun(@(n) max(conv(x,ones(n,1))),1:nmax);\r\n\r\n<span class=\"comment\">% Use the function to compare snowfall totals<\/span>\r\nfigure\r\nsemilogx(maxXinNdays(snowPre2015))\r\nhold <span class=\"string\">on<\/span>\r\nsemilogx(maxXinNdays(snow2015))\r\nsemilogx(maxXinNdays(snow1978))\r\nxlim([1 200])\r\n\r\n<span class=\"comment\">% Annotate the graph<\/span>\r\nxlabel(<span class=\"string\">'Number of days'<\/span>)\r\nylabel(<span class=\"string\">'Total snowfall (in)'<\/span>)\r\ntitle(<span class=\"string\">'Max total snowfall in an n-day period'<\/span>)\r\nlegend(<span class=\"string\">'pre-2015 record'<\/span>,<span class=\"string\">'2015'<\/span>,<span class=\"string\">'1978'<\/span>,<span class=\"string\">'Location'<\/span>,<span class=\"string\">'southeast'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2015\/winterOfOurVectorization_01.png\" alt=\"\"> <h4>Spring into Action!<a name=\"bbc77f33-8100-443c-a5c2-70d269bd865c\"><\/a><\/h4><p>You can <a href=\"https:\/\/github.com\/MattTearle\/snowfall\">download the data and code<\/a> and try this out for yourself. (The data also includes temperature readings.) What other statistics or metrics can you calculate with this technique of vectorizing anonymous functions with <tt>arrayfun<\/tt>? Was 2015 really that bad? Let us know what you discover about the winter of our discontent <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=1162#respond\">here<\/a>.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_843fc250ab9d4a019130d2835153d012() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='843fc250ab9d4a019130d2835153d012 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 843fc250ab9d4a019130d2835153d012';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2015 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_843fc250ab9d4a019130d2835153d012()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2015a<br><\/p><p class=\"footer\"><br>\r\n      Published with MATLAB&reg; R2015a<br><\/p><\/div><!--\r\n843fc250ab9d4a019130d2835153d012 ##### SOURCE BEGIN #####\r\n%% The Winter of Our Vectorization\r\n\r\n%%\r\n% Today I'd like to introduce guest blogger Matt Tearle who works on our\r\n% MATLAB Product Training materials here at MathWorks. Matt is on a mission\r\n% to teach the world MATLAB, but this winter is testing his resolve.\r\n% Annoyed that 22\" of snow forced him to reschedule a training, today he\r\n% shows just how bad this winter had been.\r\n\r\n%% Going the Whole 9 Feet\r\n% It has been a brutal winter at <https:\/\/goo.gl\/maps\/gKw0r MathWorks'\r\n% headquarters in Natick, MA> (a little west of Boston). I recently read a\r\n% <http:\/\/web.mit.edu\/bletham\/www\/winter2015.html fascinating article> by\r\n% MIT doctoral student Ben Letham who created a convincing visualization of\r\n% just how relentless the snowfall has been this year. It was so convincing\r\n% that I had to reproduce it myself (in MATLAB, of course!). In the process\r\n% I had some fun with vectorization.\r\n%\r\n% It's common to look at statistics like the largest amount of snowfall in\r\n% a single day (23.5\" on Feb 17, 2003) or the highest total for the entire\r\n% winter (108.5\" for 2014\/15). But 6\" of snow every day for a week would\r\n% also cause chaos, even though it would seem mild by either of those\r\n% measures. Hence, Mr. Letham's idea is to calculate the largest total of\r\n% snow that fell during a range of given time periods REPLACE_WITH_DASH_DASH three days, a\r\n% week, two weeks, ... For example, 40.5\" of snow fell during the worst\r\n% week-long period this winter, compared to just over 31\" in the worst week\r\n% of 1995\/96 (which had the most snow of an entire winter until this year)\r\n% and only(!) 28.5\" in 2002\/03 (23.5\" of which fell in one day).\r\n%\r\n% When you plot these \"worsts\" for a range of time periods, you get:\r\n%%\r\n% \r\n% <<worstNDays.png>>\r\n% \r\n%%\r\n% Sure enough, 2014\/15 had the worst week ever, the worst two weeks ever,\r\n% the worst month ever, the worst of just about *any* time period ever! And\r\n% by some distance, too. The winter of 1977\/78 is included because Boston\r\n% locals talk with nostalgic reverence about the\r\n% <http:\/\/en.wikipedia.org\/wiki\/Northeastern_United_States_blizzard_of_1978\r\n% blizzard of 1978> that brought the area to a standstill. 1978 was pretty\r\n% bad, but nothing compared to the snowfall we've had this winter.\r\n%\r\n%% Enter MATLAB\r\n% Intrigued by Mr. Letham's analysis, I wanted to play with the numbers\r\n% myself. Like any good researcher, Mr. Letham documented his sources,\r\n% which turned out to be <http:\/\/www.ncdc.noaa.gov\/ NOAA's National\r\n% Climatic Data Center>. This meant that I could\r\n% <http:\/\/www.ncdc.noaa.gov\/cdo-web\/search?datasetid=GHCND retrieve the\r\n% data for myself> and start playing (and so can you!). Before doing my own\r\n% thing (which turned out to be fun, but a story for another day), I wanted\r\n% to check that I could reproduce Mr. Letham's results. As it turned out,\r\n% this helped me remember why I love MATLAB and gain appreciation for a\r\n% function I had previously overlooked.\r\n%\r\n% I downloaded the data from NOAA, used the\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/import_export\/import-data-interactively.html\r\n% Import Tool> to import it into MATLAB, and did some basic preprocessing.\r\n% As I often do, at this point I saved the cleaned-up data as a MAT-file,\r\n% so I could now play with impugnity.\r\n\r\nload snowfalls\r\n\r\n%% Time for Some Logical Indexing\r\n% To compare different winters, I need to be able to extract the snowfall\r\n% data for the given time periods. This is easily achieved with a\r\n% combination of the\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/matlab_prog\/represent-date-and-times-in-MATLAB.html\r\n% date and time variable types that were introduced in R2014b> and my\r\n% all-time favorite MATLAB construction:\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/math\/matrix-indexing.html#bq7egb6-1\r\n% logical indexing>. I used the\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/import_export\/import-formatted-dates-and-times.html\r\n% Import Tool to read the dates> as\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/ref\/datetime.html |datetime|>\r\n% variables. Now I needed to separate 2014\/15 from the rest of the data.\r\n% Because the Northern Hemisphere winter crosses the new year, I decided to\r\n% break years by the summer solstice (June 21):\r\n\r\npre2015 = (t <= datetime(2014,6,21));\r\nsnowPre2015 = snow(pre2015);\r\nsnow2015 = snow(~pre2015);\r\n\r\n%%\r\n% Numerical comparisons work naturally with |datetime| variables, so\r\n% |pre2015| is a logical variable which I can use to extract the snowfall\r\n% values. To extract the year of the legendary blizzard, I used the handy\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/ref\/isbetween.html |isbetween|>\r\n% function:\r\n\r\nblizz = isbetween(t,datetime(1977,6,22),datetime(1978,6,21));\r\nsnow1978 = snow(blizz);\r\n\r\n%% Convoluted, Yet Simpler\r\n% So now the crux of the problem: how to calculate the largest total over\r\n% an $n$-day period. The typical programming language approach would be to\r\n% loop over the array, taking $n$ elements starting with the $k^{th}$\r\n% element. But thinking about this process made me realize that I was\r\n% looking for something very similar to a moving average. I knew that you\r\n% can calculate moving averages in MATLAB by using convolution, implemented\r\n% with the <https:\/\/www.mathworks.com\/help\/matlab\/ref\/conv.html |conv|>\r\n% function. In this case, I want the sum (rather than the sum divided by\r\n% $n$, as I'd want for an average). Hence, I can calculate all the $n$-day\r\n% totals in a vectorized way:\r\n\r\nn = 7;\r\nndaytotals = conv(snow2015,ones(n,1));\r\nmax7day = max(ndaytotals)\r\n\r\n%% Fun with |arrayfun|\r\n% Now I need to loop over $n$. I couldn't see any neat way to vectorize\r\n% this, so maybe I should use a loop (remembering to preallocate):\r\n\r\nnmax = 180;\r\nndaysmax = zeros(nmax,1);\r\nfor n = 1:nmax\r\n    ndaysmax(n) = max(conv(snow2015,ones(n,1)));\r\nend\r\n\r\n%%\r\n% I could also use <https:\/\/www.mathworks.com\/help\/matlab\/ref\/arrayfun.html\r\n% |arrayfun|> to perform this same operation without explicitly writing a\r\n% |for|-loop:\r\n\r\nnmax = 180;\r\ncalcndaymax = @(n) max(conv(snow2015,ones(n,1)));\r\nndaysmax = arrayfun(calcndaymax,1:nmax);\r\n\r\n%%\r\n% But is this really achieving anything? The code is very slightly more\r\n% compact, but arguably more complicated (at least to someone who isn't a\r\n% hard-core MATLAB-phile). I admit that, for this very reason, I tend to\r\n% look down a bit on |arrayfun| as \"how to cheat at vectorizing\". But in\r\n% this case I realized I had actually stumbled upon a great reason to use\r\n% |arrayfun|...\r\n%\r\n% Having applied this clever convolution operation to the 2015 data, I now\r\n% need to do the same to the pre-2015 data and the 1978 data (and,\r\n% presumably, any others that might strike me as interesting). Should I\r\n% copy and paste code? No! This sounds like a job for a function. The\r\n% \"classical\" approach would be to write a function that took an array of\r\n% data and a maximum value of $n$ as inputs. It would return the maximum\r\n% total for each $n$ from 1 to $n_{max}$. This wouldn't be a particularly\r\n% difficult function to write, given that the code is basically done\r\n% already (above). But do I really need to write a function file just for\r\n% this? No:\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/matlab_prog\/anonymous-functions.html\r\n% anonymous functions> provide a way to create functions without writing a\r\n% separate function file.\r\n%\r\n% The one important caveat, though, is that anonymous functions must have a\r\n% one-line definition (i.e., no loops!). So here's a compelling reason to\r\n% use |arrayfun|: I can define my operation in a single line, enabling me\r\n% to create an anonymous function, which can then be applied to any data\r\n% set of interest:\r\n\r\nnmax = 180;\r\nmaxXinNdays = @(x) arrayfun(@(n) max(conv(x,ones(n,1))),1:nmax);\r\n\r\n%%\r\n% Now I can apply this function to the snowfall data:\r\n\r\nndaysmaxPre2015 = maxXinNdays(snowPre2015);\r\nndaysmax2015 = maxXinNdays(snow2015);\r\nndaysmax1978 = maxXinNdays(snow1978);\r\n\r\n%% Putting it All Together\r\n% Now I have a fairly short script that performs the whole analysis,\r\n% leveraging the power of an anonymous function vectorized with |arrayfun|:\r\n\r\nload snowfalls\r\n\r\n% Extract data for 2014\/15 winter and for everything else until then\r\npre2015 = (t <= datetime(2014,6,21));\r\nsnowPre2015 = snow(pre2015);\r\nsnow2015 = snow(~pre2015);\r\n% Extract data for the infamous winter of 1977\/78\r\nblizz = isbetween(t,datetime(1977,6,22),datetime(1978,6,21));\r\nsnow1978 = snow(blizz);\r\n\r\n% Define a function that calculates the maximum total value in n\r\n% consecutive elements of the array x, for n from 1 to nmax\r\nnmax = 180;\r\nmaxXinNdays = @(x) arrayfun(@(n) max(conv(x,ones(n,1))),1:nmax);\r\n\r\n% Use the function to compare snowfall totals\r\nfigure\r\nsemilogx(maxXinNdays(snowPre2015))\r\nhold on\r\nsemilogx(maxXinNdays(snow2015))\r\nsemilogx(maxXinNdays(snow1978))\r\nxlim([1 200])\r\n\r\n% Annotate the graph\r\nxlabel('Number of days')\r\nylabel('Total snowfall (in)')\r\ntitle('Max total snowfall in an n-day period')\r\nlegend('pre-2015 record','2015','1978','Location','southeast')\r\n\r\n%% Spring into Action!\r\n% You can <https:\/\/github.com\/MattTearle\/snowfall download the data and\r\n% code> and try this out for yourself. (The data also includes temperature\r\n% readings.) What other statistics or metrics can you calculate with this\r\n% technique of vectorizing anonymous functions with |arrayfun|? Was 2015\r\n% really that bad? Let us know what you discover about the winter of our\r\n% discontent <https:\/\/blogs.mathworks.com\/loren\/?p=1162#respond here>.\r\n\r\n##### SOURCE END ##### 843fc250ab9d4a019130d2835153d012\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2015\/winterOfOurVectorization_01.png\" onError=\"this.style.display ='none';\" \/><\/div><p>ContentsGoing the Whole 9 FeetEnter MATLABTime for Some Logical IndexingConvoluted, Yet SimplerFun with arrayfunPutting it All TogetherSpring into Action!Today I'd like to introduce guest blogger... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2015\/04\/01\/the-winter-of-our-vectorization\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[62,57,33,3,12],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/1162"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=1162"}],"version-history":[{"count":4,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/1162\/revisions"}],"predecessor-version":[{"id":1968,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/1162\/revisions\/1968"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=1162"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=1162"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=1162"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}