Computer Science 121

Scientific Computing Winter 2012 Chapter 9 Computer Science 121 … and (un)related issues

9.1 Environments and Scope • There are rules expressing how we link variable names to values • Script : names are linked in workspace • Function: names are linked only inside function • In general, the “place” in which a variable is linked is called the variable’s environment. • The rules that determine how the linking takes place are called scoping rules.

9.1 Environments and Scope • When a function is invoked, a new environment is created for its variables – both the parameters (input vars) and the internal variables (“local” vars) • This environment contains a copy of each input value, which is why input values cannot be modified. • This environment is destroyed when the function returns (exits). • Scripts do not create a new environment, which is why they mess with your workspace variables.

9.2 The Debugger • How do we find bugs in our programs? • Leave off semicolons so we see intermediate values • Use the debugger • Debugger has two basic features • Stop your program’s execution at a crucial “breakpoint” • Examine values of variables in the current environment

9.2 Debugger Commands (Basic) • dbstop : set a breakpoint >> dbstop primeFactors >> dbstop in FindLargest at 6 • dbstep : step to next line • dbcont : continue running, till next breakpoint • dbquit : bail out of debugger

9.2 Debugger Commands (Environmental) • dbstack : show the “stack” of function invocations that got us where we are. • dbup : go up one level in the stack • dbdown : go down one level

9.3 Shared Environments • Recall “two views” of computation: transformation from • Input → Output • State → State • Functions (so far) limit us to the input → output view

9.3 Shared Environments • What if we wanted to add state to a function? • E.g., a simple counter (“clicker”): >> keepcount ans = 1 >> keepcount ans = 2 >> keepcount ans = 3

9.3 Shared Environments • We can use the special “declaration” persistent to tell Matlab to maintain the previous value of a local variable: function res = counter persistent count if isempty(count) % first time around count = 0; end count = count + 1; res = count;

9.3 Shared Environments: The Evil global declaration • With persistent, only the function containing the variable can “see” the variable. • Instead of persistent, you can use global, which makes the variable visible to all functions that declare it this way, and makes things easier. • But this undermines the whole purpose of functions – i.e., variables are no longer hidden! • So don’t use global– Matlab may not even support it much longer.

9.4 Scoping of Functions • Functions themselves (as opposed to their parameters and local variables) are available everywhere • In workspace • To other functions • Therefore, if two functions have the same name (in two different folders), we get a conflict. • How to resolve this?

9.4 Scoping of Functions • Recall the path concept from Chapter 5: to find a file (.m function, data, etc.), Matlab first looks in the working directory, then successively in the list of directories in your path (which you can view/modify by using Set Path...) • Therefore, name conflicts can be resolved by using the version of the function first encountered in the path list. This is what Matlab does. • But it is a BAD IDEA to rely on this mechanism – i.e, DO NOT re-use function names if you can avoid that!

When / Why to Write a Function • Biggest motivation is always : Do not repeat code! • But there are other reasons • “Segregation” : If a variable is used only once in your code, maybe you should put that part of your code into a function, so the rest of the code doesn't see the variable: recall amortized smallestfactor function....

function res = smallestfactor(n) persistent primes res = n; if isempty(primes) primes = [2 3]; end if primes(end) < sqrt(abs(n)) % we need more primes for k = primes(end)+2:sqrt(abs(n)) if isprime(k) primes(end+1) = k; end end end for k = 1:length(primes) % a different k if rem(n, primes(k)) == 0 res = primes(k); break end end

function res = smallestfactor(n) persistent primes res = n; if isempty(primes) primes = [2 3]; end if primes(end) < sqrt(abs(n)) primes = more_primes(primes, n); end for k = 1:length(primes) if rem(n, primes(k)) == 0 res = primes(k); break end end

function newprimes = moreprimes(oldprimes, n) newprimes = oldprimes; for k = oldprimes(end)+2:sqrt(abs(n)) if isprime(k) newprimes(end+1) = k; end end

When / Why to Write a Function • Biggest motivation is always : Do not repeat code! • But there are other reasons • Segregation • Abstraction: • Maybe you will want to use a different algorithm in the future . Putting your current algorithm into a separate function makes it easier to find and replace your old algorithm quickly.

function newstate = caiter(oldstate) newstate=zeros(size(oldstate)); % commented-out : yucky! %for i = 2:(length(oldstate)-1) % newstate(i)=rule22(oldstate(i-1),oldstate(i),oldstate(i+1)); %end newstate(2:length(oldstate)-1) = rule22(oldstate(1:end-2), ... oldstate(2:end-1), ... oldstate(3:end)); newstate(1)=rule22(oldstate(end),oldstate(1),oldstate(2)); newstate(end)=rule22(oldstate(end-1),oldstate(end),oldstate(1));

function newstate = update_innerV1(oldstate) for i = 2:(length(oldstate)-1) newstate(i)=rule22(oldstate(i-1),oldstate(i),oldstate(i+1)); end function newstate = update_innerV2(oldstate) newstate(2:length(oldstate)-1) = rule22(oldstate(1:end-2), ... oldstate(2:end-1), ... oldstate(3:end)); function newstate = caiter(oldstate) newstate=zeros(size(oldstate)); newstate = update_innerV1(oldstate); newstate(1)=rule22(oldstate(end),oldstate(1),oldstate(2)); newstate(end)=rule22(oldstate(end-1),oldstate(end),oldstate(1));

When / Why to Write a Function • Biggest motivation is always : Do not repeat code! • But there are other reasons • Segregation • Abstraction: • Maybe you will want to use a different algorithm in the future . Putting your current algorithm into a separate function makes it easier to find and replace your old algorithm quickly. • Also supports division of labor: Bob writes primeFactors; Alice writes smallestfactor .

9.6 Warnings and Errors (Skip 9.5) • Sometimes we want to let the programmer do something “wrong” without causing the program to crash. • E.g., divide by zero: >> x = 1/0; Warning: Divide by zero x = Inf • When writing your own functions, the choice is up to you: function res = average(vec) res = sum(vec)/ length(vec); >> average([]) % ???

9.6 Warnings and Errors function res = averageV2(vec) if length(vec) < 1 error('Empty vector, Einstein!') end res = sum(vec) / length(vec); function res = averageV3(vec) if length(vec) < 1 warning('Empty vector, ain''t no thing') res = NaN; % not a number else res = sum(vec) / length(vec); end

Catching Errors • Instead of just quitting, we can have our program respond to an error by doing some kind of “repair” or “panic mode” • We try to do the operation, and then catch the error: % compute means of vectors in cell array c for i = 1:size(c, 1) try means(i) = averageV2(c{i}); catch warning('Empty vector; using NaN mean') means(i) = NaN; end end

9.7 Testing Functions • eXtreme Programming • “More eyes make shallower bugs” : pair programming (think about it for your projects!) • Write your tests before you write your solution (otherwise, what problem are you solving?) • Test-cases • Compare to known answers • Look at “edge cases”

9.7 Testing Functions • Compare to known answers : e.g., product of prime factors of N should equal N: >> num = 83723423424323; >> f = primeFactors(num) f = 31 109211 2479703 prod(f) == num ans = 1

9.7 Testing Functions • Edge cases: what you usually don't consider a typical input to the function >> primeFactors(2) ans = 2 >> primeFactors(1) ans = [] % why? >> primeFactors(0) ans = [] >> primeFactors(-10) ans = [] >> primeFactors(4.1) ans = 4.1 % ???

%primeFactors(n) - prime factors of n • function factors = primeFactors(n) • factors = []; • remaining = n; • while remaining > 1 % Condition for continuing • sf = smallestfactor(remaining); • factors(end+1) = sf; % Update accumulator • remaining = remaining/sf; • end

% smallestfactor(n) - smallest integer factor • function res = smallestfactor(n) • res = n; • for k=2:sqrt(abs(n)) • if rem(n,k) == 0 • res = k; • break • end • end

9.7 Testing Functions • What do edge cases tell us? • Situations (like primeFactors) where we don't get an error and probably should (0, negative, non-integer) • Situations where we do get an error and probably shouldn't: function res = findSock(placesToLook) for i = 1:length(placesToLook) if contains(placesToLook{i}, 'sock') res = placesToLook{i}; return end end

9.7 Testing Functions • Large projects will often involve a “test harness” that puts the code through its paces automatically: function res = testPrimeFactors(ncases) % test on random #'s, returning problem cases res = []; edgecases = [-1 0 1 2 4.1]; randomcases = floor(1000000000000*rand(1,ncases)); for k = [edgecases randomcases] try factors = primefactors(k); if prod(factors) ~= k res = [res k]; end catch res = [res k]; end end

9.7 Testing: Final Thoughts • Testing is annoying, but it's much cheaper than paying the cost of untested code later. • “If there is no way to check the output of your program, ... you have left the realm of scientific computation and entered that of mysticism, numerology, and the occult.” • Never “program around the problem” (like Captain Crunch).

Test Cases: amortized primeFactors revisited >> n = [1:50000]; >> p = n(find(isprime(n))); >> compare(p(end)*p(end-1)) Original took 0.028368 seconds Amortized took 8.804079 seconds >> compare(p(end)*p(end-1)) Original took 0.022651 seconds Amortized took 0.005601 seconds >> compare(p(end)*p(end-1)) Original took 0.025047 seconds Amortized took 0.005447 seconds >> compare(p(end)*p(end-1)) Original took 0.025387 seconds Amortized took 0.005321 seconds

9.8 Optional and Default Arguments • Recall plot function: >> plot(x, y, 'ro') >> plot(x, y, 'b-’) >> plot(x, y) % same as above • In general, Matlab functions (should) support a quoted option which defaults to some sensible value if unspecified • This requires use of special local variable nargin - automagically ignore missing arguments

function res = nicecos(angle, units) • % “nice” cosine supporting degrees and • % radians (default radians) • ifnargin== 1 % one input argument • units = 'radians'; • end • switch(lower(units)) % ??? • case {'radians', 'rad', 'r'} • res = cos(angle); • case {'degrees', 'deg', 'd'} • res = cos(pi.*angle./180); • otherwise • error('Dammit, Beavis!') • end

function res = nicecos(angle, units) • % “nice” cosine supporting degrees and • % radians (default radians); less redundant • if nargin== 1 % one input argument • units = 'radians'; • end • switch(lower(units)) % ??? • case {'radians', 'rad', 'r'} • a = angle; • case {'degrees', 'deg', 'd'} • a = pi.*angle./180 • otherwise • error('Dammit, Beavis!') • end • res = cos(a);

function res = nicetrig(fcn, angle, units) • % “nice” trig function supporting degrees and • % radians (default radians), multiple functions • if nargin< 3 % two input arguments • units = 'radians'; • end • switch(lower(units)) % ??? • case {'radians', 'rad', 'r'} • a = angle; • case {'degrees', 'deg', 'd'} • a = pi.*angle./180 • otherwise • error('Dammit, Beavis!') • end • res = feval(fcn, a);

9.8 Optional and Default Arguments: Named Arguments • axes function lets us adjust the axes in a figure: >> plot(x, y) >> axes('Position', [.5 .5 2 3], ... 'Units', 'inches') >> axes('Units', 'inches' , … 'Position', [.5 .5 2 3])

9.8 Optional and Default Arguments: Named Arguments • axes function lets us adjust the axes in a figure: >> plot(x, y) >> axes('Position', [.5 .5 2 3], ... 'Units', 'inches' ) • Requires use of special varargin variable: function axes( varargin) for i = 1:2:length(varargin) param = varargin{i}; value = varargin{i+1}; switch (param) case 'Position': % etc.

9.8 Optional and Default Arguments: Named Arguments • varargin can follow required arguments: function plot( x, varargin ) if isempty(varargin) y = x; x = 1:length(y); % etc.

Computer Science 121