<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40"><head><meta http-equiv=Content-Type content="text/html; charset=us-ascii"><meta name=Generator content="Microsoft Word 14 (filtered medium)"><style><!--
/* Font Definitions */
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
span.EmailStyle17
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:#1F497D;}
.MsoChpDefault
{mso-style-type:export-only;
font-family:"Calibri","sans-serif";}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]--></head><body lang=EN-US link=blue vlink=purple><div class=WordSection1><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Hi,<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>In ADMB the indices of vector operations must be consistent. I.e., if you subtract a (sub)vector of X which has dimension 1,10 e.g., :<o:p></o:p></span></p><p class=MsoNormal style='text-indent:.5in'><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>X(2,10) – X(1,9); <o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>will return an error. A shift operator is required (page 15-10 in my manual). A shorthand way to do this is:<o:p></o:p></span></p><p class=MsoNormal style='text-indent:.5in'><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>X(2,10) – ++X(1,9);<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>A long hand way is:<o:p></o:p></span></p><p class=MsoNormal style='text-indent:.5in'><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>X(2,10) – X(1,9).shift(1);<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Cheers,<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Jim<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'> <o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'> <o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>James Ianelli<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><a href="http://www.afsc.noaa.gov/REFM/stocks/assessments.htm"><span style='color:blue'>REFM Division, Alaska Fisheries Science Center</span></a><o:p></o:p></span></p><p class=MsoNormal><i><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>National Marine Fisheries Service, National Oceanic and Atmospheric Administration<o:p></o:p></span></i></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>7600 Sand Point Way NE<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Seattle WA 98115<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'>Ph. 206 526 6510<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:11.0pt;font-family:"Calibri","sans-serif";color:#1F497D'><o:p> </o:p></span></p><div style='border:none;border-left:solid blue 1.5pt;padding:0in 0in 0in 4.0pt'><div><div style='border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0in 0in 0in'><p class=MsoNormal><b><span style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'>From:</span></b><span style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'> users-bounces@admb-project.org [mailto:users-bounces@admb-project.org] <b>On Behalf Of </b>Luis Ridao<br><b>Sent:</b> Thursday, October 07, 2010 7:45 AM<br><b>To:</b> users@admb-project.org<br><b>Subject:</b> [ADMB Users] error when imposing penalty on objective function<o:p></o:p></span></p></div></div><p class=MsoNormal><o:p> </o:p></p><p class=MsoNormal>ADMB-help,<o:p></o:p></p><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>This is my first e-mail to this mailing-list. I'm new to ADMB and <o:p></o:p></p></div><div><p class=MsoNormal>trying to learn little by little.<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>I'm working on a catch-at-age model provided in the admb manual.<o:p></o:p></p></div><div><p class=MsoNormal>The first step was to import my personal stock data set. As a second step <o:p></o:p></p></div><div><p class=MsoNormal>I wish to impose a penalty on the objective function based on the differences in the fishing mortality<o:p></o:p></p></div><div><p class=MsoNormal>coefficient ("log_fy_coff")<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>ADMB fails to run the model with the following message:<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>"Incompatible bounds in prevariable operator * (_CONST dvar_vector& v1,_CONST dvar_vector& v2)"<o:p></o:p></p></div><div><p class=MsoNormal>The line causing the error is the FUNCTION "evaluate_the_objective_function" section at the of the file<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal># originally (from the manual admb.pdf)<o:p></o:p></p></div><div><div><p class=MsoNormal>FUNCTION evaluate_the_objective_function<o:p></o:p></p></div><div><p class=MsoNormal> // penalty functions to ``regularize '' the solution<o:p></o:p></p></div><div><p class=MsoNormal> f+=.01*norm2(log_relpop);<o:p></o:p></p></div><div><p class=MsoNormal> avg_F=sum(F)/double(size_count(F));<o:p></o:p></p></div><div><p class=MsoNormal> if (last_phase())<o:p></o:p></p></div><div><p class=MsoNormal> {<o:p></o:p></p></div><div><p class=MsoNormal> // a very small penalty on the average fishing mortality<o:p></o:p></p></div><div><p class=MsoNormal> f+= .001*square(log(avg_F/.2));<o:p></o:p></p></div><div><p class=MsoNormal> }<o:p></o:p></p></div><div><p class=MsoNormal> else<o:p></o:p></p></div><div><p class=MsoNormal> {<o:p></o:p></p></div><div><p class=MsoNormal> f+= 1000.*square(log(avg_F/.2));<o:p></o:p></p></div><div><p class=MsoNormal> }<o:p></o:p></p></div><div><p class=MsoNormal> f+=0.5*double(size_count(C)+size_count(log_fy_coff)) <o:p></o:p></p></div><div><p class=MsoNormal> * log( sum(elem_div(square(C-obs_catch_at_age),.01+C))<o:p></o:p></p></div><div><p class=MsoNormal> + 100*norm2(log_fy_coff))); // ORIGINAL LINE <o:p></o:p></p></div></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal># the modification of the original code (last line)<o:p></o:p></p></div><div><div><p class=MsoNormal>FUNCTION evaluate_the_objective_function<o:p></o:p></p></div><div><p class=MsoNormal> // penalty functions to ``regularize '' the solution<o:p></o:p></p></div><div><p class=MsoNormal> f+=.01*norm2(log_relpop);<o:p></o:p></p></div><div><p class=MsoNormal> avg_F=sum(F)/double(size_count(F));<o:p></o:p></p></div><div><p class=MsoNormal> if (last_phase())<o:p></o:p></p></div><div><p class=MsoNormal> {<o:p></o:p></p></div><div><p class=MsoNormal> // a very small penalty on the average fishing mortality<o:p></o:p></p></div><div><p class=MsoNormal> f+= .001*square(log(avg_F/.2));<o:p></o:p></p></div><div><p class=MsoNormal> }<o:p></o:p></p></div><div><p class=MsoNormal> else<o:p></o:p></p></div><div><p class=MsoNormal> {<o:p></o:p></p></div><div><p class=MsoNormal> f+= 1000.*square(log(avg_F/.2));<o:p></o:p></p></div><div><p class=MsoNormal> }<o:p></o:p></p></div><div><p class=MsoNormal> f+=0.5*double(size_count(C)+size_count(log_fy_coff)) <o:p></o:p></p></div><div><p class=MsoNormal> * log( sum(elem_div(square(C-obs_catch_at_age),.01+C))<o:p></o:p></p></div><div><p class=MsoNormal> + 100*norm2(log_fy_coff(2,nyrs)-log_fy_coff(1,nyrs-1))); // MODIFIED LINE<o:p></o:p></p></div></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>The .tpl file is provided below.<o:p></o:p></p></div><div><p class=MsoNormal>I'm running on a Windows machine.<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>########################################################################################<o:p></o:p></p></div><div><div><p class=MsoNormal>DATA_SECTION<o:p></o:p></p></div><div><p class=MsoNormal> init_int nyrs // the number of years odf data<o:p></o:p></p></div><div><p class=MsoNormal> init_int nages // the number of age classess in the population<o:p></o:p></p></div><div><p class=MsoNormal> init_matrix obs_catch_at_age(1,nyrs,1,nages) // observed catch-at-age data<o:p></o:p></p></div><div><p class=MsoNormal> init_number M // estimate of natural mortality rate<o:p></o:p></p></div><div><p class=MsoNormal> init_vector relwt(2,nages); // need to have relative weight-at-age to calculate B2+<o:p></o:p></p></div><div><p class=MsoNormal> vector ages(1,nages); // ages of data<o:p></o:p></p></div><div><p class=MsoNormal> vector ages4plus(1,nages-1); // non-recruiting ages in the population<o:p></o:p></p></div><div><p class=MsoNormal> vector years(1,nyrs); // years of data<o:p></o:p></p></div><div><p class=MsoNormal> int pred_year; // prediction year<o:p></o:p></p></div><div><p class=MsoNormal>INITIALIZATION_SECTION<o:p></o:p></p></div><div><p class=MsoNormal> //log_q -1 // original -1<o:p></o:p></p></div><div><p class=MsoNormal> log_popscale 5 // original 5<o:p></o:p></p></div><div><p class=MsoNormal>PARAMETER_SECTION<o:p></o:p></p></div><div><p class=MsoNormal> //init_number log_q(1) // log-catchability<o:p></o:p></p></div><div><p class=MsoNormal> init_number log_popscale(1) // overall population scaling parameter<o:p></o:p></p></div><div><p class=MsoNormal> init_bounded_dev_vector log_sel_coff(1,nages-1,-15.,15.,2) // original log_sel_coff(1,nages-1,-15.,15.,2) <o:p></o:p></p></div><div><p class=MsoNormal> init_bounded_dev_vector log_relpop(1,nyrs+nages-1,-15.,15.,2) // original log_relpop(1,nyrs+nages-1,-15.,15.,2) <o:p></o:p></p></div><div><p class=MsoNormal> init_bounded_dev_vector log_fy_coff(1,nyrs,-.3,.3,3) // original log_fy_coff(1,nyrs,-2.,2.,3) <o:p></o:p></p></div><div><p class=MsoNormal> vector log_sel(1,nages)<o:p></o:p></p></div><div><p class=MsoNormal> vector log_fy(1,nyrs)<o:p></o:p></p></div><div><p class=MsoNormal> vector log_initpop(1,nyrs+nages-1);<o:p></o:p></p></div><div><p class=MsoNormal> matrix F(1,nyrs,1,nages) // instantaneous fishing mortality<o:p></o:p></p></div><div><p class=MsoNormal> matrix Z(1,nyrs,1,nages) // instantaneous total mortality<o:p></o:p></p></div><div><p class=MsoNormal> matrix S(1,nyrs,1,nages) // survival rate<o:p></o:p></p></div><div><p class=MsoNormal> matrix N(1,nyrs,1,nages) // predicted numbers-at-age<o:p></o:p></p></div><div><p class=MsoNormal> matrix C(1,nyrs,1,nages) // predicted catch-at-age<o:p></o:p></p></div><div><p class=MsoNormal> objective_function_value f<o:p></o:p></p></div><div><p class=MsoNormal> number recsum<o:p></o:p></p></div><div><p class=MsoNormal> number initsum<o:p></o:p></p></div><div><p class=MsoNormal> sdreport_number avg_F<o:p></o:p></p></div><div><p class=MsoNormal> sdreport_vector predicted_N(2,nages)<o:p></o:p></p></div><div><p class=MsoNormal> sdreport_vector ratio_N(2,nages)<o:p></o:p></p></div><div><p class=MsoNormal> // changed from the manual because adjusted likelihood routine doesn't<o:p></o:p></p></div><div><p class=MsoNormal> // work<o:p></o:p></p></div><div><p class=MsoNormal> likeprof_number pred_B<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>PRELIMINARY_CALCS_SECTION<o:p></o:p></p></div><div><p class=MsoNormal> ages.fill_seqadd(3,1); // vector of ages<o:p></o:p></p></div><div><p class=MsoNormal> ages4plus.fill_seqadd(4,1); // vector of non recruiting ages<o:p></o:p></p></div><div><p class=MsoNormal> years.fill_seqadd(1975,1); // fill vector of years with years<o:p></o:p></p></div><div><p class=MsoNormal> pred_year=years[nyrs]+1; // year of prediction = last_year +1<o:p></o:p></p></div><div><p class=MsoNormal>PROCEDURE_SECTION<o:p></o:p></p></div><div><p class=MsoNormal> // example of using FUNCTION to structure the procedure section<o:p></o:p></p></div><div><p class=MsoNormal> get_mortality_and_survivial_rates();<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal> get_numbers_at_age();<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal> get_catch_at_age();<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal> evaluate_the_objective_function();<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>FUNCTION get_mortality_and_survivial_rates<o:p></o:p></p></div><div><p class=MsoNormal> int i, j;<o:p></o:p></p></div><div><p class=MsoNormal> //calculate the selectivity from the sel_coffs ---------------------<o:p></o:p></p></div><div><p class=MsoNormal> for (j=1;j<nages;j++)<o:p></o:p></p></div><div><p class=MsoNormal> {<o:p></o:p></p></div><div><p class=MsoNormal> log_sel(j)=log_sel_coff(j);<o:p></o:p></p></div><div><p class=MsoNormal> }<o:p></o:p></p></div><div><p class=MsoNormal> //the selectivity is the same for the last two age classes<o:p></o:p></p></div><div><p class=MsoNormal> log_sel(nages)=log_sel_coff(nages-1);<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal> for (i=1;i<=nyrs;i++)<o:p></o:p></p></div><div><p class=MsoNormal> {<o:p></o:p></p></div><div><p class=MsoNormal> log_fy(i)=log_fy_coff(i);<o:p></o:p></p></div><div><p class=MsoNormal> }<o:p></o:p></p></div><div><p class=MsoNormal> F=outer_prod(mfexp(log_fy),mfexp(log_sel)); // F=outer_prod(mfexp(log_fy),mfexp(log_sel)) ; F=outer_prod(mfexp(log_q)*effort,mfexp(log_sel));<o:p></o:p></p></div><div><p class=MsoNormal> Z=F+M;<o:p></o:p></p></div><div><p class=MsoNormal> // get the survival rate<o:p></o:p></p></div><div><p class=MsoNormal> S=mfexp(-1.0*Z);<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>FUNCTION get_numbers_at_age<o:p></o:p></p></div><div><p class=MsoNormal> int i, j;<o:p></o:p></p></div><div><p class=MsoNormal> <o:p></o:p></p></div><div><p class=MsoNormal> log_initpop=log_relpop+log_popscale;<o:p></o:p></p></div><div><p class=MsoNormal> <o:p></o:p></p></div><div><p class=MsoNormal> for (i=1;i<=nyrs;i++)<o:p></o:p></p></div><div><p class=MsoNormal> {<o:p></o:p></p></div><div><p class=MsoNormal> N(i,1)=mfexp(log_initpop(i));<o:p></o:p></p></div><div><p class=MsoNormal> }<o:p></o:p></p></div><div><p class=MsoNormal> for (j=2;j<=nages;j++)<o:p></o:p></p></div><div><p class=MsoNormal> {<o:p></o:p></p></div><div><p class=MsoNormal> N(1,j)=mfexp(log_initpop(nyrs+j-1));<o:p></o:p></p></div><div><p class=MsoNormal> }<o:p></o:p></p></div><div><p class=MsoNormal> for (i=1;i<nyrs;i++)<o:p></o:p></p></div><div><p class=MsoNormal> {<o:p></o:p></p></div><div><p class=MsoNormal> for (j=1;j<nages;j++)<o:p></o:p></p></div><div><p class=MsoNormal> {<o:p></o:p></p></div><div><p class=MsoNormal> N(i+1,j+1)=N(i,j)*S(i,j);<o:p></o:p></p></div><div><p class=MsoNormal> }<o:p></o:p></p></div><div><p class=MsoNormal> }<o:p></o:p></p></div><div><p class=MsoNormal> // calculated predicted numbers at age for next year<o:p></o:p></p></div><div><p class=MsoNormal> for (j=1;j<nages;j++)<o:p></o:p></p></div><div><p class=MsoNormal> {<o:p></o:p></p></div><div><p class=MsoNormal> predicted_N(j+1)=N(nyrs,j)*S(nyrs,j);<o:p></o:p></p></div><div><p class=MsoNormal> ratio_N(j+1)=predicted_N(j+1)/N(1,j+1);<o:p></o:p></p></div><div><p class=MsoNormal> }<o:p></o:p></p></div><div><p class=MsoNormal> // calculated predicted Biomass for next year for<o:p></o:p></p></div><div><p class=MsoNormal> // adjusted profile likelihood<o:p></o:p></p></div><div><p class=MsoNormal> pred_B=(predicted_N * relwt);<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>FUNCTION get_catch_at_age<o:p></o:p></p></div><div><p class=MsoNormal> C=elem_prod(elem_div(F,Z),elem_prod(1.-S,N));<o:p></o:p></p></div><div><p class=MsoNormal><o:p> </o:p></p></div><div><p class=MsoNormal>FUNCTION evaluate_the_objective_function<o:p></o:p></p></div><div><p class=MsoNormal> // penalty functions to ``regularize '' the solution<o:p></o:p></p></div><div><p class=MsoNormal> f+=.01*norm2(log_relpop);<o:p></o:p></p></div><div><p class=MsoNormal> avg_F=sum(F)/double(size_count(F));<o:p></o:p></p></div><div><p class=MsoNormal> if (last_phase())<o:p></o:p></p></div><div><p class=MsoNormal> {<o:p></o:p></p></div><div><p class=MsoNormal> // a very small penalty on the average fishing mortality<o:p></o:p></p></div><div><p class=MsoNormal> f+= .001*square(log(avg_F/.2));<o:p></o:p></p></div><div><p class=MsoNormal> }<o:p></o:p></p></div><div><p class=MsoNormal> else<o:p></o:p></p></div><div><p class=MsoNormal> {<o:p></o:p></p></div><div><p class=MsoNormal> f+= 1000.*square(log(avg_F/.2));<o:p></o:p></p></div><div><p class=MsoNormal> }<o:p></o:p></p></div><div><p class=MsoNormal> f+=0.5*double(size_count(C)+size_count(log_fy_coff)) <o:p></o:p></p></div><div><p class=MsoNormal> * log( sum(elem_div(square(C-obs_catch_at_age),.01+C))<o:p></o:p></p></div><div><p class=MsoNormal> + 100*norm2(log_fy_coff(2,nyrs)-log_fy_coff(1,nyrs-1))); <o:p></o:p></p></div><div><p class=MsoNormal>########################################################################################<o:p></o:p></p></div></div></div></div></body></html>