Differences

This shows you the differences between two versions of the page.

--- mfiles:mfiles [2010/05/04 08:03] – lennartfr
+++ mfiles:mfiles [2010/12/13 15:34] (current) – lennartfr
@@ Line 1: / Line 1: @@
 ====== 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 generation of .pin files //may// be a bit faulty, in particular when running other models than the simple model. These files can be modified manually. Adding zeros to the end of the .pin file until ADMB runs is usually sufficient.
+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.
+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,
@@ Line 18: / Line 22: @@
 % Enter the proper directory
-cd ~lennartfr/common_files/mfiles
+%cd ~lennartfr/common_files/mfiles
 % Define data files
@@ Line 41: / Line 45: @@
 % Run ADMB
-unix('./cod');
+unix('unset LD_LIBRARY_PATH; ./cod');
 % Return to mfile directory and plot
@@ Line 47: / Line 51: @@
 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 ====
@@ Line 89: / Line 143: @@
 ==== datpar2plot.m ====
 <file txt datpar2plot.m>
+function datpar2plot(filename,new_figure);
 % datpar2plot(filename,make_new_figure)
 %
@@ Line 111: / Line 166: @@
 %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);
@@ Line 142: / Line 198: @@
   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(firstyear+(i-1),N(i,1),std_coeff * std_N0(i,2),colors{i});
     %
-    % ERRORBAR doesn't quite work in octave, code below makes our own bar
+    % 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{2}]);
+    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));i
+    %fprintf(2,'Errorbar plot failed for cohort born in %4i -- wrong variable name?\n',firstyear+(i-1));
-    disp(lasterr);
+    %disp(lasterr);
     errorbar(firstyear+(i-1),N(i,1),std_coeff * std_N0(i,2));
   end;
@@ Line 160: / Line 216: @@
 for (i=1:size(survey,1)),
-  % Check if we have one or several qs
+  % Check if q_alpha is present, which indicates sigmoid function
-  if (exist('q','var')),
+  if (exist('q_alpha','var')),
-    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)}]);
+    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)/q0,[colors{survey(i,1)} dotstyles{survey(i,3)}]);
+    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;
@@ Line 173: / Line 242: @@
 ylabel('Volume in billions');
-curr_axis = axis;
+%curr_axis = axis;
-axis([curr_axis(1) curr_axis(2) 0 2500/N0scale curr_axis(5) curr_axis(6)]);
+%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;
@@ Line 250: / Line 338: @@
   % Send the variable name and values to the command line
   eval(sprintf('%s = [ %s ];',varname,l));
+  %sprintf('%s = [ %s ];',varname,l)
 end;
@@ Line 270: / Line 358: @@
   output = fgetl(fid);
-  N_cell{i} = str2double(output);
+  N_cell{i} = str2num(output);
   max_width_so_far = max(max_width_so_far,length(N_cell{i}));
@@ Line 354: / Line 442: @@
 This file is used for generating synthetic data.
-:!: **NOTE: Does not generate .pin file!**
+:!: **NOTE: Generates .pin file for simple model.**
 FIXME **NOTE: The line**
@@ Line 589: / Line 677: @@
 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>
@@ Line 694: / Line 792: @@
  l = fgetl(fid);
- candidate_number = str2double(l);
+ candidate_number = str2num(l);
  if (~isempty(candidate_number) & isfinite(norm(candidate_number)) ),
@@ Line 742: / Line 841: @@
 end;
+%fprintf(1,'Reading finished!\n');
 fclose(fid);
 </file>
@@ Line 1011: / Line 1110: @@
   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>