diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 9830dd1..fdf7ba1 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -564,16 +564,16 @@ BEGIN_RCPP END_RCPP } // EM_LBFGS_nhmm_singlechannel -Rcpp::List EM_LBFGS_nhmm_singlechannel(arma::mat& eta_pi, const arma::mat& X_pi, arma::cube& eta_A, const arma::cube& X_A, arma::cube& eta_B, const arma::cube& X_B, const arma::umat& obs, const arma::uvec& Ti, const bool icpt_only_pi, const bool icpt_only_A, const bool icpt_only_B, const bool iv_A, const bool iv_B, const bool tv_A, const bool tv_B, const arma::uword n_obs, const arma::uword maxeval, const double ftol_abs, const double ftol_rel, const double xtol_abs, const double xtol_rel, const arma::uword print_level, const arma::uword maxeval_m, const double ftol_abs_m, const double ftol_rel_m, const double xtol_abs_m, const double xtol_rel_m, const arma::uword print_level_m, const double lambda, const double pseudocount); +Rcpp::List EM_LBFGS_nhmm_singlechannel(const arma::mat& eta_pi, const arma::mat& X_pi, const arma::cube& eta_A, const arma::cube& X_A, const arma::cube& eta_B, const arma::cube& X_B, const arma::umat& obs, const arma::uvec& Ti, const bool icpt_only_pi, const bool icpt_only_A, const bool icpt_only_B, const bool iv_A, const bool iv_B, const bool tv_A, const bool tv_B, const arma::uword n_obs, const arma::uword maxeval, const double ftol_abs, const double ftol_rel, const double xtol_abs, const double xtol_rel, const arma::uword print_level, const arma::uword maxeval_m, const double ftol_abs_m, const double ftol_rel_m, const double xtol_abs_m, const double xtol_rel_m, const arma::uword print_level_m, const double lambda, const double pseudocount); RcppExport SEXP _seqHMM_EM_LBFGS_nhmm_singlechannel(SEXP eta_piSEXP, SEXP X_piSEXP, SEXP eta_ASEXP, SEXP X_ASEXP, SEXP eta_BSEXP, SEXP X_BSEXP, SEXP obsSEXP, SEXP TiSEXP, SEXP icpt_only_piSEXP, SEXP icpt_only_ASEXP, SEXP icpt_only_BSEXP, SEXP iv_ASEXP, SEXP iv_BSEXP, SEXP tv_ASEXP, SEXP tv_BSEXP, SEXP n_obsSEXP, SEXP maxevalSEXP, SEXP ftol_absSEXP, SEXP ftol_relSEXP, SEXP xtol_absSEXP, SEXP xtol_relSEXP, SEXP print_levelSEXP, SEXP maxeval_mSEXP, SEXP ftol_abs_mSEXP, SEXP ftol_rel_mSEXP, SEXP xtol_abs_mSEXP, SEXP xtol_rel_mSEXP, SEXP print_level_mSEXP, SEXP lambdaSEXP, SEXP pseudocountSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::mat& >::type eta_pi(eta_piSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type eta_pi(eta_piSEXP); Rcpp::traits::input_parameter< const arma::mat& >::type X_pi(X_piSEXP); - Rcpp::traits::input_parameter< arma::cube& >::type eta_A(eta_ASEXP); + Rcpp::traits::input_parameter< const arma::cube& >::type eta_A(eta_ASEXP); Rcpp::traits::input_parameter< const arma::cube& >::type X_A(X_ASEXP); - Rcpp::traits::input_parameter< arma::cube& >::type eta_B(eta_BSEXP); + Rcpp::traits::input_parameter< const arma::cube& >::type eta_B(eta_BSEXP); Rcpp::traits::input_parameter< const arma::cube& >::type X_B(X_BSEXP); Rcpp::traits::input_parameter< const arma::umat& >::type obs(obsSEXP); Rcpp::traits::input_parameter< const arma::uvec& >::type Ti(TiSEXP); @@ -604,16 +604,16 @@ BEGIN_RCPP END_RCPP } // EM_LBFGS_nhmm_multichannel -Rcpp::List EM_LBFGS_nhmm_multichannel(arma::mat& eta_pi, const arma::mat& X_pi, arma::cube& eta_A, const arma::cube& X_A, arma::field& eta_B, const arma::cube& X_B, const arma::ucube& obs, const arma::uvec& Ti, const bool icpt_only_pi, const bool icpt_only_A, const bool icpt_only_B, const bool iv_A, const bool iv_B, const bool tv_A, const bool tv_B, const arma::uword n_obs, const arma::uword maxeval, const double ftol_abs, const double ftol_rel, const double xtol_abs, const double xtol_rel, const arma::uword print_level, const arma::uword maxeval_m, const double ftol_abs_m, const double ftol_rel_m, const double xtol_abs_m, const double xtol_rel_m, const arma::uword print_level_m, const double lambda, const double pseudocount); +Rcpp::List EM_LBFGS_nhmm_multichannel(const arma::mat& eta_pi, const arma::mat& X_pi, const arma::cube& eta_A, const arma::cube& X_A, const arma::field& eta_B, const arma::cube& X_B, const arma::ucube& obs, const arma::uvec& Ti, const bool icpt_only_pi, const bool icpt_only_A, const bool icpt_only_B, const bool iv_A, const bool iv_B, const bool tv_A, const bool tv_B, const arma::uword n_obs, const arma::uword maxeval, const double ftol_abs, const double ftol_rel, const double xtol_abs, const double xtol_rel, const arma::uword print_level, const arma::uword maxeval_m, const double ftol_abs_m, const double ftol_rel_m, const double xtol_abs_m, const double xtol_rel_m, const arma::uword print_level_m, const double lambda, const double pseudocount); RcppExport SEXP _seqHMM_EM_LBFGS_nhmm_multichannel(SEXP eta_piSEXP, SEXP X_piSEXP, SEXP eta_ASEXP, SEXP X_ASEXP, SEXP eta_BSEXP, SEXP X_BSEXP, SEXP obsSEXP, SEXP TiSEXP, SEXP icpt_only_piSEXP, SEXP icpt_only_ASEXP, SEXP icpt_only_BSEXP, SEXP iv_ASEXP, SEXP iv_BSEXP, SEXP tv_ASEXP, SEXP tv_BSEXP, SEXP n_obsSEXP, SEXP maxevalSEXP, SEXP ftol_absSEXP, SEXP ftol_relSEXP, SEXP xtol_absSEXP, SEXP xtol_relSEXP, SEXP print_levelSEXP, SEXP maxeval_mSEXP, SEXP ftol_abs_mSEXP, SEXP ftol_rel_mSEXP, SEXP xtol_abs_mSEXP, SEXP xtol_rel_mSEXP, SEXP print_level_mSEXP, SEXP lambdaSEXP, SEXP pseudocountSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< arma::mat& >::type eta_pi(eta_piSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type eta_pi(eta_piSEXP); Rcpp::traits::input_parameter< const arma::mat& >::type X_pi(X_piSEXP); - Rcpp::traits::input_parameter< arma::cube& >::type eta_A(eta_ASEXP); + Rcpp::traits::input_parameter< const arma::cube& >::type eta_A(eta_ASEXP); Rcpp::traits::input_parameter< const arma::cube& >::type X_A(X_ASEXP); - Rcpp::traits::input_parameter< arma::field& >::type eta_B(eta_BSEXP); + Rcpp::traits::input_parameter< const arma::field& >::type eta_B(eta_BSEXP); Rcpp::traits::input_parameter< const arma::cube& >::type X_B(X_BSEXP); Rcpp::traits::input_parameter< const arma::ucube& >::type obs(obsSEXP); Rcpp::traits::input_parameter< const arma::uvec& >::type Ti(TiSEXP); diff --git a/src/mnhmm_EM.cpp b/src/mnhmm_EM.cpp index 38c2142..869e83a 100644 --- a/src/mnhmm_EM.cpp +++ b/src/mnhmm_EM.cpp @@ -32,7 +32,7 @@ double mnhmm_base::objective_omega(const arma::vec& x, arma::vec& grad) { double val = arma::dot(counts.rows(idx), log_omega.rows(idx)); if (!std::isfinite(val)) { if (!grad.is_empty()) { - grad.fill(maxval); + grad.zeros(); } return maxval; } @@ -43,7 +43,7 @@ double mnhmm_base::objective_omega(const arma::vec& x, arma::vec& grad) { diff.rows(idx) = counts(idx) - sum_eo * omega.rows(idx); grad -= arma::vectorise(tQd * diff * X_omega.col(i).t()); if (!grad.is_finite()) { - grad.fill(maxval); + grad.zeros(); return maxval; } } @@ -127,7 +127,7 @@ double mnhmm_base::objective_pi(const arma::vec& x, arma::vec& grad) { double val = arma::dot(counts.rows(idx), log_pi(current_d).rows(idx)); if (!std::isfinite(val)) { if (!grad.is_empty()) { - grad.fill(maxval); + grad.zeros(); } return maxval; } @@ -135,10 +135,10 @@ double mnhmm_base::objective_pi(const arma::vec& x, arma::vec& grad) { // Only update grad if it's non-empty (i.e., for gradient-based optimization) if (!grad.is_empty()) { diff.zeros(); - diff = counts.rows(idx) - sum_epi * pi(current_d).rows(idx); + diff.rows(idx) = counts.rows(idx) - sum_epi * pi(current_d).rows(idx); grad -= arma::vectorise(tQs * diff * X_pi.col(i).t()); if (!grad.is_finite()) { - grad.fill(maxval); + grad.zeros(); return maxval; } } @@ -241,7 +241,7 @@ double mnhmm_base::objective_A(const arma::vec& x, arma::vec& grad) { double val = arma::dot(counts.rows(idx), log_A1.rows(idx)); if (!std::isfinite(val)) { if (!grad.is_empty()) { - grad.fill(maxval); + grad.zeros(); } return maxval; } @@ -251,7 +251,7 @@ double mnhmm_base::objective_A(const arma::vec& x, arma::vec& grad) { diff.rows(idx) = counts.rows(idx) - sum_ea * A1.rows(idx); grad -= arma::vectorise(tQs * diff * X_A.slice(i).col(t).t()); if (!grad.is_finite()) { - grad.fill(maxval); + grad.zeros(); return maxval; } } @@ -365,7 +365,7 @@ double mnhmm_sc::objective_B(const arma::vec& x, arma::vec& grad) { double val = e_b * log_B1(obs(t, i)); if (!std::isfinite(val)) { if (!grad.is_empty()) { - grad.fill(maxval); + grad.zeros(); } return maxval; } @@ -374,7 +374,7 @@ double mnhmm_sc::objective_B(const arma::vec& x, arma::vec& grad) { grad -= arma::vectorise(tQm * e_b * (I.col(obs(t, i)) - B1) * X_B.slice(i).col(t).t()); if (!grad.is_finite()) { - grad.fill(maxval); + grad.zeros(); return maxval; } } @@ -493,7 +493,7 @@ double mnhmm_mc::objective_B(const arma::vec& x, arma::vec& grad) { double val = e_b * log_B1(obs(current_c, t, i)); if (!std::isfinite(val)) { if (!grad.is_empty()) { - grad.fill(maxval); + grad.zeros(); } return maxval; } @@ -502,7 +502,7 @@ double mnhmm_mc::objective_B(const arma::vec& x, arma::vec& grad) { grad -= arma::vectorise(tQm * e_b * (I.col(obs(current_c, t, i)) - B1) * X_B.slice(i).col(t).t()); if (!grad.is_finite()) { - grad.fill(maxval); + grad.zeros(); return maxval; } } diff --git a/src/mnhmm_base.h b/src/mnhmm_base.h index 2a5b085..db51c47 100644 --- a/src/mnhmm_base.h +++ b/src/mnhmm_base.h @@ -72,9 +72,9 @@ struct mnhmm_base { const bool iv_B_, const bool tv_A_, const bool tv_B_, - arma::mat& eta_omega_, - arma::field& eta_pi_, - arma::field& eta_A_, + const arma::mat& eta_omega_, + const arma::field& eta_pi_, + const arma::field& eta_A_, const arma::uword n_obs_ = 0, const double lambda_ = 0, double maxval_ = 1e6) diff --git a/src/mnhmm_mc.h b/src/mnhmm_mc.h index 286f46a..a4df4fa 100644 --- a/src/mnhmm_mc.h +++ b/src/mnhmm_mc.h @@ -39,10 +39,10 @@ struct mnhmm_mc : public mnhmm_base { const bool tv_A_, const bool tv_B_, const arma::ucube& obs_, - arma::mat& eta_omega_, - arma::field& eta_pi_, - arma::field& eta_A_, - arma::field& eta_B_, + const arma::mat& eta_omega_, + const arma::field& eta_pi_, + const arma::field& eta_A_, + const arma::field& eta_B_, const arma::uword n_obs_ = 0, const double lambda_ = 0) : mnhmm_base(S_, D_, X_d_, X_pi_, X_s_, X_o_, Ti_, icpt_only_omega_, diff --git a/src/mnhmm_sc.h b/src/mnhmm_sc.h index 989b631..d6f5607 100644 --- a/src/mnhmm_sc.h +++ b/src/mnhmm_sc.h @@ -36,10 +36,10 @@ struct mnhmm_sc : public mnhmm_base { const bool tv_A_, const bool tv_B_, const arma::umat& obs_, - arma::mat& eta_omega_, - arma::field& eta_pi_, - arma::field& eta_A_, - arma::field& eta_B_, + const arma::mat& eta_omega_, + const arma::field& eta_pi_, + const arma::field& eta_A_, + const arma::field& eta_B_, const arma::uword n_obs_ = 0, const double lambda_ = 0) : mnhmm_base(S_, D_, X_d_, X_pi_, X_s_, X_o_, Ti_, icpt_only_omega_, diff --git a/src/nhmm_EM.cpp b/src/nhmm_EM.cpp index c26ea8f..e4a9b8c 100644 --- a/src/nhmm_EM.cpp +++ b/src/nhmm_EM.cpp @@ -14,7 +14,6 @@ double nhmm_base::objective_pi(const arma::vec& x, arma::vec& grad) { mstep_iter++; double value = 0; - eta_pi = arma::mat(x.memptr(), S - 1, K_pi); gamma_pi = sum_to_zero(eta_pi, Qs); arma::mat tQs = Qs.t(); @@ -31,27 +30,18 @@ double nhmm_base::objective_pi(const arma::vec& x, arma::vec& grad) { double sum_epi = arma::accu(counts.rows(idx)); // this is != 1 if pseudocounts are used double val = arma::dot(counts.rows(idx), log_pi.rows(idx)); if (!std::isfinite(val)) { - if (!grad.is_empty()) { - grad.fill(maxval); - } - return maxval; + grad.zeros(); + return arma::datum::inf; } value -= val; - // Only update grad if it's non-empty (i.e., for gradient-based optimization) - if (!grad.is_empty()) { - diff.zeros(); - diff = counts.rows(idx) - sum_epi * pi.rows(idx); - grad -= arma::vectorise(tQs * diff * X_pi.col(i).t()); - if (!grad.is_finite()) { - grad.fill(maxval); - return maxval; - } - } + + diff.zeros(); + diff.rows(idx) = counts.rows(idx) - sum_epi * pi.rows(idx); + grad -= arma::vectorise(tQs * diff * X_pi.col(i).t()); + } } - if (!grad.is_empty()) { - grad += lambda * x; - } + grad += lambda * x; return value + 0.5 * lambda * std::pow(arma::norm(x, 2), 2); } void nhmm_base::mstep_pi(const double xtol_abs, const double ftol_abs, @@ -72,18 +62,12 @@ void nhmm_base::mstep_pi(const double xtol_abs, const double ftol_abs, auto objective_pi_wrapper = [](unsigned n, const double* x, double* grad, void* data) -> double { auto* self = static_cast(data); arma::vec x_vec(const_cast(x), n, false, true); - if (grad) { - arma::vec grad_vec(grad, n, false, true); - return self->objective_pi(x_vec, grad_vec); - } else { - arma::vec grad_dummy; - return self->objective_pi(x_vec, grad_dummy); - } + arma::vec grad_vec(grad, n, false, true); + return self->objective_pi(x_vec, grad_vec); }; - arma::vec x_pi = arma::vectorise(eta_pi); + arma::vec x_pi(eta_pi.memptr(), eta_pi.n_elem, false, true); nlopt_opt opt_pi = nlopt_create(NLOPT_LD_LBFGS, x_pi.n_elem); - //nlopt_opt opt_pi = nlopt_create(NLOPT_LN_SBPLX, x_pi.n_elem); nlopt_set_min_objective(opt_pi, objective_pi_wrapper, this); nlopt_set_xtol_abs1(opt_pi, xtol_abs); nlopt_set_ftol_abs(opt_pi, ftol_abs); @@ -102,7 +86,6 @@ void nhmm_base::mstep_pi(const double xtol_abs, const double ftol_abs, nlopt_destroy(opt_pi); return; } - eta_pi = arma::mat(x_pi.memptr(), S - 1, K_pi); nlopt_destroy(opt_pi); } @@ -138,27 +121,21 @@ double nhmm_base::objective_A(const arma::vec& x, arma::vec& grad) { } double val = arma::dot(counts.rows(idx), log_A1.rows(idx)); if (!std::isfinite(val)) { - if (!grad.is_empty()) { - grad.fill(maxval); - } - return maxval; + Rcpp::Rcout<<"nonfinite val "< double { auto* self = static_cast(data); arma::vec x_vec(const_cast(x), n, false, true); - if (grad) { - arma::vec grad_vec(grad, n, false, true); - return self->objective_A(x_vec, grad_vec); - } else { - arma::vec grad_dummy; - return self->objective_A(x_vec, grad_dummy); - } + arma::vec grad_vec(grad, n, false, true); + return self->objective_A(x_vec, grad_vec); }; - arma::vec x_A(eta_A.slice(0).n_elem); - nlopt_opt opt_A = nlopt_create(NLOPT_LD_LBFGS, x_A.n_elem); - //nlopt_opt opt_A = nlopt_create(NLOPT_LN_SBPLX, x_A.n_elem); + nlopt_opt opt_A = nlopt_create(NLOPT_LD_LBFGS, eta_A.slice(0).n_elem); nlopt_set_min_objective(opt_A, objective_A_wrapper, this); nlopt_set_xtol_abs1(opt_A, xtol_abs); nlopt_set_ftol_abs(opt_A, ftol_abs); @@ -209,7 +179,7 @@ void nhmm_base::mstep_A(const double ftol_abs, const double ftol_rel, for (arma::uword s = 0; s < S; s++) { current_s = s; - x_A = arma::vectorise(eta_A.slice(s)); + arma::vec x_A(eta_A.slice(s).memptr(), eta_A.slice(s).n_elem, false, true); mstep_iter = 0; return_code = nlopt_optimize(opt_A, x_A.memptr(), &minf); if (print_level > 2 && return_code > 0) { @@ -222,7 +192,6 @@ void nhmm_base::mstep_A(const double ftol_abs, const double ftol_rel, nlopt_destroy(opt_A); return; } - eta_A.slice(s) = arma::mat(x_A.memptr(), S - 1, K_A); } nlopt_destroy(opt_A); } @@ -256,26 +225,21 @@ double nhmm_sc::objective_B(const arma::vec& x, arma::vec& grad) { } double val = e_b * log_B1(obs(t, i)); if (!std::isfinite(val)) { - if (!grad.is_empty()) { - grad.fill(maxval); - } - return maxval; + grad.zeros(); + return arma::datum::inf; } value -= val; - if (!grad.is_empty()) { - grad -= arma::vectorise(tQm * e_b * (I.col(obs(t, i)) - B1) * X_B.slice(i).col(t).t()); - if (!grad.is_finite()) { - grad.fill(maxval); - return maxval; - } + grad -= arma::vectorise(tQm * e_b * (I.col(obs(t, i)) - B1) * X_B.slice(i).col(t).t()); + if (!grad.is_finite()) { + grad.zeros(); + return arma::datum::inf; } } } } } - if (!grad.is_empty()) { - grad += lambda * x; - } + grad += lambda * x; + return value + 0.5 * lambda * std::pow(arma::norm(x, 2), 2); } void nhmm_sc::mstep_B(const double ftol_abs, const double ftol_rel, @@ -306,17 +270,10 @@ void nhmm_sc::mstep_B(const double ftol_abs, const double ftol_rel, auto objective_B_wrapper = [](unsigned n, const double* x, double* grad, void* data) -> double { auto* self = static_cast(data); arma::vec x_vec(const_cast(x), n, false, true); - if (grad) { - arma::vec grad_vec(grad, n, false, true); - return self->objective_B(x_vec, grad_vec); - } else { - arma::vec grad_dummy; - return self->objective_B(x_vec, grad_dummy); - } + arma::vec grad_vec(grad, n, false, true); + return self->objective_B(x_vec, grad_vec); }; - arma::vec x_B(eta_B.slice(0).n_elem); - nlopt_opt opt_B = nlopt_create(NLOPT_LD_LBFGS, x_B.n_elem); - //nlopt_opt opt_B = nlopt_create(NLOPT_LN_SBPLX, x_B.n_elem); + nlopt_opt opt_B = nlopt_create(NLOPT_LD_LBFGS, eta_B.slice(0).n_elem); nlopt_set_min_objective(opt_B, objective_B_wrapper, this); nlopt_set_xtol_abs1(opt_B, xtol_abs); nlopt_set_ftol_abs(opt_B, ftol_abs); @@ -327,7 +284,7 @@ void nhmm_sc::mstep_B(const double ftol_abs, const double ftol_rel, int return_code; for (arma::uword s = 0; s < S; s++) { current_s = s; - x_B = arma::vectorise(eta_B.slice(s)); + arma::vec x_B(eta_B.slice(s).memptr(), eta_B.slice(s).n_elem, false, true); mstep_iter = 0; return_code = nlopt_optimize(opt_B, x_B.memptr(), &minf); if (print_level > 2 && return_code > 0) { @@ -340,7 +297,6 @@ void nhmm_sc::mstep_B(const double ftol_abs, const double ftol_rel, nlopt_destroy(opt_B); return; } - eta_B.slice(s) = arma::mat(x_B.memptr(), M - 1, K_B); } nlopt_destroy(opt_B); } @@ -375,27 +331,21 @@ double nhmm_mc::objective_B(const arma::vec& x, arma::vec& grad) { } double val = e_b * log_B1(obs(current_c, t, i)); if (!std::isfinite(val)) { - if (!grad.is_empty()) { - grad.fill(maxval); - } - return maxval; + grad.zeros(); + return arma::datum::inf; } value -= val; - if (!grad.is_empty()) { - grad -= arma::vectorise(tQm * e_b * (I.col(obs(current_c, t, i)) - B1) * - X_B.slice(i).col(t).t()); - if (!grad.is_finite()) { - grad.fill(maxval); - return maxval; - } + grad -= arma::vectorise(tQm * e_b * (I.col(obs(current_c, t, i)) - B1) * + X_B.slice(i).col(t).t()); + if (!grad.is_finite()) { + grad.zeros(); + return arma::datum::inf; } } } } } - if (!grad.is_empty()) { - grad += lambda * x; - } + grad += lambda * x; return value + 0.5 * lambda * std::pow(arma::norm(x, 2), 2); } @@ -429,19 +379,13 @@ void nhmm_mc::mstep_B(const double ftol_abs, const double ftol_rel, auto objective_B_wrapper = [](unsigned n, const double* x, double* grad, void* data) -> double { auto* self = static_cast(data); arma::vec x_vec(const_cast(x), n, false, true); - if (grad) { - arma::vec grad_vec(grad, n, false, true); - return self->objective_B(x_vec, grad_vec); - } else { - arma::vec grad_dummy; - return self->objective_B(x_vec, grad_dummy); - } + arma::vec grad_vec(grad, n, false, true); + return self->objective_B(x_vec, grad_vec); }; double minf; int return_code; for (arma::uword c = 0; c < C; c++) { - arma::vec x_B(eta_B(c).slice(0).n_elem); - nlopt_opt opt_B = nlopt_create(NLOPT_LD_LBFGS, x_B.n_elem); + nlopt_opt opt_B = nlopt_create(NLOPT_LD_LBFGS, eta_B(c).slice(0).n_elem); nlopt_set_min_objective(opt_B, objective_B_wrapper, this); nlopt_set_xtol_abs1(opt_B, xtol_abs); nlopt_set_ftol_abs(opt_B, ftol_abs); @@ -451,7 +395,7 @@ void nhmm_mc::mstep_B(const double ftol_abs, const double ftol_rel, current_c = c; for (arma::uword s = 0; s < S; s++) { current_s = s; - x_B = arma::vectorise(eta_B(c).slice(s)); + arma::vec x_B(eta_B(c).slice(s).memptr(), eta_B(c).slice(s).n_elem, false, true); mstep_iter = 0; return_code = nlopt_optimize(opt_B, x_B.memptr(), &minf); if (print_level > 2 && return_code > 0) { @@ -471,9 +415,9 @@ void nhmm_mc::mstep_B(const double ftol_abs, const double ftol_rel, } // [[Rcpp::export]] Rcpp::List EM_LBFGS_nhmm_singlechannel( - arma::mat& eta_pi, const arma::mat& X_pi, - arma::cube& eta_A, const arma::cube& X_A, - arma::cube& eta_B, const arma::cube& X_B, + const arma::mat& eta_pi, const arma::mat& X_pi, + const arma::cube& eta_A, const arma::cube& X_A, + const arma::cube& eta_B, const arma::cube& X_B, const arma::umat& obs, const arma::uvec& Ti, const bool icpt_only_pi, const bool icpt_only_A, const bool icpt_only_B, const bool iv_A, const bool iv_B, const bool tv_A, const bool tv_B, @@ -678,9 +622,9 @@ Rcpp::List EM_LBFGS_nhmm_singlechannel( // [[Rcpp::export]] Rcpp::List EM_LBFGS_nhmm_multichannel( - arma::mat& eta_pi, const arma::mat& X_pi, - arma::cube& eta_A, const arma::cube& X_A, - arma::field& eta_B, const arma::cube& X_B, + const arma::mat& eta_pi, const arma::mat& X_pi, + const arma::cube& eta_A, const arma::cube& X_A, + const arma::field& eta_B, const arma::cube& X_B, const arma::ucube& obs, const arma::uvec& Ti, const bool icpt_only_pi, const bool icpt_only_A, const bool icpt_only_B, const bool iv_A, const bool iv_B, const bool tv_A, const bool tv_B, @@ -890,3 +834,4 @@ Rcpp::List EM_LBFGS_nhmm_multichannel( Rcpp::Named("relative_x_change") = relative_x_change ); } + diff --git a/src/nhmm_base.h b/src/nhmm_base.h index de8eacc..d37465f 100644 --- a/src/nhmm_base.h +++ b/src/nhmm_base.h @@ -59,8 +59,8 @@ struct nhmm_base { const bool iv_B_, const bool tv_A_, const bool tv_B_, - arma::mat& eta_pi_, - arma::cube& eta_A_, + const arma::mat& eta_pi_, + const arma::cube& eta_A_, const arma::uword n_obs_ = 0, const double lambda_ = 0, double maxval_ = 1e6) diff --git a/src/nhmm_mc.h b/src/nhmm_mc.h index bf79964..caf1277 100644 --- a/src/nhmm_mc.h +++ b/src/nhmm_mc.h @@ -36,9 +36,9 @@ struct nhmm_mc : public nhmm_base { const bool tv_A_, const bool tv_B_, const arma::ucube& obs_, - arma::mat& eta_pi_, - arma::cube& eta_A_, - arma::field& eta_B_, + const arma::mat& eta_pi_, + const arma::cube& eta_A_, + const arma::field& eta_B_, const arma::uword n_obs_ = 0, const double lambda_ = 0) : nhmm_base(S_, X_pi_, X_s_, X_o_, Ti_, icpt_only_pi_, icpt_only_A_, diff --git a/src/nhmm_sc.h b/src/nhmm_sc.h index e14f277..8c2cd77 100644 --- a/src/nhmm_sc.h +++ b/src/nhmm_sc.h @@ -33,9 +33,9 @@ struct nhmm_sc : public nhmm_base { const bool tv_A_, const bool tv_B_, const arma::umat& obs_, - arma::mat& eta_pi_, - arma::cube& eta_A_, - arma::cube& eta_B_, + const arma::mat& eta_pi_, + const arma::cube& eta_A_, + const arma::cube& eta_B_, const arma::uword n_obs_ = 0, const double lambda_ = 0) : nhmm_base(S_, X_pi_, X_s_, X_o_, Ti_, icpt_only_pi_, icpt_only_A_,