====== Matlab / Octave Utilities ======

This page list the various Matlab and octave utilities which interact with ADMB, typically generating input for and interpreting output from ADMB. It is meant to be as complete as possible. This means that it is likely to contain many non-essential files as well as the essential files, for defining and running data sets.

====== Files written spring 2010 ======

Here various .m-files which interact with ADMB will be listed.

FIXME **NOTE:** The generated .pin files are for the //simple// model, but can of course be modified manually. Adding zeros to the end of the .pin file until ADMB runs is usually sufficient.

==== example_real.m ====

This is an example of how to run and plot real data for the model with two Ms.
<file txt example_real.m>
% EXAMPLE script which runs the model with two Ms on real data,
% and then plots the results.

% If on Octave, turn off the pausing of text output
try,
  more off;
end;

% Enter the proper directory
%cd ~lennartfr/common_files/mfiles

% Define data files
catchall = '../ICES_AFWG_2009/ices_afwg_2009_table37_1_2_zero.txt';
surveyfiles = '../ICES_AFWG_2009/ices_afwg_2009_tableA13.txt';

% Define directories and file names
cod_dir = '../cod_two_Ms';
base_dir = pwd;
filename = [cod_dir '/cod'];

% Generate data from tables. This will create (or overwrite)
% the files cod.dat and cod.pin in the directory defined by cod_dir.
table2admb(catchall,surveyfiles,1985,1995,9,3,filename);

% Enter directory where admb program is found
cd(cod_dir);

% Note that table2admb generates pin file for the simple model, so we need to add one zero
% since we are running the model with two Ms
unix('echo 0 >> cod.pin');

% Run ADMB
unix('unset LD_LIBRARY_PATH; ./cod');

% Return to mfile directory and plot
cd(base_dir);
datpar2plot(filename);
</file>

==== example_synthetic.m ====

This is an example of how to run and plot synthetic data for the simple model.
<file txt example_synthetic.m>
% EXAMPLE script which runs the simple model on synthetic data,
% and then plots the results.

% If on Octave, turn off the pausing of text output
try,
  more off;
end;

% Enter the proper directory
% cd ~lennartfr/common_files/mfiles

% Define parameters for data generation
N0 = [200 400 600];
q = [0.4 0.6];
s = 0.1;
M = 0.15;
startyear = 1980;
num_survey_indices = inf;
num_years = 15;
catch_factor = 2;

% Define directories and file names
cod_dir = '../cod_admb';
base_dir = pwd;
filename = [cod_dir '/cod'];

% Generate data This will create (or overwrite)
% the files cod.dat and cod.pin in the directory defined by cod_dir.
% The return variable Nt is the true trajectory (transposed)
Nt = generate_data(N0,q,s,M,filename,length(N0),num_years, ...
     length(q),num_survey_indices,startyear,catch_factor);

% Enter directory where admb program is found
cd(cod_dir);

% Run ADMB
unix('unset LD_LIBRARY_PATH; ./cod');

% Return to mfile directory and plot
cd(base_dir);
datpar2plot(filename);
% Add true trajectory to plot 
add_N_to_plot(Nt,startyear,gcf)
</file>


==== add_N_to_plot.m ====
Adds the population trajectory defined in N_transpose to the current plot. Meant to add the true trajectory from generate_data.m to the plot made by datpar2plot.m, so you would do, for instance:
<code>
Nt = generate_data([800],[0.8 0.8],1e-1,0.15,'../cod_admb/cod',1,15,2,1*14*2,1980);
% (...run ADMB on the data set somehow, either from matlab/octave or in a separate terminal window...)
datpar2plot('../cod_admb/cod'); 
add_N_to_plot(Nt,1980,gcf)
</code>
<file txt add_N_to_plot.m>
function add_N_to_plot(N_transpose,first_year,figure_handle);
%  add_N_to_plot(N_transpose,first_year,figure_handle)
%
%  Given a figure made by datpar2plot, adds the population
%  trajectory defined in N_transpose. N_transpose has one
%  cohort per ROW, the oldest in the first row.
%
%  Cohorts are assumed to be born in consecutive years,
%  first_year being the first.
%
%  Written by Lennart Frimannslund, April 2010.

if (exist('figure_handle','var') == 1),
  % Set figure
  figure(figure_handle);

end;

N0scale = 1000;

%N = N_transpose';
colors = {'r' 'g' 'b' 'c' 'm' 'y' 'k'};
hold on;
for (i=1:size(N_transpose,2)),

  plot(first_year+(i-1)+[0:size(N_transpose,1)-1],N_transpose(:,i)/N0scale,sprintf('%s-^',colors{i}))

end;
</file>

==== datpar2plot.m ====
<file txt datpar2plot.m>
function datpar2plot(filename,new_figure);
% datpar2plot(filename,make_new_figure)
%
% Reads filename.par, filename.rep and filename.dat and creates a 
% plot of the population trajectory, with one color per cohort. In
% Addition, the survey data divided by q is plotted with the same color as the 
% corresponding cohort, and with the marker for the survey data
% being the same for data from the same survey (but potentially 
% from different cohorts).
%
% Written by Lennart Frimannslund, April 2010

% Get N matrix from filename.rep
extractpar;
% Scale back N to match scaling in cod.tpl
N0scale = 1000;
N = N/N0scale;

% Get stds from filename.std
%clear std_*;
std_filename = [filename '.std'];
%who('std_*');
extractstd;

%whos('std_*')
std_coeff = 1.96;

% Get first year from dat file. 
[filename '.dat']
[m,startyears,endyearsy,C,num_surveys,survey] = read_dat([filename '.dat']);
firstyear = startyears(1);

% Helper arrays for plotting
colors = {'r' 'g' 'b' 'c' 'm' 'y' 'k'};
styles = {'-' '--' '-.' ':'};
dotstyles = {'+' 'o' '*' 'x'};

width = size(N,2);

num_cohorts = size(N,1);

if (exist('new_figure') == 1 && new_figure == 0),
 % No new figure
 fprintf(1,'%s debug: No new figure requested.\n',mfilename);
else,
 figure;
end;
hold on;

% Plot trajectories
fprintf(1,'%s debug: num_cohorts = %i\n',mfilename,num_cohorts);

for (i=1:num_cohorts),
  plot(firstyear+(i-1):firstyear+(i-1) + (width-1) - (i-1),N(i,1: width - (i-1)),colors{i});
  try,
    %errorbar(firstyear+(i-1),N(i,1),std_coeff * std_N0(i,2),colors{i});
    %
    % ERRORBAR doesn't quite work in octave, so make our own bar
    plot([ firstyear+(i-1) ; firstyear+(i-1)],[N(i,1) - std_coeff * std_N0(i,2); N(i,1) + std_coeff * std_N0(i,2)], [colors{i} styles{1}]);

  catch,
    %fprintf(2,'Errorbar plot failed for cohort born in %4i -- wrong variable name?\n',firstyear+(i-1));
    %disp(lasterr);
    errorbar(firstyear+(i-1),N(i,1),std_coeff * std_N0(i,2));
  end;
end;


% Plot observations divided by q
startyears = firstyear:firstyear+num_cohorts-1;

for (i=1:size(survey,1)),
  
  % Check if q_alpha is present, which indicates sigmoid function
  if (exist('q_alpha','var')),
    q_sigmoid = q0 * ((exp(q_alpha*survey(i,2)+q_beta))/(1+exp(q_alpha*survey(i,2)+q_beta)));
    
    plot(startyears(survey(i,1))+survey(i,2), ...
        (1/N0scale)*survey(i,4)/q_sigmoid, ...
        [colors{survey(i,1)} dotstyles{survey(i,3)}]);
  elseif (exist('q02','var')),
    if (survey(i,2) >= 3),
       plot(startyears(survey(i,1))+survey(i,2),(1/N0scale)*survey(i,4)/q3plus(survey(i,3)),...
        [colors{survey(i,1)} dotstyles{survey(i,3)}]);
    else,
       plot(startyears(survey(i,1))+survey(i,2),(1/N0scale)*survey(i,4)/q02(survey(i,3)),...
        [colors{survey(i,1)} dotstyles{survey(i,3)}]);
    end;

  else,
    plot(startyears(survey(i,1))+survey(i,2),(1/N0scale)*survey(i,4)/q(survey(i,3)),[colors{survey(i,1)} dotstyles{survey(i,3)}]); 
  end;

end;

title(sprintf('Cod population trajectory. File: %s(.dat,.rep)',string_convert_underscore(filename)));
xlabel('Year');
ylabel('Volume in billions');

%curr_axis = axis;
%try,
%  axis([curr_axis(1) curr_axis(2) 0 2500/N0scale curr_axis(5) curr_axis(6)]);
%catch,
%  axis([curr_axis(1) curr_axis(2) 0 2500/N0scale]);
%end;
hold off;

% If q is a sigmoid function, plot this function
if (exist('q_alpha','var')),

sf = figure;
xx = 0:1:20;
yy = q0 * (exp(q_alpha*xx+q_beta)./(1+exp(q_alpha*xx+q_beta)));
plot(xx,yy);
title('Sigmoid function q0 * exp(q_a x + q_b)/(1+exp(q_a x + q_b))');
xlabel('Age');
ylabel('Catchability');

end;


return;

% Convert _ to \_ so that latex interpretation doesn't garble the string
function [out] = string_convert_underscore(input);

 num_underscores = length(find(input == '_'));
 if (num_underscores == 0),
   out = input;
 else,

   out = char(32*ones(length(input)+num_underscores,1));
   counter = 1;
   for (i=1:length(input)),
     if (input(i) == '_'),
       out(counter:counter+1) = '\_';
       counter = counter + 2;
     else,
       out(counter) = input(i);
       counter = counter + 1;
     end;
   end;
 end;

return;
</file>

==== extractpar.m ====
This file, despite its name, extracts the information in the .par **and** the .rep file to the current workspace.
<file txt extractpar.m>
% extractpar - Extraction of data in .par and .rep file to Matlab/Octave environment
%
% Input: filename, a string variable in the current workspace
%        with the path to the .par and .rep files without the
%        extensions, e.g. use filename = '../../cod' to indicate
%        ../../cod.par and ../../cod.rep.
%
% Output: N, matrix with population trajectory, read from .rep file.
%
%         The variables in the .par file, with the same names as in the 
%         .par file.
%
% NOTE: The .par file must be EXACTLY on the form provided by ADMB.
%       The .rep file must contain only numbers (defining N) first.
%
% Written by Lennart Frimannslund, April 2010.
%

if (exist('filename','var') == 0),

  fprintf(1,'%s: No filename specified!\n',mfilename);
  return;

end;


% Get the .par file
fid = fopen([filename '.par'],'r');

% Discard first line, which contains information about the objective function, etc.
l = fgetl(fid);
l = 0;

% Process remaining lines, in pairs. The first line of the pair has the name, the second the value.
while (l ~= -1),

  l = fgetl(fid);
  if (l == -1),
    break;
  end;
  varname = strip_hash_and_colon(l);
  l = fgetl(fid); 

  % Send the variable name and values to the command line
  eval(sprintf('%s = [ %s ];',varname,l));
  %sprintf('%s = [ %s ];',varname,l)
end;


fclose(fid);

% Get the .rep file

% Number of lines in the rep file
n_lines = length(N0);

N_cell = cell(n_lines,1);
max_width_so_far = 0;

fid = fopen([filename '.rep']);

% Get one line at a time
for (i=1:n_lines),

  output = fgetl(fid);
  N_cell{i} = str2num(output);  
  max_width_so_far = max(max_width_so_far,length(N_cell{i}));

end;

fclose(fid);

N = zeros(n_lines,max_width_so_far);

for (i=1:n_lines),
  N(i,1:length(N_cell{i})) = N_cell{i}(:)';
end;


% ...and clean up our temporary variables
clear N_cell output max_width_so_far n_lines i

return;
</file>

==== extractstd.m ====
<file txt extractstd.m>
% extractstd - script which extracts variable values and standard deviations from
%              the file specified in the variable std_filename.
% 
% Output:      Variables given in the std file, prepended with std_, e.g. std_N0.
%
% NOTE:        In the std file, vectors are simply given as several scalar variables
%              with the same name. In this case, these scalars are simply appended
%              to form a vector, if a variable name already exists. For this reason
%              the std-variables should be cleared before calling this function!!!
%
% Written by Lennart Frimannslund, April 2010

% Check for already existing variables starting with std_ (except std_filename, which should be present).
who_cell = who('std_*');
if (length(who_cell)>1),
  fprintf(2,'WARNING!!! Variables starting with %sstd_%s detected. %s might malfunction.\n',char(39),char(39),mfilename);
  fprintf(2,'Consider calling %sclear std_*%s\n',char(39),char(39));
end;

if (exist('std_filename','var') == 0),
  fprintf(1,'%s: Filename variable %sstd_filename%s not defined!\n',mfilename,char(39),char(39));
  return;
end;

fid = fopen(std_filename,'r');
if (fid == -1),
  fprintf(1,'%s: Unable to open file %s.\n',mfilename,filename);
  return;
end;

% Get rid of first line
l = fgetl(fid);

% Now get all variable names, values and stddevs, in columns 2,3 and 4, respectively.
while (l ~= -1),

  l = fgetl(fid);
  if (l==-1),
    break;
  end;

  [token, remainder] = strtok(l);
  [currname, remainder] = strtok(remainder);
  [currval, remainder] = strtok(remainder);
  [currstd, remainder] = strtok(remainder);

  if (exist(['std_' currname],'var') == 0),
    eval(sprintf('std_%s = [ %s %s ];',currname,currval,currstd));
  else,
    eval(sprintf('std_%s = [ std_%s; %s %s ];',currname,currname,currval,currstd));
  end;

end;


fclose(fid);
clear fid token l currname currval currstd remainder who_cell;
</file>

==== generate_data.m ====
This file is used for generating synthetic data. 

:!: **NOTE: Generates .pin file for simple model.**

FIXME **NOTE: The line**
<code>
  CC(i-1,j) = rand(1,1) * cfactor/n * N(i-1,j);
</code>
**...should probably be changed to be independent of n and have less variance. **
<file txt generate_data.m>
function N = generate_data(N0,q,s,M, filename, cohorts, years, surveys, requested_observations,year_offset,in_cfactor);
% filegenerator(N0,q,s,M, filename, cohorts, years, surveys, observations, year_offset, cfactor)
%
% Generates .dat file for cod problem. 
%
% VERSION WITH SIMPLE MODEL - no random effects, and just one M!
%
% Written by Lennart Frimannslund
 
if (exist('year_offset','var')==0 | length(year_offset) == 0),
  year_offset = 1980;
end;
 
% Generate filename if not specified
if (exist('filename','var') == 0 | length(filename) == 0),
  c = clock;
  filename = sprintf('herring-%s-%02d.%02d.%02d',date,c(4),c(5), ...
                     floor(c(6)));
  clear c;
end;
 
 
% Get number of surveys
if (exist('q','var') ~= 0 & length(q) ~= 0),
  surveys = length(q);
end;
 
% s - Generate random number between 0 and 1
if (exist('s','var') == 0 | length(s) == 0),
  s = rand;
end;
 
% M - Generate random number between 0 and 0.2
if (exist('M','var') == 0 | length(M) == 0),
  M = rand(1,1) * 0.2;
end;
 
 
% If no of cohorts not specified - generate random number between 1 and 6
if (exist('cohorts','var') == 0 | length(cohorts) == 0),
  cohorts = ceil((rand*5)) + 1;
end;
 
% years - if not specified, generate random number between 5 and 25
if (exist('years','var') == 0 | length(years) == 0),
  years = ceil((rand*20)) + 5;
end;
 
% surveys - if not specified, generate random number between 1 and 6
if (exist('surveys','var') == 0 | length(surveys) == 0),
  surveys = ceil((rand*5)) + 1;
end;
 
% if not specified - generate which fraction of maximum possible
% observations available
%if (exist('observations','var') == 0 | length(observations) == 0),
%  observations = ceil(rand * surveys * years * cohorts);
%end;
 
 
% Define N0 if not supplied 
% exp(N0) between 0 and 3
%if (exist('N0','var') == 0 | length(N0) == 0),
%  N0 = 3 * rand(1,1);
%end;
 
% Generate q - between 0.1 and 0.8 
if (exist('q','var') == 0 | length(q) == 0),
  q = 0.7 * rand(1,surveys) + 0.1;
end;
 
 
 
% Input variables now in place - generate data set itself
%
% For the most part translated line-by-line from r/s-code.
 
n = years;
m = cohorts;
%sM = t;
 
logs = log(s);
 
if (exist('in_cfactor','var')),
  cfactor = in_cfactor;
else,
  cfactor = 0.65;
end;
 
fprintf(1,'Using cfactor = %s. (This number, best left between 0 and 1\n',num2str(cfactor));
fprintf(1,'determines how much of the stock should be subject to catch. A larger\n');
fprintf(1,'value means more catch.)\n\n');
 
 
%alpha = 0;
%Z = (rand(1,cohorts)) - 0.5;
 
% TRANSPOSE OF ITS USUAL FORM!
N = zeros(n,m);
 
%N(1,:) = N0;
% Added "* 100" 01032010
N(1,:) = N0(:)'; 
 
% SAME HERE, TRANSPOSE!
CC = zeros(size(N));
 
%u = randn(n+m-1,1) * sM;
 
%for (i=2:n),
%  
%  if (i>10),
%    CC(i-1,:) = rand(1,m) * cfactor/n  .* N(i-1,:);
%  end;
%  
%  N(i,:)  = (N(i-1,:) - CC(i-1,:)) * exp(-(M+u(i)));
%  
%end;
 
% NEW loop, replaces the one above
for (j=1:m),
  for (i=2:n),
 
    if (i>2),
      CC(i-1,j) = rand(1,1) * cfactor/n * N(i-1,j);
    end;
 
    N(i,j) = (N(i-1,j) - CC(i-1,j)) * exp(-M);
  end;
end;
 
%fprintf(1,'Population trajectory:\n');
%N'
 
II = cell(surveys,1);
for (i=1:surveys),
 
  % Updated on 26Apr2010 to not include 0 year group
  II{i} = q(i) * N(2:n,:) .* exp(randn(n-1,m) * exp(logs));
end;
 
surveyraw = [];
% now create surveyraw matrix, to be put in file
% Updated on 26Apr2010 to not include 0 year group
for (l=1:surveys),
 
  surveyraw_l = zeros((years-1)*cohorts,4);
  III = II{l};
 
  for (i = 1:cohorts),
    for (j=1:years-1),
 
      surveyraw_l((i-1)*(years-1)+j,:) = [i j l III(j,i)];
 
 
    end;
  end;
 
  surveyraw = [surveyraw; surveyraw_l];
end;
 
% Surveyraw matrix now in place, extract requested  number of rows from
% it randomly
max_possible_observations = (years-1)*cohorts*surveys;
 
if (requested_observations >= max_possible_observations),
  observations = max_possible_observations;
  if (requested_observations > max_possible_observations),
    fprintf(1,'Requested number of observations %s larger than maximum number of\n',num2str(requested_observations))
    fprintf(1,'observations possible (%s). Creating %s observations.\n',num2str(max_possible_observations), ...
            num2str(max_possible_observations));
    fprintf(1,'(No observations of zero age group in this version.)\n');
  end;
  permutation = 1:max_possible_observations;
else,  
  observations = requested_observations;
  fprintf(1,'Including %s of %s possible observations, randomly chosen.\n',num2str(requested_observations), ...
          num2str(max_possible_observations));
  initial_permutation = randperm(max_possible_observations);
  permutation = sort(initial_permutation(1:requested_observations));
end;
 
% Write to file
fid = fopen([filename '.dat'],'w');
% disp(fopen('all'));
fprintf(fid,'# Randomly generated herring file\n');
fprintf(fid,'# True values:\n');
fprintf(fid,'#\tN0: ');
fprintf(fid,'%s ',num2str(N0));
fprintf(fid,'\n');
fprintf(fid,'#\tq: ');
fprintf(fid,'%s ',num2str(q));
fprintf(fid,'\n');
fprintf(fid,'#\ts: %8.5f \n#\tM: %8.5f\n#\n',s,M);
fprintf(fid,'# True population trajectory, transpose(N):\n#\n');
for (i=1:size(N,1)),
  fprintf(fid,'# ');
  for (j=1:size(N,2)),
   fprintf(fid,'%s ',num2str(N(i,j)));
  end;
  fprintf(fid,'\n');
end;
fprintf(fid,'#\n');
fprintf(fid,'\n# m (Antall kohorter)\n %3d \n\n',cohorts);
 
% Endret 01032010 -- gammel kode hadde et sluttår for mye
fprintf(fid,'# Startaar for kohorter\n');
fprintf(fid,'%i ',[0:cohorts-1]+year_offset);
fprintf(fid,'\n# Sluttaar for kohorter\n');
fprintf(fid,'%i ',[years-1:years-1+cohorts-1]+year_offset);
fprintf(fid,'\n# Catches (one cohort per line)\n');
 
for (j=1:cohorts),
  fprintf(fid,'%5.3f ',CC(:,j)');
  fprintf(fid,'\n');
end;
 
fprintf(fid,'\n# Number of surveys\n %3d \n',surveys);
fprintf(fid,'# Number of survey indices\n %4d \n',observations);
fprintf(fid,'# Survey table\n');
 
 
for (i=1:observations),
  fprintf(fid,'%3d %3d %3d %6.3f\n',surveyraw(permutation(i),:));
end;
 
fclose(fid);

%  Now generate start values for optimisation, namely
%  true values as above

% ADMB format
fid = fopen([filename '.pin'],'w');
fprintf(fid,'# N0\n%s\n# q\n%s\n# logs\n%s\n# M\n%s\n', ...
        num2str(N0/1000),num2str(q),num2str(logs),num2str(M));
fclose(fid);


return;
</file>

==== montecarlo_one_cohort_cfactor.m ====
:!: This file is intended only to document how a monte-carlo simulation **can** be done, it will probably need modification to suit most needs.

:!: If matlab gives an error message about gcc libs being missing, try replacing
<code>
unix('./cod');
</code>
with
<code>
unix('unset LD_LIBRARY_PATH; ./cod');
</code>
<file txt montecarlo_one_cohort_cfactor.m>
try,
 more off;
end;

% number of replications per choice for cfactor
numreps = 1;

% cfs contains the cfactor choices
%cfs = [0.3 0.6 0.9 1.2 1.5];
cfs = [5];

ns = length(cfs);

cod_dir = '../cod_admb';
base_dir = pwd;
filename = [cod_dir '/cod'];

results = cell(ns,1);

for (ii=1:ns),
  results{ii} = zeros(4,numreps);

  for (jj=1:numreps),

    % Generate data
    Nt = generate_data([500 ],0.5,0.1,0.15,'../cod_admb/cod',1,15,1,1*14*1,1980,cfs(ii));

    % Run the exp
    cd(cod_dir);
    pause(0.1);
    unix('./cod');
    pause(0.1);
    cd(base_dir);

    % save the data
    try,
      extractpar;
      results{ii}(1,jj) = N0(1); 
      results{ii}(2,jj) = q(1);
      results{ii}(3,jj) = M;
      results{ii}(4,jj) = exp(logs);
    end;
  end;

end;

names = {'N0','q','M','s'};

% Plotting command:
for (i=1:length(cfs)),
  figure; 
  plot(sort(results{i}(1,:))); 
  hold on; plot([1 100], [1 1]*median(sort(results{i}(1,:))),'r'); 
  title(sprintf('cfactor = %s',num2str(cfs(i)))); 
  xlabel('Replication number (sorted)'); 
  ylabel('Estimated N0 and median'); 
  plot([1 100],[0.500 0.500],'g'); 
  legend('Estimates','Median','True N0','Location','NorthWest')
end;

% Example command for saving image
% print(gcf,'-dpng','../images_for_wiki/montecarlo_cfactor_test.png')
</file>

==== read_dat.m ====

<file txt read_dat.m>
function [m,sy,ey,C,ns,I] = read_dat(filename);
% [num_cohorts,start_years,end_years,catch,num_surveys,survey_indices] = read_dat(filename);
%
% Reads the ADMB data file filename ( = e.g. 'cod.dat'); 
% and returns the data as output variables.
%
% Written by Lennart Frimannslund, April 2010.

fid = fopen(filename,'r');
l = 0;
m_done =0;
sy_done = 0;
ey_done = 0;
C_done = 0;
ns_done = 0;
nsi_done = 0;
I_done = 0;
I_counter = 1;

while (l~= -1),

 l = fgetl(fid);

 candidate_number = str2num(l);
 
 if (~isempty(candidate_number) & isfinite(norm(candidate_number)) ),

   if (m_done == 0),
     m = candidate_number;
     m_done = 1;
   elseif (sy_done == 0),
     sy = candidate_number;
     sy_done = 1;
   elseif (ey_done == 0),
     ey = candidate_number;
     ey_done = 1;
   elseif (C_done == 0),
     C = [candidate_number; zeros(m-1,length(candidate_number))];
     for (i = 2:m),
       l = fgetl(fid);
       candidate_number = str2double(l);
       C(i,1:length(candidate_number)) = candidate_number;
     end;
     C_done = 1;
   elseif (ns_done == 0),
     ns = candidate_number;
     ns_done = 1;
   elseif (nsi_done == 0),
     num_survey_indices = candidate_number;
     if (num_survey_indices == 0),
       I = [];
       break;
     else,
       I = zeros(candidate_number,4);
       nsi_done = 1;
     end;
   elseif (I_done == 0),
     I(I_counter,:) = candidate_number;
     if (I_counter == size(I,1)),
       break;
     else,
       I_counter = I_counter + 1;
     end;
   end;
 end;

 if (isempty(l)),
   l = 0;
 end;

end;
%fprintf(1,'Reading finished!\n');
fclose(fid);
</file>

==== strip_hash_and_colon.m ====

This is a helper function called by extractpar. The user should not need to call this function directly.
<file txt strip_hash_and_colon.m>
function [out] = strip_hash_and_colon(input);
% [out] = strip_hash_and_colon(input)
%
% Helper function for extractpar.
%
% Input  - string
% Output - the same string, with # and : converted to spaces.
%
% Written by Lennart Frimannslund, April 2010

 out = input;
 for (i=1:length(out)), 

    if (out(i) == '#' | out(i) == ':'),
      out(i) = ' ';
    end;
 end;

return;
</file>

==== table2admb.m ====
:!: **NOTE: Generates .pin file for simple model!**
<file txt table2admb.m>
% table2admb(catchfile,surveyfiles,firstyear,lastyear, maxage, num_cohorts, outfilename)
%
% Conversion of ICES table data to the .dat format used by AD Model builder.
%
% Input:
%
% catchfile   - string with path to ascii file containing catch data
% surveyfiles - either string with path to ascii file containing
%               survey data, or cell array of such strings.
% firstyear   - first year of observation/estimation (e.g 1995)
% lastyear    - last year of observation (e.g 2001)
% maxage      - maximum age for cohort (after maxage is reached no more
%               computation is done for the cohort in question).
% num_cohorts - number of cohorts to follow  (e.g. 4)
%
% NOTE: first year, last year, num_cohorts and maxage should have
%       consistent values. We always include the cohort born in 
%       firstyear, and at the moment always have cohorts born in
%       consecutive years.
%
% Output:
%
% outfilename  - (Optional) string with filename of output file without suffix,
%                that is, outfilename = '../cod' will create ../cod.dat and
%                ../cod.pin. .pin-file generation may be omitted in the future.
%                If outfilename is not defined, then output is written to the 
%                screen.
%
%     
% Written by Lennart Frimannslund, April 2010 

% Test input variables
if (in_maxage > in_lastyear-firstyear),
  fprintf(1,'Maxage = %i, but lastyear-firstyear = %i. Setting maxage to %i.\n',maxage,in_lastyear-firstyear,in_lastyear-firstyear);
  maxage = lastyear-firstyear;
else,
  maxage = in_maxage;
end;

if (in_lastyear > firstyear + (num_cohorts-1) + maxage),
  fprintf(1,'Lastyear = %i, but firstyear = %i, num_cohorts = %i and maxage = %i.\nSetting lastyear to %i.\n', ...
          in_lastyear,firstyear,num_cohorts,maxage,firstyear + (num_cohorts-1) + maxage);
  lastyear = firstyear + (num_cohorts-1) + maxage;
else,
 lastyear = in_lastyear;
end;


if (num_cohorts > lastyear-firstyear),
  fprintf(2,'Error: firstyear = %i, lastyear = %i, maxage = %i and num_cohorts = %i not consistent.\nExiting...\n', ...
          firstyear, lastyear, maxage,num_cohorts);

end;

% Read input files
C = load(catchfile,'-ascii');
if (iscell(surveyfiles)),
  S = cell(length(surveyfiles),1);
  for (i=1:length(surveyfiles)),
    S{i} = load(surveyfiles{i},'-ascii');
  end;
else,
  S = load(surveyfiles,'-ascii');
end;

% Create output file, or OVERWRITE it if it exists.
if (exist('outfilename','var') == 0),
 fid = 1;
else,
 try,
  fid = fopen([outfilename '.dat'],'w');
 catch,
  disp(lasterr);
  fid = 1;
 end;
end;

% Print data to file, item by item

% Headeer and number of cohorts
fprintf(fid,'#Generated .dat-file with catch and acoustic/trawl data from ICES report.\n#\n#Number of cohorts\n');
fprintf(fid,'%i\n',num_cohorts);

% Birth years
fprintf(fid,'\n# Cohort birth years (Startaar)\n');
for (i=1:num_cohorts),
 fprintf(fid,'%4i ',firstyear+i-1);
end;

% Last year for each cohorst
fprintf(fid,'\n\n# Last year of observation (Sluttaar)\n');
for (i=1:num_cohorts),
 fprintf(fid,'%4i ',min(firstyear+i-1+maxage,lastyear));
end;

% Catch. Note that catch in the table is in thousands, but here we divide the volume by 1000 so that it is in millions
fprintf(fid,'\n\n# Catch data IN MILLIONS. One cohort per row, one year per column.\n# First column is the birth year. Last column is not used by .tpl file\n# but must be present.\n');

youngest_age_caught = C(1,2);
max_age_in_catch_table = max(C(1,:));

% Populate catch matrix, one cohort at a time
for (i=1:num_cohorts),
  [this_birthyear_index] = find(C(:,1) == firstyear + i - 1);

  % Set catch to zero for age groups too young to be in the table
  for (j=0:youngest_age_caught - 1),
    fprintf(fid,'0.0 ');
  end;

  % Go diagonally in the catch table until either the end is reached or maxage is reached
  for (j=youngest_age_caught:min(max_age_in_catch_table,min(maxage,lastyear-firstyear - i + 1))),
    % Make sure we don't exceed size of table
    desired_row = this_birthyear_index+j-youngest_age_caught;
    if (desired_row <= size(C,1)),
      % Divide by 1000 since catch table data is in thousands, not millions
      fprintf(fid,'%f ',C(desired_row,j-youngest_age_caught+2)/1000);
    end;
  end;
  fprintf(fid,'\n');
end;

% Surveys
fprintf(fid,'\n# Survey section. Surveys included:\n#\n');
if (iscell(S)),
  for (i=1:length(S)),
    fprintf(fid,'# %s\n',surveyfiles{i});
  end;
  fprintf(fid,'#\n\n# Number of surveys per year\n%i\n',length(S));
else,
  fprintf(fid,'# %s\n',surveyfiles);
  fprintf(fid,'#\n\n# Number of surveys per year\n1\n');
end;

% For the surveys we don't yet know how many survey indices there will be,
% even though the data format requires this number to come first.
% We first write the survey data to a matrix of maximum possible size, and then to file
if (iscell(S)),
  num_surveys = length(S);
else,
  num_surveys = 1;
end;


I = zeros(num_cohorts*num_surveys*(lastyear-firstyear+1),4);


num_survey_indices = 0;

for (k=1:num_surveys),
  
   if (num_surveys == 1),
     this_survey = S;
   else, 
     this_survey = S{k};
   end;

  % Populate one cohort at a time
  for (i=1:num_cohorts),
    % Find the row in column 2 which corresponds to the current cohort
    this_cohort_row_in_col_2 = find( (this_survey(:,1) - this_survey(1,2) - (firstyear + (i-1)) ) == 0);

    last_legal_row = find(this_survey(:,1) == lastyear);
    last_legal_col = find(this_survey(1,:) == lastyear-firstyear - i + 1);
    max_j = min(last_legal_row-this_cohort_row_in_col_2,last_legal_col-2);

    % go diagonally, same as for catch.
    for (j=0:max_j),
      try,
       candidate_volume = this_survey(this_cohort_row_in_col_2+j,2+j);
        if (isfinite(candidate_volume)),
          % If we encounter an observation of zero (which will make the 
          % objective function infinity), don't include it
          if (candidate_volume > 0.0),
             num_survey_indices = num_survey_indices + 1;
             I(num_survey_indices,:) = [ i this_survey(1,2+j) k candidate_volume];
          end;
        end;

      % For debug purposes, report if our indices are out of bounds. Shouldn't happen anymore.
      catch,
        fprintf(1,'Survey table index exceeded. (i,j) = (%i,%i), but size S{%i} = (%i,%i).\n', ...
           this_cohort_row_in_col_2+j,this_cohort_row_in_col_2+1+j,k,size(this_survey,1),size(this_survey,2));
      end;
    end;
   
  end;
 
end;

% Print everything to file
fprintf(fid,'\n# Number of survey indices\n');
fprintf(fid,'%i \n',num_survey_indices);
fprintf(fid,'\n# Survey table: (Cohort, age, survey_number, observed_volume)\n');
for (i=1:num_survey_indices),
  fprintf(fid,'%3i %3i %3i %f\n',I(i,1),I(i,2),I(i,3),I(i,4));
end;

fprintf(fid,'# End of file');

% Done, close file if applicable
if (fid ~= 1),
  fclose(fid);
end;

%return;
% Make .pin-file too
if (exist('outfilename','var') == 0),
 fid = 1;
else,
 try,
  fid = fopen([outfilename '.pin'],'w');
 catch,
  disp(lasterr);
  fid = 1;
 end;
end;

fprintf(fid,'\n\n# Autogenerated .pin file\n# N0\n');
for (i=1:num_cohorts),

   fprintf(fid,'%f  ',1e3);

end;
fprintf(fid,'\n\n# q\n');

for (i=1:num_surveys),

  fprintf(fid,'%4.2f ',0.4);

end;

fprintf(fid,'\n\n# log(s)\n0.1\n\n# M\n0.2\n');

if (fid ~= 1),
  fclose(fid);
end;
</file>

====== Files written autumn 2010 ======

==== table2allCohortsModel.m ====

<file txt table2allCohortsModel.m>
function table2allCohortsModel(startyear, endyear, minage, maxage, catchfile, ...
		    surveyfile, datfile, include_weight_maturity);
% table2allCohortsModel(startyear, endyear, minage, maxage, catchfile, surveyfile, datfile, include_weight_and_maturity);


fid = fopen(datfile,'w');

try,
  fprintf(fid,'# Auto-generated datfile\n\n');

  fprintf(fid,'# Startyear\n%i\n\n# Endyear\n%i\n\n',startyear,endyear);

  fprintf(fid,'# Min age\n%i\n\n# Max age\n%i\n\n',minage,maxage);

  fprintf(fid,'# Catch\n');

  C = load(catchfile);

  startyear_row = find(C(:,1) == startyear);
  endyear_row = find(C(:,1) == endyear);
  
  minage_col = find(C(1,:) == minage);
  maxage_col = find(C(1,:) == maxage);
  
  for (i=startyear_row:endyear_row),
    for (j=minage_col:maxage_col),
      fprintf(fid,'%10.6f ',C(i,j));
    end;
    fprintf(fid,'\n');
  end;
  
  
  fprintf(fid,'# Survey data\n');
  S = load(surveyfile);
  
  startyear_row = find(S(:,1) == startyear);
  endyear_row = find(S(:,1) == endyear);
  
  minage_col = find(S(1,:) == minage);
  maxage_col = find(S(1,:) == maxage);
  
  for (i=startyear_row:endyear_row),
    for (j=minage_col:maxage_col),
      fprintf(fid,'%10.6f ',S(i,j));
    end;
    fprintf(fid,'\n');
  end;
  
  if (exist('include_weight_maturity','var') && include_weight_maturity == 1),
    
    mature_3_12_2009;
    weight_3_11_2009;
    mature = mature';
    weight = weight';
    
    fprintf(fid,'# Proportion mature at age\n');
    startyear_row = find(mature(:,1) == startyear);
    endyear_row = find(mature(:,1) == endyear);
    
    minage_col = find(mature(1,:) == minage);
    maxage_col = find(mature(1,:) == maxage);
    
    for (i=startyear_row:endyear_row),
      for (j=minage_col:maxage_col),
	fprintf(fid,'%10.6f ',mature(i,j));
      end;
      fprintf(fid,'\n');
    end;
    
    
    
    fprintf(fid,'# Weight at age\n');
    startyear_row = find(weight(:,1) == startyear);
    endyear_row = find(weight(:,1) == endyear);
    
    minage_col = find(weight(1,:) == minage);
    maxage_col = find(weight(1,:) == maxage);
    
    for (i=startyear_row:endyear_row),
      for (j=minage_col:maxage_col),
	fprintf(fid,'%10.6f ',weight(i,j));
      end;
      fprintf(fid,'\n');
    end;
    
    get_sondre_data;
    csd = [nan sondre_ages; sondre_years' sondre_sigma];
    
    fprintf(fid,'# Catch standard deviations\n');
    startyear_row = find(csd(:,1) == startyear);
    endyear_row = find(csd(:,1) == endyear);
    
    minage_col = find(csd(1,:) == minage);
    maxage_col = find(csd(1,:) == maxage);
    
    for (i=startyear_row:endyear_row),
      for (j=minage_col:maxage_col),
	fprintf(fid,'%10.6f ',csd(i,j));
      end;
      fprintf(fid,'\n');
    end;
    
    
  end;
  
  
  fclose(fid);

catch,
  fprintf(2,'%s\n',lasterr);
  fclose(fid);
end;


</file>