[ADMB Users] no assignment operator for init_bounded_number_vector
Larry Jacobson
larry.jacobson at noaa.gov
Fri Jul 19 10:43:01 PDT 2013
1) The idea is to ask the user for an arithmetic starting value which is
used to initialize a log scale parameter (unbeknownst to the user who
does not think in logs).
2) I have a vector of doubles declared in the DATA_SECTION called
srx_init(1,2). See line 193 in attached tpl file.
3) I declare a bounded_number_vector in the PARAMETER_SECTION called
logsrx(1,2). See line 283 in tpl file.
4) In the PRELIMINARY_CALCS_SECTION (line 410) I want to initialize the
bounded_number_vector with code like:
logsrx=log(srx_init);
5) However, the compiler complains:
C:\Users\ljacobso\Documents\Assessments\Shrimp_2013\CSA-v4\july18_2013>cl /c
/EHsc /DUSE_LAPLACE /DWIN32 /DSAFE_ALL /D
__MSVC32__=8 /I. /I"C:\Program Files\ADMB"\include /I"C:\Program
Files\ADMB"\contrib\include csav41.cpp
Microsoft (R) 32-bit C/C++ Optimizing Compiler Version 16.00.30319.01
for 80x86
Copyright (C) Microsoft Corporation. All rights reserved.
csav41.cpp
csav41.cpp(362) : error C2679: binary '=' : no operator found which
takes a right-hand operand of type 'dvector' (or the
re is no acceptable conversion)
C:\Program Files\ADMB\include\admodel.h(2700): could be
'param_init_bounded_number_vector ¶m_init_bounded_nu
mber_vector::operator =(const param_init_bounded_number_vector &)'
while trying to match the argument list
'(param_init_bounded_number_vector, dvector)'
6) I've tried doing the assignment in a loop like:
for (i=1;i<=nsurveys;i++) logsrx(i)=log(srx_init(i));
using an intermediate dvar_vector like:
dvar_vector temp(1,2)=srx_init;
logsrxi=srxinit;
and many other hacks (trying to find a combination of data types with an
assignment operator).
6) Any ideas or suggestions?
Cheers and thanks!
--
**********************
Opinions that may be expressed herein are my own
and not official positions of NOAA or the National
Marine Fisheries Service.
**********************
Larry Jacobson
National Marine Fisheries Service
Northeast Fisheries Science Center
166 Water Street
Woods Hole, MA 02543-1026
Voice: 508-495-2317
Fax: 508-495-2393
E-mail: larry.jacobson at noaa.gov
**********************
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.admb-project.org/pipermail/users/attachments/20130719/30f3a6eb/attachment.html>
-------------- next part --------------
// CSA Version 4.1
//
// 18 July 2013 - Version 4.1 - A. Seaver
// 5 November 2012 - Version 4.0 - A.Seaver
//
GLOBALS_SECTION
#include <admodel.h>
#include <time.h>
#include <iostream>
using namespace std;
// a macro that writes the name and value of any variable to the screen
#define EO(x) cout << #x << "=" << x << endl;
char dtstring[12];
char tmstring[6];
ofstream basicMCMC("csaMCMC.txt");
DATA_SECTION
// First & Last year in Catch
init_int nfyear;
init_int nxyear;
int nyears;
!! nyears = nxyear - nfyear + 1;
// Number of Surveys
init_int nsurvey;
// Catch Data (Landings & Discards)
init_matrix total_catch(1,nyears,1,2);
// Number of Observations
init_int nobs;
// Survey Observations
//
// column 1 - Survey index
// column 2 - year
// column 3 - observation type 1 = Recruit, 2 = Post-Recruit, 3 = Aggregate Surveys 4 = CPUE Surveys
// column 4 - Observed Value
// column 5 - Coefficient of variation
// column 6 - time of survey obervation ts >= 0.0 ts < 1.0
// column 7 - Selectivity of recruitment - (recruits only)
//
init_matrix survey_data(1,nobs,1,7);
// Natural Mortality
init_vector nmort(1,nyears);
// CV on Catch
init_number cvcatch;
// Weighting Factor on Catch
init_number lambda_catch;
// Survey Specific weight factor - use to turn off surveys or to give a survey greater emphasis
init_vector lambda_survey(1,nsurvey);
// Catchability lower & upper bounds
init_number qlbound;
init_number qubound;
// Catchability Initial Guesses
init_vector qinit(1,nsurvey);
// Survey Options - 1
//
// column 1 - Estimated Selectivity - Initial Value
// column 2 - Estimated Selectivity - Lower Bound
// column 3 - Estimated Selectivity - Upper Bound
// column 4 - Estimated Selectivity - Phase
init_matrix survey_opt1(1,nsurvey,1,4);
// Survey Options - 2
//
// column 1 - CPUE Exponent - Initial Value
// column 2 - CPUE Exponent - Lower Bound
// column 3 - CPUE Exponent - Upper Bound
// column 4 - CPUE Exponent - Phase
init_matrix survey_opt2(1,nsurvey,1,4);
// Weights per fish
//
// Column 1 - Landings
// Column 2 - Discards
// Column 3 - Recruits
// Column 4 - Post-Recruits
//
init_matrix wtmat(1,nyears,1,4);
// MCMC stuff
init_int doMCMC;
init_int MCMCnboot;
init_int MCMCnthin;
init_int MCMCseed;
LOCAL_CALCS
if (doMCMC == 1)
{
for (int ix=1; ix<= nyears; ix++)
{
basicMCMC << "F_" << nfyear+ix-1 << " ";
}
for (int ix=1; ix<= nyears; ix++)
{
basicMCMC << "RHAT_" << nfyear+ix-1 << " ";
}
for (int ix=1; ix<= nyears; ix++)
{
basicMCMC << "PHAT_" << nfyear+ix-1 << " ";
}
for (int jx=1; jx<= nsurvey; jx++)
{
basicMCMC << "QHAT_" << jx << " ";
}
for (int jx=1; jx<= nsurvey; jx++)
{
basicMCMC << "SRX_" << jx << " ";
}
for (int jx=1; jx<= nsurvey; jx++)
{
basicMCMC << "Beta_" << jx << " ";
}
basicMCMC << endl;
}
END_CALCS
vector obs_land(1,nyears);
vector obs_disc(1,nyears);
vector obs_catch(1,nyears);
vector obs_survey(1,nobs);
vector cv(1,nobs);
vector ts(1,nobs);
vector sr(1,nobs);
ivector kyear(1,nobs);
ivector ksurv(1,nobs);
ivector ktype(1,nobs);
vector wt_land(1,nyears);
vector wt_disc(1,nyears);
vector wt_r(1,nyears);
vector wt_pr(1,nyears);
vector srx_init(1,nsurvey);
vector srx_lb(1,nsurvey);
vector srx_ub(1,nsurvey);
ivector srx_phase(1,nsurvey);
vector beta_init(1,nsurvey);
vector beta_lb(1,nsurvey);
vector beta_ub(1,nsurvey);
ivector beta_phase(1,nsurvey);
int ks;
int ny;
int ixx;
LOCAL_CALCS
struct tm *Today;
time_t xstart;
time(&xstart);
Today = localtime(&xstart);
strftime(dtstring,12,"%d_%b_%Y",Today);
strftime(tmstring,6,"%H:%M",Today);
srx_init = column(survey_opt1,1);
srx_lb = column(survey_opt1,2);
srx_ub = column(survey_opt1,3);
beta_init = column(survey_opt2,1);
beta_lb = column(survey_opt2,2);
beta_ub = column(survey_opt2,3);
for (ixx = 1; ixx <= nsurvey; ixx++)
{
srx_phase[ixx] = (int) survey_opt1[ixx][4];
beta_phase[ixx] = (int) survey_opt2[ixx][4];
}
END_CALCS
PARAMETER_SECTION
// All intrinsically positive parameters estimated as logs (LDJ, 7/18/2013).
// Read in log scale parameter values and make space for their transformed
// values in caclulations. For each parameter, determine the specified bounds,
// declare the log parameter and make space for the arithmetic value that will
// be set in the PROCEDURE_SECTION. Initial values are set in the
// PRELIMINARY_CALCS section. LDJ, 19 July 2013
// Indirect declaration for bounds seems to work best here
!! double lob=1e-9;double upb=1e9;
!! lob=log(lob);upb=log(upb);
init_bounded_vector logrhat(1,nyears,lob,upb);
vector rhat(1,nyears);
!! lob=1e-9;upb=3;
!! lob=log(lob);upb=log(upb);
init_bounded_vector logf_calc(1,nyears,lob,upb);
vector f_calc(1,nyears);
!! lob=1e-9;upb=1e3;
!! lob=log(lob);upb=log(upb);
init_bounded_vector logqhat(1,nsurvey,lob,upb);
vector qhat(1,nsurvey);
!! lob=1e-9;upb=1e9;
!! lob=log(lob);upb=log(upb);
init_bounded_number logphat1(lob,upb);
//indirect declaration that seems to work (LDJ, 19 July 2013)
!!dvector vlob(1,nsurvey);
!!dvector vupb(1,nsurvey);
!!vlob=log(srx_lb);vupb=log(srx_ub);
init_bounded_number_vector logsrx(1,nsurvey,vlob,vupb,srx_phase);
vector srx(1,nsurvey);
!!EO("in parameter section");EO(logsrx);
!!vlob=log(beta_lb);vupb=log(beta_ub);
init_bounded_number_vector logbeta(1,nsurvey,vlob,vupb,beta_phase);
vector beta(1,nsurvey);
!!EO(logbeta);
number phat1;
sdreport_vector phat(1,nyears+1);
sdreport_vector totnum(1,nyears);
sdreport_vector totbio(1,nyears);
number c;
number m;
number n;
number z;
number f;
number cx;
number nx;
number lx;
number lxc;
number lxr;
number lxp;
number lxt;
number lxw;
number xxx;
vector z_calc(1,nyears);
vector c_calc(1,nyears);
vector pred_survey(1,nobs);
vector resid_log(1,nobs);
vector resid_catch(1,nyears);
vector l_calc(1,nobs);
vector l_calc_wt(1,nobs);
vector l_resid(1,nobs);
vector survey_negloglike(1,nsurvey);
vector catch_negloglike(1,nyears);
vector catch_likeresid(1,nyears);
vector bio_r(1,nyears);
vector bio_p(1,nyears);
vector survey_r_like(1,nsurvey);
vector survey_p_like(1,nsurvey);
number p_est;
number r_est;
number t_est;
number w_est;
number c_est;
number d_est;
number sd;
objective_function_value negloglikelihood;
PRELIMINARY_CALCS_SECTION
obs_land = column(total_catch,1);
obs_disc = column(total_catch,2);
for (int ix=1; ix <= nobs; ix++)
{
ksurv[ix] = (double) survey_data[ix][1];
kyear[ix] = (double) survey_data[ix][2];
ktype[ix] = (double) survey_data[ix][3];
}
obs_survey = column(survey_data,4);
cv = column(survey_data,5);
ts = column(survey_data,6);
sr = column(survey_data,7);
wt_land = column(wtmat,1);
wt_disc = column(wtmat,2);
wt_r = column(wtmat,3);
wt_pr = column(wtmat,4);
// Very last thing-supply initial log scale parameter estimates (LDJ)
for (int i=1;i<=nyears;i++) logf_calc(i)=log(0.5);
logqhat=log(qinit);
for (int i=1;i<=nyears;i++) logrhat(i)=log(0.0001);
for (int ix=1; ix <= nobs; ix++)
{
if (ktype[ix] == 1)
{
ny = kyear[ix] - nfyear+1;
ks = ksurv[ix];
logrhat[ny] = log(obs_survey[ix]) - log(sr[ix]) - logqhat[ks];
}
}
for (int ix=1; ix <= nobs; ix++)
{
ks = ksurv[ix];
if (ktype[ix] == 2 && kyear[ix] == nfyear)
{
phat1 = obs_survey[ix] / value(mfexp(logqhat[ks]));
}
}
logphat1 = log(phat1);
// PROBLEM HERE-cannot assign to an init_bounded_number_vector
logsrx=log(srx_init);
logbeta=log(beta_init);
PROCEDURE_SECTION
int i; // loop counter (LDJ)
negloglikelihood=0.0;
lxc = 0.0;
lxr = 0.0;
lxp = 0.0;
lxt = 0.0;
lxw = 0.0;
// convert parameters for intrinsically positive parameters
// to their nautral arithmetic values (LDJ, 19 July 2013)
rhat(1,nyears)=mfexp(logrhat);
f_calc=mfexp(logf_calc);
qhat=mfexp(logqhat);
phat1=mfexp(logphat1);
EO("procedure section");
for (int i=1; i<= nsurvey; i++)
{
dvariable temp=logsrx(i);
srx(i)=mfexp(temp);
}
EO(logsrx);EO(srx);
beta=mfexp(logbeta);
EO(logbeta);EO(beta);
phat[1] = mfexp(logphat1);
// for (double i = 1; i <= nyears; i++)
for (int i = 1; i <= nyears; i++)
{
obs_catch[i] = obs_land[i] + obs_disc[i];
m = nmort[i];
n = phat[i] + rhat[i];
f = f_calc[i];
calc_catch();
z_calc[i] = z;
c_calc[i] = cx;
phat[i+1] = nx;
totnum[i] = n;
}
sd = sqrt(log(1.0 + square(cvcatch)));
// for (double i=1; i <= nyears; i++)
for (int i=1; i <= nyears; i++)
{
if (obs_catch[i] > 0.0 && c_calc[i] > 0.0)
{
resid_catch[i] = log(obs_catch[i]) - log(c_calc[i]);
lx = log(sd) + 0.5 * square(resid_catch[i] / sd);
catch_negloglike[i] = lx;
catch_likeresid[i] = lx - log(sd);
negloglikelihood += lx * lambda_catch;
lxc += lx * lambda_catch;
}
}
// for (double i=1; i <= nsurvey; i++)
for (int i=1; i <= nsurvey; i++)
{
survey_negloglike[i] = 0.0;
survey_r_like[i] = 0.0;
survey_p_like[i] = 0.0;
}
// for (double j=1; j <= nobs; j++)
for (int j=1; j <= nobs; j++)
{
ks = ksurv[j];
ny = kyear[j] - nfyear+1;
if (ktype[j] == 1)
pred_survey[j] = qhat[ks] * rhat[ny] * sr[j] * exp(-ts[j]*z_calc[ny]);
else if (ktype[j] == 2)
pred_survey[j] = qhat[ks] * phat[ny] * exp(-ts[j]*z_calc[ny]);
else if (ktype[j] == 3)
pred_survey[j] = qhat[ks] * exp(-ts[j]*z_calc[ny]) * (rhat[ny] * srx[ks] + phat[ny]);
else
{
xxx = (rhat[ny] * srx[ks] + phat[ny]) * exp(-ts[j]*z_calc[ny]);
pred_survey[j] = qhat[ks] * pow(xxx,beta[ks]);
}
sd = sqrt(log(1.0 + square(cv[j])));
if (obs_survey[j] > 0.0 && pred_survey[j] > 0.0)
resid_log[j] = log(obs_survey[j]) - log(pred_survey[j]);
else
resid_log[j] = 0.0;
l_calc[j] = log(sd) + 0.5 * square(resid_log[j] / sd);
l_calc_wt[j] = l_calc[j] * lambda_survey[ks];
l_resid[j] = l_calc[j] - log(sd);
if (ktype[j] == 1)
survey_r_like[ks] += l_resid[j];
else if (ktype[j] == 2)
survey_p_like[ks] += l_resid[j];
survey_negloglike[ks] += l_calc_wt[j];
negloglikelihood += l_calc_wt[j];
if (ktype[j] == 1)
lxr += l_calc_wt[j];
else if (ktype[j] == 2)
lxp += l_calc_wt[j];
else if (ktype[j] == 3 || ktype[j] == 4)
lxt += l_calc_wt[j];
else
lxw += l_calc_wt[j];
}
// ADMB (Microsoft compiler) does not like the way double
// and int variables are mixed here (LDJ, July 19, 2013)
// for (double i = 1; i <=nyears; i++)
for(i = 1; i <=nyears; i++)
{
bio_r[i] = wt_r[i] * rhat[i];
bio_p[i] = wt_pr[i] * phat[i];
totbio[i] = bio_r[i] + bio_p[i];
}
// cout << last_phase() << endl;
if (mceval_phase())
{
write_MCMC();
}
EO(rhat);
EO(f_calc);
EO(qhat);
EO(phat1);
EO(srx);
EO(beta);
EO(rhat);
EO(phat)
EO(negloglikelihood);
exit(1);
FUNCTION calc_catch
z = f + m;
cx = f * n * (1.0 - exp(-z)) / z;
nx = n * exp(-z);
FUNCTION write_MCMC
int i;
for (i= 1; i <= nyears; i++)
{
basicMCMC << f_calc[i] << " ";
}
for (i= 1; i <= nyears; i++)
{
basicMCMC << rhat[i] << " ";
}
for (i= 1; i <= nyears; i++)
{
basicMCMC << phat[i] << " ";
}
for (i= 1; i <= nsurvey; i++)
{
basicMCMC << qhat[i] << " ";
}
for (i= 1; i <= nsurvey; i++)
{
basicMCMC << srx[i] << " ";
}
for (i= 1; i <= nsurvey; i++)
{
basicMCMC << beta[i] << " ";
}
basicMCMC << endl;
REPORT_SECTION
int i;
report << "CSA Version 4.1 - 18 July 2013" << endl << endl;
report << "Date_of_Run: " << dtstring << endl;
report << "Time_of_Run: " << tmstring << endl << endl;
report << "First Year in Data = " << nfyear << endl;
report << "Last Year in Data = " << nxyear << endl << endl;
report << "Number of Years = " << nyears << endl << endl;
report << "Number of Surveys = " << nsurvey << endl << endl;
report << "Number of Observations = " << nobs << endl << endl;
report << "Negative Log-Likelihood = " << negloglikelihood << endl << endl;
report << " - Catch = " << lxc << endl;
report << " - Recruits = " << lxr << endl;
report << " - Post-Recruits = " << lxp << endl;
report << " - Aggregate Surveys = " << lxt << endl;
report << " - CPUE = " << lxw << endl << endl;
report << " - Surveys " << endl;
for (i = 1; i <= nsurvey; i++)
report << " - Survey - " << i << " = " << survey_negloglike[i] << endl;
report << endl << endl;
report << "Log-Likelihood Residuals " << endl;
report << " Recruits " << endl;
for (i = 1; i <= nsurvey; i++)
report << " - Survey - " << i << " = " << survey_r_like[i] << endl;
report << endl;
report << " Post-Recruits " << endl;
for (i = 1; i <= nsurvey; i++)
report << " - Survey - " << i << " = " << survey_p_like[i] << endl;
report << endl << endl;
report << "Catchabilty" << endl;
for (i = 1; i <= nsurvey; i++)
report << "Survey - " << i << " = " << qhat[i] << endl;
report << endl;
report << "Estimable Survey Selectivity" << endl;
for (i = 1; i <= nsurvey; i++)
report << "Survey - " << i << " = " << srx[i] << endl;
report << endl;
report << "CPUE Exponent" << endl;
for (i = 1; i <= nsurvey; i++)
report << "Survey - " << i << " = " << beta[i] << endl;
report << endl;
report << scientific << setw(18) << setprecision(8);
report << endl;
report << "Catch Observed Predicted Log-Resid neg-log_likehood like_resid" << endl;
// ADMB (Microsoft compiler) does not like the way double
// and int variables are mixed here (LDJ, July 19, 2013)
// for (i = 1; i <=nyears; i++)
for(i = 1; i <=nyears; i++)
{
report << (int) nfyear+i-1 << " " << obs_catch[i] << " " << c_calc[i] << " " << resid_catch[i] << " " << catch_negloglike[i] << " " << catch_likeresid[i] << endl;
}
report << endl;
report << "Survey Year Type Observed Predicted Log-Resid neg-log_likehood like-resid likehood_wt" << endl << endl;
for (i=1; i <= nobs; i++)
{
report << ksurv[i] << " " << kyear[i] << " " << ktype[i] << " ";
report << obs_survey[i] << " " << pred_survey[i] << " ";
report << resid_log[i] << " " << l_calc[i] << " ";
report << l_resid[i] << " " << l_calc_wt[i] << endl;
}
report << endl;
report << fixed << setw(15) << setprecision(6);
report << endl;
report << "Mortality M F Z" << endl;
for (i = 1; i <=nyears; i++)
{
report << nfyear+i-1 << " " << nmort[i] << " " << f_calc[i] << " " << z_calc[i] << endl;
}
report << endl;
report << scientific << setw(18) << setprecision(8);
report << endl;
report << "Population Estimate Recruits Post-Recruits Total" << endl;
for (i = 1; i <=nyears; i++)
{
r_est = rhat[i];
p_est = phat[i];
t_est = r_est + p_est;
report << nfyear+i-1 << " " << r_est << " " << p_est << " " << t_est << endl;
}
report << endl;
report << "Catch Biomass Estimate Landings Discards Total" << endl;
for (i = 1; i <=nyears; i++)
{
c_est = c_calc[i] * wt_land[i];
d_est = obs_disc[i] * wt_disc[i];
w_est = c_est + d_est;
report << nfyear+i-1 << " " << c_est << " " << d_est << " " << w_est << endl;
}
report << endl;
report << "Stock Biomass Estimate Recruits Post-Recruits Total" << endl;
for (i = 1; i <=nyears; i++)
{
report << nfyear+i-1 << " " << bio_r[i] << " " << bio_p[i] << " " << totbio[i] << endl;
}
report << endl;
// END of MODEL
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Lucky1.DAT
Type: application/ms-tnef
Size: 6772 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/users/attachments/20130719/30f3a6eb/attachment.bin>
More information about the Users
mailing list