diff --git a/src/dcrws_ssm_argos.hpp b/src/dcrws_ssm_argos.hpp new file mode 100644 index 0000000000000000000000000000000000000000..8d9cbf30a2710bdc6b2f0507b28bfb01dc4079ec --- /dev/null +++ b/src/dcrws_ssm_argos.hpp @@ -0,0 +1,136 @@ +#undef TMB_OBJECTIVE_PTR +#define TMB_OBJECTIVE_PTR obj + +using namespace density; + + +template +Type dcrwSSMargos(objective_function * obj) { + + // data + + // names of each track + vector tracknames = StringList("tracknames",obj); + int ntracks = tracknames.size(); + vector > tracks(ntracks); // vector< kind of element> name(size) + + // tracks + for(int i = 0; i < ntracks; i++){ + track indtrack(tracknames[i],obj); + tracks[i] = indtrack; + } + + + + REPORT(ntracks); + + + // parameters + + // Input parameters - i.e. parameters to estimate. + PARAMETER_VECTOR(working_theta); + vector theta(2); + theta(0) = 2*M_PI/(1.0+exp(-working_theta(0))) - M_PI; // -pi to pi + theta(1) = 2*M_PI/(1.0+exp(-working_theta(1))); // 0 to 2pi + ADREPORT(theta); + + PARAMETER_VECTOR(working_gamma); + vector gamma(2); + gamma(0) = 1.0/(1.0+exp(-working_gamma(0))); + gamma(1) = gamma(0)/(1.0+exp(-working_gamma(1))); + ADREPORT(gamma); + + PARAMETER(working_sigma_lon); + Type sigma_lon = exp(working_sigma_lon); + ADREPORT(sigma_lon); + PARAMETER(working_sigma_lat); + Type sigma_lat = exp(working_sigma_lat); + ADREPORT(sigma_lat); + // Variance-covariance matrix for the process equation. + matrix Sigma(2,2); + Sigma << sigma_lon*sigma_lon, 0.0, + 0.0, sigma_lat*sigma_lat; + + PARAMETER_MATRIX(working_A); //matrix of switching probabilities + int m = working_A.rows(); // number of states + matrix A(m, m); + for(int i = 0; i < m; ++i){ + for(int j = 0; j < (m-1); ++j){ + A(i,j) = exp(working_A(i,j)); //ensures the probabilities we estimate are >0 + } + A(i,m-1) = 1.0; + A.row(i) *= 1.0/A.row(i).sum(); + } + ADREPORT(A); + // system of equations stat dist + matrix system(m, m); + for(int i=0; i ones(m); + for(int i=0; i sys_inverse = system.inverse(); + vector delta = sys_inverse*ones; + + PARAMETER(working_psi); //working value for measurement error + Type psi = exp(working_psi); + ADREPORT(psi); + + + + + + // likelihood calculations + + //Calculate logLikelihood + Type nll = 0.0; + + // behaviour equation + vector behav(ntracks); + Type tmpb = 0.0; + for(int i=0; i < ntracks; ++i){ + tmpb = behaviour(tracks[i].b, delta, A); + behav(i) = tmpb; + } + REPORT(behav); + nll -= behav.sum(); + + // movement equation + Type tmpp = 0.0; + vector move(ntracks); + for(int i=0; i < ntracks; ++i){ + tmpp = movement(tracks[i].x, tracks[i].b, gamma, theta, Sigma); + move(i) = tmpp; + } + REPORT(move); + nll += move.sum(); + + // measurement equation + Type tmpm = 0.0; + vector meas(ntracks); + for(int i=0; i < ntracks; ++i){ + tmpm = measurement_argos(tracks[i].y, tracks[i].x, tracks[i].idx, tracks[i].jidx, tracks[i].ae, psi); + meas(i) = tmpm; + } + REPORT(meas); + nll += meas.sum(); + + for(int i=0; i -Type dcrwSSM(objective_function * obj) { +Type dcrwSSMcombo(objective_function * obj) { // data @@ -82,15 +82,15 @@ Type dcrwSSM(objective_function * obj) { matrix sys_inverse = system.inverse(); vector delta = sys_inverse*ones; - // PARAMETER(working_psi); //working value for measurement error - // Type psi = exp(working_psi); - // ADREPORT(psi); - // - // PARAMETER_VECTOR(working_gps_err); // working value for gps measurement error - // vector gps_err = exp(working_gps_err); - // matrix gpsSigma(2,2); - // gpsSigma << gps_err(0)*gps_err(0), 0.0, - // 0.0, gps_err(1)*gps_err(1); + PARAMETER(working_psi); //working value for measurement error + Type psi = exp(working_psi); + ADREPORT(psi); + + PARAMETER_VECTOR(working_gps_err); // working value for gps measurement error + vector gps_err = exp(working_gps_err); + matrix gpsSigma(2,2); + gpsSigma << gps_err(0)*gps_err(0), 0.0, + 0.0, gps_err(1)*gps_err(1); @@ -124,9 +124,13 @@ Type dcrwSSM(objective_function * obj) { // measurement equation Type tmpm = 0.0; vector meas(ntracks); + // for(int i=0; i < ntracks; ++i){ + // int measdatatype=tracks[i].datatype(0); + // tmpm = measurement(obj, tracks[i].y, tracks[i].x, tracks[i].kind, tracks[i].idx, tracks[i].jidx, tracks[i].ae, measdatatype); + // meas(i) = tmpm; + // } for(int i=0; i < ntracks; ++i){ - int measdatatype=tracks[i].datatype(0); - tmpm = measurement(obj, tracks[i].y, tracks[i].x, tracks[i].kind, tracks[i].idx, tracks[i].jidx, tracks[i].ae, measdatatype); + tmpm = measurement_combo(tracks[i].y, tracks[i].x, tracks[i].kind, tracks[i].idx, tracks[i].jidx, tracks[i].ae, psi, gpsSigma); meas(i) = tmpm; } // for(int i=0; i < ntracks; ++i){ diff --git a/src/dcrws_ssm_gps.hpp b/src/dcrws_ssm_gps.hpp new file mode 100644 index 0000000000000000000000000000000000000000..c95285c5abb4dbb27999090295c7ba713f2d0649 --- /dev/null +++ b/src/dcrws_ssm_gps.hpp @@ -0,0 +1,140 @@ +#undef TMB_OBJECTIVE_PTR +#define TMB_OBJECTIVE_PTR obj + +using namespace density; + + +template +Type dcrwSSMgps(objective_function * obj) { + + // data + + // names of each track + vector tracknames = StringList("tracknames",obj); + int ntracks = tracknames.size(); + vector > tracks(ntracks); // vector< kind of element> name(size) + + // tracks + for(int i = 0; i < ntracks; i++){ + track indtrack(tracknames[i],obj); + tracks[i] = indtrack; + } + + + + REPORT(ntracks); + + + // parameters + + // Input parameters - i.e. parameters to estimate. + PARAMETER_VECTOR(working_theta); + vector theta(2); + theta(0) = 2*M_PI/(1.0+exp(-working_theta(0))) - M_PI; // -pi to pi + theta(1) = 2*M_PI/(1.0+exp(-working_theta(1))); // 0 to 2pi + ADREPORT(theta); + + PARAMETER_VECTOR(working_gamma); + vector gamma(2); + gamma(0) = 1.0/(1.0+exp(-working_gamma(0))); + gamma(1) = gamma(0)/(1.0+exp(-working_gamma(1))); + ADREPORT(gamma); + + PARAMETER(working_sigma_lon); + Type sigma_lon = exp(working_sigma_lon); + ADREPORT(sigma_lon); + PARAMETER(working_sigma_lat); + Type sigma_lat = exp(working_sigma_lat); + ADREPORT(sigma_lat); + // Variance-covariance matrix for the process equation. + matrix Sigma(2,2); + Sigma << sigma_lon*sigma_lon, 0.0, + 0.0, sigma_lat*sigma_lat; + + PARAMETER_MATRIX(working_A); //matrix of switching probabilities + int m = working_A.rows(); // number of states + matrix A(m, m); + for(int i = 0; i < m; ++i){ + for(int j = 0; j < (m-1); ++j){ + A(i,j) = exp(working_A(i,j)); //ensures the probabilities we estimate are >0 + } + A(i,m-1) = 1.0; + A.row(i) *= 1.0/A.row(i).sum(); + } + ADREPORT(A); + // system of equations stat dist + matrix system(m, m); + for(int i=0; i ones(m); + for(int i=0; i sys_inverse = system.inverse(); + vector delta = sys_inverse*ones; + + + PARAMETER_VECTOR(working_gps_err); // working value for gps measurement error + vector gps_err = exp(working_gps_err); + matrix gpsSigma(2,2); + gpsSigma << gps_err(0)*gps_err(0), 0.0, + 0.0, gps_err(1)*gps_err(1); + + + + + + // likelihood calculations + + //Calculate logLikelihood + Type nll = 0.0; + + // behaviour equation + vector behav(ntracks); + Type tmpb = 0.0; + for(int i=0; i < ntracks; ++i){ + tmpb = behaviour(tracks[i].b, delta, A); + behav(i) = tmpb; + } + REPORT(behav); + nll -= behav.sum(); + + // movement equation + Type tmpp = 0.0; + vector move(ntracks); + for(int i=0; i < ntracks; ++i){ + tmpp = movement(tracks[i].x, tracks[i].b, gamma, theta, Sigma); + move(i) = tmpp; + } + REPORT(move); + nll += move.sum(); + + // measurement equation + Type tmpm = 0.0; + vector meas(ntracks); + for(int i=0; i < ntracks; ++i){ + tmpm = measurement_gps(tracks[i].y, tracks[i].x, tracks[i].idx, tracks[i].jidx, tracks[i].ae, gpsSigma); + meas(i) = tmpm; + } + + REPORT(meas); + nll += meas.sum(); + + for(int i=0; i x, vector b, vector gamma, vector th return rslt; } -// // measurement process -// template -// Type measurement(array y, matrix x, vector idx, vector jidx, matrix ae, Type psi){ -// -// vector measnll(y.cols()); -// vector xhat(2); -// vector tmp(2); -// vector tmp2(2); -// -// for(int i = 0; i < y.cols(); ++i){ -// xhat = (1.0-jidx(i)) * x.col(idx(i)-1) + jidx(i)*x.col(idx(i)); -// tmp = y.col(i); -// tmp2 = tmp-xhat; -// // nll += nll_meas(tmp-xhat); //MVNORM_t returns the negative log likelihood -// // taken from Marie's DCRW_Argos.cpp file -// // dt returns the log likelihood -// measnll(i) = -(0.5*log(psi) - log(ae(i,0)) + dt(sqrt(psi)*tmp2(0)/ae(i,0),ae(i,1),true) ); // Longitude -// measnll(i) -= (0.5*log(psi) - log(ae(i,2)) + dt(sqrt(psi)*tmp2(1)/ae(i,2),ae(i,3),true) ); // Latitude -// } -// -// Type rslt = measnll.sum(); -// return rslt; // negative log-lik -// } +// measurement process +template +Type measurement_argos(array y, matrix x, vector idx, vector jidx, matrix ae, Type psi){ + + vector measnll(y.cols()); + vector xhat(2); + vector tmp(2); + vector tmp2(2); + + for(int i = 0; i < y.cols(); ++i){ + xhat = (1.0-jidx(i)) * x.col(idx(i)-1) + jidx(i)*x.col(idx(i)); + tmp = y.col(i); + tmp2 = tmp-xhat; + // nll += nll_meas(tmp-xhat); //MVNORM_t returns the negative log likelihood + // taken from Marie's DCRW_Argos.cpp file + // dt returns the log likelihood + measnll(i) = -(0.5*log(psi) - log(ae(i,0)) + dt(sqrt(psi)*tmp2(0)/ae(i,0),ae(i,1),true) ); // Longitude + measnll(i) -= (0.5*log(psi) - log(ae(i,2)) + dt(sqrt(psi)*tmp2(1)/ae(i,2),ae(i,3),true) ); // Latitude + } + + Type rslt = measnll.sum(); + return rslt; // negative log-lik +} + +template +Type measurement_gps(array y, matrix x, vector idx, vector jidx, matrix ae, matrix gpsSigma){ + + vector measnll(y.cols()); + vector xhat(2); + vector tmp(2); + vector tmp2(2); + density::MVNORM_t gps_nll(gpsSigma); + + for(int i = 0; i < y.cols(); ++i){ + xhat = (1.0-jidx(i)) * x.col(idx(i)-1) + jidx(i)*x.col(idx(i)); + tmp = y.col(i); + tmp2 = tmp-xhat; + measnll(i) = gps_nll(tmp2); //MVNORM_t returns the negative log likelihood + } + + Type rslt = measnll.sum(); + return rslt; // negative log-lik +} + + +template +Type measurement_combo(array y, matrix x, vector kind, vector idx, vector jidx, matrix ae, Type psi, matrix gpsSigma){ + + vector measnll(y.cols()); + vector xhat(2); + vector tmp(2); + vector tmp2(2); + density::MVNORM_t gps_nll(gpsSigma); + + for(int i = 0; i < y.cols(); ++i){ + xhat = (1.0-jidx(i)) * x.col(idx(i)-1) + jidx(i)*x.col(idx(i)); + tmp = y.col(i); + tmp2 = tmp-xhat; + + if(kind(i) == 0){ + measnll(i) = gps_nll(tmp2); //MVNORM_t returns the negative log likelihood + } else if(kind(i) == 1) { + // taken from Marie's DCRW_Argos.cpp file + // dt returns the log likelihood + measnll(i) = -(0.5*log(psi) - log(ae(i,0)) + dt(sqrt(psi)*tmp2(0)/ae(i,0),ae(i,1),true) ); // Longitude + measnll(i) -= (0.5*log(psi) - log(ae(i,2)) + dt(sqrt(psi)*tmp2(1)/ae(i,2),ae(i,3),true) ); // Latitude + } + } + + Type rslt = measnll.sum(); + return rslt; // negative log-lik +} // diff --git a/src/swim.cpp b/src/swim.cpp index fad2e7752efbad1b81969c6240f7cfdcd83a3ae0..8df08068ffe93d1c36d12e1bc26b991dcf415d07 100644 --- a/src/swim.cpp +++ b/src/swim.cpp @@ -5,13 +5,17 @@ #include "likelihoods.hpp" #include "measlikelihood.hpp" #include "dcrws_hmm.hpp" -#include "dcrws_ssm.hpp" +#include "dcrws_ssm_combo.hpp" +#include "dcrws_ssm_argos.hpp" +#include "dcrws_ssm_gps.hpp" // map the data_integer to a character enum model_choices { dcrwhmm = 0, - dcrwssm = 1, + dcrwssmargos = 1, + dcrwssmgps = 2, + dcrwssmcombo = 3, }; template @@ -24,8 +28,14 @@ switch(model){ case dcrwhmm: return dcrwHMM(this); break; - case dcrwssm: - return dcrwSSM(this); + case dcrwssmargos: + return dcrwSSMargos(this); + break; + case dcrwssmgps: + return dcrwSSMgps(this); + break; + case dcrwssmcombo: + return dcrwSSMcombo(this); break; } } diff --git a/src/trackstruct.hpp b/src/trackstruct.hpp index dea40af5178a0ca3a810c301f86fb03cdd51a0a5..ff6cc4f75b7219ce4134442a14b49a0b905441f4 100644 --- a/src/trackstruct.hpp +++ b/src/trackstruct.hpp @@ -9,7 +9,7 @@ struct track { array y; vector b; vector kind; - vector datatype; + // vector datatype; vector idx; vector jidx; matrix ae; @@ -30,7 +30,7 @@ track::track(std::string name, objective_function *obj) { // objecti y = tmbutils::asArray(getListElement(tracking_data, "y")); b = asVector(getListElement(tracking_data, "b")); kind = asVector(getListElement(tracking_data, "kind")); - datatype = asVector(getListElement(tracking_data, "datatype")) ; + // datatype = asVector(getListElement(tracking_data, "datatype")) ; idx = asVector(getListElement(tracking_data, "idx")); jidx = asVector(getListElement(tracking_data, "jidx")); ae = asMatrix(getListElement(tracking_data, "ae"));