Hi:
I need to optimize the final output (a real number) of a long program
w.r.t. three parameters that determine how the program works. I have
prepared the smallest toy example I could think of below. My purpose is
to compare the program's prediction of the final number of individuals
in a cohort (as a result of the parameters that determine egg
development and egg mortality in the early life) with the same number
obtained from another sampling-model process, to make these two numbers
as close as possible.
I have found information to optimize a function (a mathematical
expression), for example with samin and bfgsmin, but I am not sure if
this is the right strategy to optimize the output of a long program. In
my program, there are a few places where random numbers are called and
used to create variety among the individuals.
The toy example is the minimum implementation of an individual-based
model; in this case I have 1000 super-individuals and each one starts
containing the same starting number of individuals but they differ
because they have a different potential of temperature-dependent growth
and wind-dependent mortality.
Shall I treat this problem as a function optimization problem, turning
my program into a function to be called by the optimizer?
Toy example follows. Real program takes about 3 minutes to finish on a
Windows Vista Core Duo at 1.73 Ghz and 1 GB RAM.
Thanks in advance
Ruben
%par.par, text file with five real numbers
par = load('par.par');
N0 = par(1); % start number of eggs, very large, ca. 10^12
b = par(2); % development parameter 1
c = par(3); % development parameter 2
d = par(4); % egg mortality parameter
e = par(5); % juvenile mortality parameter
f = par(6); % day when egg turn to juvenile
g = par(7); % final day of simulation
Nt = par(8); % number of individuals in final day, known from another
model-sampling process
%data, text file with two numeric columns and n rows, one column for
temperature, one for wind speed.
data = load('data.dat');
temp = data(:,1);
wind = data(:,2);
ind = 100000; number of super-individuals
winter = zeros(ind,6); %matrix containing the information about
super-individuals, overwritten at each time step.
%Column 1: Start number;
%Column 2: Temp-dependent development;
%Column 3: Surviving numbers;
for i=1:ind
winter(i,1) = a/ind; % starting number of individuals. Does not change
with time
winter(i,2) = random growth potential, uniform number between 0.5 and 1.5;
winter(i,3) = random susceptibility to wind-related mortality, normal
between 0.75 and 1.25;
winter(i,4) = development(a,b)*winter(i,2); turns to 1 when
development of egg into juvenile is completed
winter(i,5) = day when egg turns into juvenile; when winter(i,4) turns
to 1
winter(i,6) = number of individuals, at first increase due to eggs
transforming into juveniles (but controlled by wind-dependent egg
mortality), then decrease due to egg and juvenile tomortality
end
for day=1:g
winter(i,4) = degree of development dependent on daily temperature an
individual growth potential
winter(i,5) = written with index 'day' winter (i,4) turns to 1
winter(i,6) = gets updated with temperature-dependent transformation
of eggs into juveniles and wind-
dependent mortality of eggs, until day f, and later, with juvenile
mortality rate (independent of wind)
end
compute number of individuals in g, final day, Nt-program.
compute squared difference between Nt-program and Nt observed from
another model-sampling process
optimize squared difference w.r.t. parameters b,c,d