Showing posts with label matlab. Show all posts
Showing posts with label matlab. Show all posts

Monday, 3 January 2011

You know you're not a real programming language when...

So, you want to do this:

myf = @(key, dir) if (ls('dir')) do_store(key, dir) end;

But you can't, as the syntax is only:

myf = @(args) expr


And expr can't be a conditional.

What you can do, to fake this, and display to the world that you are covered with infectious language lice, is hide the if statement inside an eval() call.

myf = @(key,dir) eval('if (ls(dir)); do_store(key, dir); end');

This, of course, only works for some use cases.

Now, most real languages have a ? : operator, so this isn't even an issue, but most languages with anonymous functions (lambdas, whatever) let you specify a whole code block, and even let you run it over multiple lines.

my $func = sub { ... };

If matlab let you define script local functions, this would also be less of a problem, as you could put the function body in the script file and reference its name, but you can't.  Nor can you access private/ functions from a script file.  Madness! And people wonder why scientists write disturbing code.  We really have no choice, sometimes.

Wednesday, 10 February 2010

Cthuloid horrors: Matlab parameter variation

When working out how much heat has flowed through the sea ice I've found a couple of good parameterisations (models) for the thermal conductivity of the ice. Now, I want to see if my choice of model seriously affects the results I'm getting, as I'm hoping that my conclusions will be robust to things I've haven't been able to measure directly myself.

To do this I could write several different scripts to run my heat flux calculations, or I can write one framework to handle the calculations then plug in different models and parameters from outside. Now, I'm using matlab, mostly because it makes very pretty graphs and makes doing maths very easy. It's not great, though, at making programming easy. Still, I've come up with the following fairly simple method for varying the parameters to a model while still making it easy to develop it with just one set to start with.

First we define our model a bit like this:

function result = hf_heat_flux(temperatures, salinity, args)
% defaults
d.RHO_SEA_ICE = 935;
d.CONDF = @hf_conductivity_mu;

% replace defaults from args
d = ag_struct_args(d, args);

% do calculation
...
conductivity = d.CONDF(d.RHO_SEA_ICE, temperatures, salinity);
...
results = make_up_numbers(conductivity, temperatures);
end

Now, ag_struct_args simply replaces structure members in d with those present in args, and returns the modified structure. This lets us have default parameters but then replace them from args easily.

function r = ag_struct_args( s, defaults )
%AG_STRUCT_ARGS fill in optional arguments
r = defaults;

n = fieldnames(s);
for k = 1:(size(n, 2))
name = n{1,k};
r.(name) = s.(name);
end

end

Now we can use our heat flux calculation for single runs easily. To do multiple runs we need another couple of tricks. First we make a function to iterate over an unknown number of varying inputs, which finally sends these inputs to our initial function in the args structure:

function r = ag_do_vary( fun, flatten, args, varyargs )
%AG_DO_VARY Run a function with varying inputs (eg. models)
if (~isempty(varyargs))
if (size(varyargs, 2) > 1)
newargs = {varyargs{1, 2:end}};
else
newargs = {};
end
fname = varyargs{1}.field;
for k = 1:size(varyargs{1}.values, 2)
args.(fname) = varyargs{1}.values{k};
r{k} = ag_do_vary(fun, flatten, args, newargs);
end
if (~isempty(flatten))
r = flatten(r');
end
return;
end;

% got the whole set of args sorted out
r.data = fun(args);
r.params = args;
end

A good flatten function is @cell2mat.

This returns an array of structures, each member having a .params entry containing the input to our model, and a .data entry letting us know the result. We call it like this, first producing an anonymous function to eventually call our model, and then calling it multiple times:

rho.field = 'RHO_SEA_ICE';
rho.values = {930 935 940};
cond.field = 'CONDF';
cond.values = {@cond_mu @cond_old};

fun = @(params) hf_heat_flux(T, S, params);
ret = ag_do_vary(fun, @cell2mat, struct(), {rho cond});