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});


No comments:

Post a Comment