diff --git a/src/nlopt_wrapper.cpp b/src/nlopt_wrapper.cpp index f30caae5..5d00a237 100644 --- a/src/nlopt_wrapper.cpp +++ b/src/nlopt_wrapper.cpp @@ -50,11 +50,18 @@ std::unique_ptr new_nlopt_optimizer(const Rcpp::List throw Rcpp::exception("nlopt_create"); } + if(config.containsElementNamed("xtol_abs")) { + if(nlopt_set_xtol_abs1(opt.get(), Rcpp::as(config["xtol_abs"])) != NLOPT_SUCCESS) { + throw Rcpp::exception("nlopt_set_xtol_abs1"); + } + } + if(config.containsElementNamed("xtol_rel")) { if(nlopt_set_xtol_rel(opt.get(), Rcpp::as(config["xtol_rel"])) != NLOPT_SUCCESS) { throw Rcpp::exception("nlopt_set_xtol_rel"); } } + if(config.containsElementNamed("ftol_abs")) { if(nlopt_set_ftol_abs(opt.get(), Rcpp::as(config["ftol_abs"])) != NLOPT_SUCCESS) { throw Rcpp::exception("nlopt_set_ftol_abs"); @@ -93,20 +100,20 @@ std::unique_ptr new_nlopt_optimizer(const Rcpp::List // throw Rcpp::exception("nlopt_set_x_weights"); // } // } - -void set_uniform_xtol_abs(NloptStruct * opt, double value) { - if(nlopt_set_xtol_abs1(opt, value) != NLOPT_SUCCESS) { - throw Rcpp::exception("nlopt_set_xtol_abs1"); - } -} -void set_per_value_xtol_abs(NloptStruct * opt, const std::vector & xtol_abs) { - if(xtol_abs.size() != nlopt_get_dimension(opt)) { - throw Rcpp::exception("set_per_value_xtol_abs: parameter size mismatch"); - } - if(nlopt_set_xtol_abs(opt, xtol_abs.data()) != NLOPT_SUCCESS) { - throw Rcpp::exception("nlopt_set_xtol_abs"); - } -} +// +// void set_uniform_xtol_abs(NloptStruct * opt, double value) { +// if(nlopt_set_xtol_abs1(opt, value) != NLOPT_SUCCESS) { +// throw Rcpp::exception("nlopt_set_xtol_abs1"); +// } +// } +// void set_per_value_xtol_abs(NloptStruct * opt, const std::vector & xtol_abs) { +// if(xtol_abs.size() != nlopt_get_dimension(opt)) { +// throw Rcpp::exception("set_per_value_xtol_abs: parameter size mismatch"); +// } +// if(nlopt_set_xtol_abs(opt, xtol_abs.data()) != NLOPT_SUCCESS) { +// throw Rcpp::exception("nlopt_set_xtol_abs"); +// } +// } // --------------------------------------------------------------------------------------- // sanity test and example @@ -137,7 +144,7 @@ bool cpp_test_nlopt() { auto optimizer = new_nlopt_optimizer(config, x.size()); - set_uniform_xtol_abs(optimizer.get(), 0); + // set_uniform_xtol_abs(optimizer.get(), 0); // set_uniform_x_weights(optimizer.get(), 1.); check(nlopt_get_algorithm(optimizer.get()) == NLOPT_LD_LBFGS, "optim algorithm"); diff --git a/src/nlopt_wrapper.h b/src/nlopt_wrapper.h index c8be5d2a..857e6529 100644 --- a/src/nlopt_wrapper.h +++ b/src/nlopt_wrapper.h @@ -25,8 +25,8 @@ std::unique_ptr new_nlopt_optimizer(const Rcpp::List // Helpers to set xtol_abs (uniform or per-parameter packed array). // This is not done by new_nlopt_optimizer as it may require packing values, which must be user specified. -void set_uniform_xtol_abs(NloptStruct * opt, double value); -void set_per_value_xtol_abs(NloptStruct * opt, const std::vector & xtol_abs); +// void set_uniform_xtol_abs(NloptStruct * opt, double value); +// void set_per_value_xtol_abs(NloptStruct * opt, const std::vector & xtol_abs); // Helpers to set x_weights (uniform or per-parameter packed array). // This is not done by new_nlopt_optimizer as it may require packing values, which must be user specified. diff --git a/src/optim_diag_cov.cpp b/src/optim_diag_cov.cpp index 5e345665..01d9558c 100644 --- a/src/optim_diag_cov.cpp +++ b/src/optim_diag_cov.cpp @@ -34,24 +34,11 @@ Rcpp::List nlopt_optimize_diagonal( metadata.map(parameters.data()) = init_S; auto optimizer = new_nlopt_optimizer(config, parameters.size()); - if(config.containsElementNamed("xtol_abs")) { - SEXP value = config["xtol_abs"]; - if(Rcpp::is(value)) { - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(value)); - } else { - auto per_param_list = Rcpp::as(value); - auto packed = std::vector(metadata.packed_size); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["B"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["M"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["S"]); - set_per_value_xtol_abs(optimizer.get(), packed); - } - } - + std::vector objective_vec ; const double w_bar = accu(w); // Optimize - auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &w_bar](const double * params, double * grad) -> double { + auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &w_bar, &objective_vec](const double * params, double * grad) -> double { const arma::mat B = metadata.map(params); const arma::mat M = metadata.map(params); const arma::mat S = metadata.map(params); @@ -64,7 +51,10 @@ Rcpp::List nlopt_optimize_diagonal( metadata.map(grad) = (X.each_col() % w).t() * (A - Y); metadata.map(grad) = diagmat(w) * ((M.each_row() / diag_sigma) + A - Y); - metadata.map(grad) = diagmat(w) * (S.each_row() % pow(diag_sigma, -1) + S % A - pow(S, -1)); + metadata.map(grad) = diagmat(w) * (S.each_row() % pow(diag_sigma, -1) + S % A - pow(S, -1)) ; + + objective_vec.push_back(objective) ; + return objective; }; OptimizerResult result = minimize_objective_on_parameters(optimizer.get(), objective_and_grad, parameters); @@ -102,6 +92,7 @@ Rcpp::List nlopt_optimize_diagonal( Rcpp::Named("monitoring", Rcpp::List::create( Rcpp::Named("status", static_cast(result.status)), Rcpp::Named("backend", "nlopt"), + Rcpp::Named("objective", objective_vec), Rcpp::Named("iterations", result.nb_iterations) )) ); @@ -135,21 +126,10 @@ Rcpp::List nlopt_optimize_vestep_diagonal( metadata.map(parameters.data()) = init_S; auto optimizer = new_nlopt_optimizer(config, parameters.size()); - if(config.containsElementNamed("xtol_abs")) { - SEXP value = config["xtol_abs"]; - if(Rcpp::is(value)) { - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(value)); - } else { - auto per_param_list = Rcpp::as(value); - auto packed = std::vector(metadata.packed_size); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["M"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["S"]); - set_per_value_xtol_abs(optimizer.get(), packed); - } - } + std::vector objective_vec ; // Optimize - auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &B, &Omega](const double * params, double * grad) -> double { + auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &B, &Omega, &objective_vec](const double * params, double * grad) -> double { const arma::mat M = metadata.map(params); const arma::mat S = metadata.map(params); @@ -162,6 +142,9 @@ Rcpp::List nlopt_optimize_vestep_diagonal( metadata.map(grad) = diagmat(w) * (M * arma::diagmat(omega2) + A - Y); metadata.map(grad) = diagmat(w) * (S.each_row() % omega2.t() + S % A - pow(S, -1)); + + objective_vec.push_back(objective) ; + return objective; }; OptimizerResult result = minimize_objective_on_parameters(optimizer.get(), objective_and_grad, parameters); @@ -186,6 +169,7 @@ Rcpp::List nlopt_optimize_vestep_diagonal( Rcpp::Named("monitoring", Rcpp::List::create( Rcpp::Named("status", static_cast(result.status)), Rcpp::Named("backend", "nlopt"), + Rcpp::Named("objective", objective_vec), Rcpp::Named("iterations", result.nb_iterations) )) ); diff --git a/src/optim_fixed_cov.cpp b/src/optim_fixed_cov.cpp index 0ae4f97f..5ef92808 100644 --- a/src/optim_fixed_cov.cpp +++ b/src/optim_fixed_cov.cpp @@ -34,23 +34,11 @@ Rcpp::List nlopt_optimize_fixed( metadata.map(parameters.data()) = init_M; metadata.map(parameters.data()) = init_S; + // Optimize auto optimizer = new_nlopt_optimizer(config, parameters.size()); - if(config.containsElementNamed("xtol_abs")) { - SEXP value = config["xtol_abs"]; - if(Rcpp::is(value)) { - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(value)); - } else { - auto per_param_list = Rcpp::as(value); - auto packed = std::vector(metadata.packed_size); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["B"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["M"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["S"]); - set_per_value_xtol_abs(optimizer.get(), packed); - } - } + std::vector objective_vec ; - // Optimize - auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &Omega](const double * params, double * grad) -> double { + auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &Omega, &objective_vec](const double * params, double * grad) -> double { const arma::mat B = metadata.map(params); const arma::mat M = metadata.map(params); const arma::mat S = metadata.map(params); @@ -63,7 +51,10 @@ Rcpp::List nlopt_optimize_fixed( metadata.map(grad) = (X.each_col() % w).t() * (A - Y); metadata.map(grad) = diagmat(w) * (M * Omega + A - Y); - metadata.map(grad) = diagmat(w) * (S.each_row() % diagvec(Omega).t() + S % A - pow(S, -1)); + metadata.map(grad) = diagmat(w) * (S.each_row() % diagvec(Omega).t() + S % A - pow(S, -1)) ; + + objective_vec.push_back(objective) ; + return objective; }; OptimizerResult result = minimize_objective_on_parameters(optimizer.get(), objective_and_grad, parameters); @@ -94,6 +85,7 @@ Rcpp::List nlopt_optimize_fixed( Rcpp::Named("monitoring", Rcpp::List::create( Rcpp::Named("status", static_cast(result.status)), Rcpp::Named("backend", "nlopt"), + Rcpp::Named("objective", objective_vec), Rcpp::Named("iterations", result.nb_iterations) )) ); diff --git a/src/optim_full_cov.cpp b/src/optim_full_cov.cpp index f69d0092..0e43698a 100644 --- a/src/optim_full_cov.cpp +++ b/src/optim_full_cov.cpp @@ -33,53 +33,28 @@ Rcpp::List nlopt_optimize( metadata.map(parameters.data()) = init_M; metadata.map(parameters.data()) = init_S; + // Optimize auto optimizer = new_nlopt_optimizer(config, parameters.size()); - if(config.containsElementNamed("xtol_abs")) { - SEXP value = config["xtol_abs"]; - if(Rcpp::is(value)) { - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(value)); - } else { - auto per_param_list = Rcpp::as(value); - auto packed = std::vector(metadata.packed_size); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["B"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["M"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["S"]); - set_per_value_xtol_abs(optimizer.get(), packed); - } - } - - // if(config.containsElementNamed("x_weights")) { - // SEXP value = config["x_weights"]; - // if(Rcpp::is(value)) { - // set_uniform_x_weights(optimizer.get(), Rcpp::as(value)); - // } else { - // auto per_param_list = Rcpp::as(value); - // auto packed = std::vector(metadata.packed_size); - // set_from_r_sexp(metadata.map(packed.data()), per_param_list["B"]); - // set_from_r_sexp(metadata.map(packed.data()), per_param_list["M"]); - // set_from_r_sexp(metadata.map(packed.data()), per_param_list["S"]); - // set_per_value_x_weights(optimizer.get(), packed); - // } - // } - + std::vector objective_vec ; const double w_bar = accu(w); - // Optimize - auto objective_and_grad = [&metadata, &Y, &X, &O, &w, &w_bar](const double * params, double * grad) -> double { + auto objective_and_grad = [&metadata, &Y, &X, &O, &w, &w_bar, &objective_vec](const double * params, double * grad) -> double { const arma::mat B = metadata.map(params); const arma::mat M = metadata.map(params); const arma::mat S = metadata.map(params); const double w_bar = accu(w); - arma::mat S2 = S % S; - arma::mat Z = O + X * B + M; + arma::mat S2 = S % S ; + arma::mat Z = O + X * B + M ; arma::mat A = exp(Z + 0.5 * S2); arma::mat Omega = w_bar * inv_sympd(M.t() * (M.each_col() % w) + diagmat(w.t() * S2)); - double objective = accu(w.t() * (A - Y % Z - 0.5 * log(S2))) - 0.5 * w_bar * real(log_det(Omega)); + double objective = accu(w.t() * (A - Y % Z - 0.5 * trunc_log(S2))) - 0.5 * w_bar * real(log_det(Omega)); - metadata.map(grad) = (X.each_col() % w).t() * (A - Y); - metadata.map(grad) = diagmat(w) * (M * Omega + A - Y); - metadata.map(grad) = diagmat(w) * (S.each_row() % diagvec(Omega).t() + S % A - pow(S, -1)); + metadata.map(grad) = (X.each_col() % w).t() * (A - Y) ; + metadata.map(grad) = diagmat(w) * (M * Omega + A - Y) ; + metadata.map(grad) = diagmat(w) * (S.each_row() % diagvec(Omega).t() + S % A - pow(S, -1)) ; + + objective_vec.push_back(objective) ; return objective; }; @@ -114,6 +89,7 @@ Rcpp::List nlopt_optimize( Rcpp::Named("monitoring", Rcpp::List::create( Rcpp::Named("status", static_cast(result.status)), Rcpp::Named("backend", "nlopt"), + Rcpp::Named("objective", objective_vec), Rcpp::Named("iterations", result.nb_iterations) )) ); @@ -145,22 +121,11 @@ Rcpp::List nlopt_optimize_vestep( metadata.map(parameters.data()) = init_M; metadata.map(parameters.data()) = init_S; + // Optimize auto optimizer = new_nlopt_optimizer(config, parameters.size()); - if(config.containsElementNamed("xtol_abs")) { - SEXP value = config["xtol_abs"]; - if(Rcpp::is(value)) { - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(value)); - } else { - auto per_param_list = Rcpp::as(value); - auto packed = std::vector(metadata.packed_size); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["M"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["S"]); - set_per_value_xtol_abs(optimizer.get(), packed); - } - } + std::vector objective_vec ; - // Optimize - auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &B, &Omega](const double * params, double * grad) -> double { + auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &B, &Omega, &objective_vec](const double * params, double * grad) -> double { const arma::mat M = metadata.map(params); const arma::mat S = metadata.map(params); @@ -173,6 +138,8 @@ Rcpp::List nlopt_optimize_vestep( metadata.map(grad) = diagmat(w) * (M * Omega + A - Y); metadata.map(grad) = diagmat(w) * (S.each_row() % diagvec(Omega).t() + S % A - pow(S, -1)); + objective_vec.push_back(objective) ; + return objective; }; OptimizerResult result = minimize_objective_on_parameters(optimizer.get(), objective_and_grad, parameters); @@ -196,6 +163,7 @@ Rcpp::List nlopt_optimize_vestep( Rcpp::Named("monitoring", Rcpp::List::create( Rcpp::Named("status", static_cast(result.status)), Rcpp::Named("backend", "nlopt"), + Rcpp::Named("objective", objective_vec), Rcpp::Named("iterations", result.nb_iterations) )) ); diff --git a/src/optim_genet_cov.cpp b/src/optim_genet_cov.cpp index af7d56b6..356cfbb1 100644 --- a/src/optim_genet_cov.cpp +++ b/src/optim_genet_cov.cpp @@ -36,20 +36,6 @@ Rcpp::List nlopt_optimize_genetic_modeling( metadata.map(parameters.data()) = init_rho; auto optimizer = new_nlopt_optimizer(configuration, parameters.size()); - if(configuration.containsElementNamed("xtol_abs")) { - SEXP value = configuration["xtol_abs"]; - if(Rcpp::is(value)) { - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(value)); - } else { - auto per_param_list = Rcpp::as(value); - auto packed = std::vector(metadata.packed_size); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["Theta"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["M"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["S"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["rho"]); - set_per_value_xtol_abs(optimizer.get(), packed); - } - } // Some fixed quantities along optimization const double w_bar = accu(w); diff --git a/src/optim_rank_cov.cpp b/src/optim_rank_cov.cpp index 986d4e69..0f09f052 100644 --- a/src/optim_rank_cov.cpp +++ b/src/optim_rank_cov.cpp @@ -37,24 +37,11 @@ Rcpp::List nlopt_optimize_rank( metadata.map(parameters.data()) = init_M; metadata.map(parameters.data()) = init_S; + // Optimize auto optimizer = new_nlopt_optimizer(config, parameters.size()); - if(config.containsElementNamed("xtol_abs")) { - SEXP value = config["xtol_abs"]; - if(Rcpp::is(value)) { - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(value)); - } else { - auto per_param_list = Rcpp::as(value); - auto packed = std::vector(metadata.packed_size); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["B"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["C"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["M"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["S"]); - set_per_value_xtol_abs(optimizer.get(), packed); - } - } + std::vector objective_vec ; - // Optimize - auto objective_and_grad = [&metadata, &O, &X, &Y, &w](const double * params, double * grad) -> double { + auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &objective_vec](const double * params, double * grad) -> double { const arma::mat B = metadata.map(params); const arma::mat C = metadata.map(params); const arma::mat M = metadata.map(params); @@ -68,7 +55,10 @@ Rcpp::List nlopt_optimize_rank( metadata.map(grad) = (X.each_col() % w).t() * (A - Y); metadata.map(grad) = (diagmat(w) * (A - Y)).t() * M + (A.t() * (S2.each_col() % w)) % C; metadata.map(grad) = diagmat(w) * ((A - Y) * C + M); - metadata.map(grad) = diagmat(w) * (S - 1. / S + A * (C % C) % S); + metadata.map(grad) = diagmat(w) * (S - 1. / S + A * (C % C) % S) ; + + objective_vec.push_back(objective) ; + return objective; }; OptimizerResult result = minimize_objective_on_parameters(optimizer.get(), objective_and_grad, parameters); @@ -101,6 +91,7 @@ Rcpp::List nlopt_optimize_rank( Rcpp::Named("monitoring", Rcpp::List::create( Rcpp::Named("status", static_cast(result.status)), Rcpp::Named("backend", "nlopt"), + Rcpp::Named("objective", objective_vec), Rcpp::Named("iterations", result.nb_iterations) )) ); @@ -133,22 +124,11 @@ Rcpp::List nlopt_optimize_vestep_rank( metadata.map(parameters.data()) = init_M; metadata.map(parameters.data()) = init_S; + // Optimize auto optimizer = new_nlopt_optimizer(config, parameters.size()); - if(config.containsElementNamed("xtol_abs")) { - SEXP value = config["xtol_abs"]; - if(Rcpp::is(value)) { - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(value)); - } else { - auto per_param_list = Rcpp::as(value); - auto packed = std::vector(metadata.packed_size); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["M"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["S"]); - set_per_value_xtol_abs(optimizer.get(), packed); - } - } + std::vector objective_vec ; - // Optimize - auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &B, &C](const double * params, double * grad) -> double { + auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &B, &C, &objective_vec](const double * params, double * grad) -> double { const arma::mat M = metadata.map(params); const arma::mat S = metadata.map(params); @@ -160,6 +140,9 @@ Rcpp::List nlopt_optimize_vestep_rank( metadata.map(grad) = diagmat(w) * ((A - Y) * C + M); metadata.map(grad) = diagmat(w) * (S - 1. / S + A * (C % C) % S); + + objective_vec.push_back(objective) ; + return objective; }; OptimizerResult result = minimize_objective_on_parameters(optimizer.get(), objective_and_grad, parameters); @@ -182,6 +165,7 @@ Rcpp::List nlopt_optimize_vestep_rank( Rcpp::Named("monitoring", Rcpp::List::create( Rcpp::Named("status", static_cast(result.status)), Rcpp::Named("backend", "nlopt"), + Rcpp::Named("objective", objective_vec), Rcpp::Named("iterations", result.nb_iterations) )) ); diff --git a/src/optim_spherical.cpp b/src/optim_spherical.cpp index 3fd0d521..784099a9 100644 --- a/src/optim_spherical.cpp +++ b/src/optim_spherical.cpp @@ -33,25 +33,12 @@ Rcpp::List nlopt_optimize_spherical( metadata.map(parameters.data()) = init_M; metadata.map(parameters.data()) = init_S; + // Optimize auto optimizer = new_nlopt_optimizer(config, parameters.size()); - if(config.containsElementNamed("xtol_abs")) { - SEXP value = config["xtol_abs"]; - if(Rcpp::is(value)) { - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(value)); - } else { - auto per_param_list = Rcpp::as(value); - auto packed = std::vector(metadata.packed_size); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["B"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["M"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["S"]); - set_per_value_xtol_abs(optimizer.get(), packed); - } - } - const double w_bar = accu(w); + std::vector objective_vec ; - // Optimize - auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &w_bar](const double * params, double * grad) -> double { + auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &w_bar, &objective_vec](const double * params, double * grad) -> double { const arma::mat B = metadata.map(params); const arma::mat M = metadata.map(params); const arma::mat S = metadata.map(params); @@ -65,7 +52,9 @@ Rcpp::List nlopt_optimize_spherical( metadata.map(grad) = (X.each_col() % w).t() * (A - Y); metadata.map(grad) = diagmat(w) * (M / sigma2 + A - Y); - metadata.map(grad) = diagmat(w) * (S / sigma2 + S % A - pow(S, -1)); + metadata.map(grad) = diagmat(w) * (S / sigma2 + S % A - pow(S, -1)) ; + + objective_vec.push_back(objective) ; return objective; }; @@ -101,6 +90,7 @@ Rcpp::List nlopt_optimize_spherical( Rcpp::Named("monitoring", Rcpp::List::create( Rcpp::Named("status", static_cast(result.status)), Rcpp::Named("backend", "nlopt"), + Rcpp::Named("objective", objective_vec), Rcpp::Named("iterations", result.nb_iterations) )) ); @@ -132,22 +122,12 @@ Rcpp::List nlopt_optimize_vestep_spherical( metadata.map(parameters.data()) = init_M; metadata.map(parameters.data()) = init_S; - auto optimizer = new_nlopt_optimizer(config, parameters.size()); - if(config.containsElementNamed("xtol_abs")) { - SEXP value = config["xtol_abs"]; - if(Rcpp::is(value)) { - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(value)); - } else { - auto per_param_list = Rcpp::as(value); - auto packed = std::vector(metadata.packed_size); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["M"]); - set_from_r_sexp(metadata.map(packed.data()), per_param_list["S"]); - set_per_value_xtol_abs(optimizer.get(), packed); - } - } // Optimize - auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &B, &Omega](const double * params, double * grad) -> double { + auto optimizer = new_nlopt_optimizer(config, parameters.size()); + std::vector objective_vec ; + + auto objective_and_grad = [&metadata, &O, &X, &Y, &w, &B, &Omega, &objective_vec](const double * params, double * grad) -> double { const arma::mat M = metadata.map(params); const arma::mat S = metadata.map(params); @@ -161,6 +141,8 @@ Rcpp::List nlopt_optimize_vestep_spherical( metadata.map(grad) = diagmat(w) * (M / omega2 + A - Y); metadata.map(grad) = diagmat(w) * (S / omega2 + S % A - pow(S, -1)); + objective_vec.push_back(objective) ; + return objective; }; OptimizerResult result = minimize_objective_on_parameters(optimizer.get(), objective_and_grad, parameters); @@ -184,6 +166,7 @@ Rcpp::List nlopt_optimize_vestep_spherical( Rcpp::Named("monitoring", Rcpp::List::create( Rcpp::Named("status", static_cast(result.status)), Rcpp::Named("backend", "nlopt"), + Rcpp::Named("objective", objective_vec), Rcpp::Named("iterations", result.nb_iterations) )) ); diff --git a/src/optim_zi-pln.cpp b/src/optim_zi-pln.cpp index 559979cb..649fed9e 100644 --- a/src/optim_zi-pln.cpp +++ b/src/optim_zi-pln.cpp @@ -95,7 +95,6 @@ Rcpp::List optim_zipln_zipar_covar( metadata.map(parameters.data()) = init_B0; auto optimizer = new_nlopt_optimizer(configuration, parameters.size()); - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(configuration["xtol_abs"])); const arma::mat Xt_R = X0.t() * R; @@ -198,8 +197,6 @@ Rcpp::List optim_zipln_M( metadata.map(parameters.data()) = init_M; auto optimizer = new_nlopt_optimizer(configuration, parameters.size()); - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(configuration["xtol_abs"])); - const arma::mat X_B = X * B; // (n,p) const arma::mat O_S2 = O + 0.5 * S % S; // (n,p) @@ -241,8 +238,6 @@ Rcpp::List optim_zipln_S( metadata.map(parameters.data()) = init_S; auto optimizer = new_nlopt_optimizer(configuration, parameters.size()); - set_uniform_xtol_abs(optimizer.get(), Rcpp::as(configuration["xtol_abs"])); - const arma::mat O_M = O + M; // Optimize diff --git a/tests/testthat/test-pln.R b/tests/testthat/test-pln.R index a32e6753..54f38e5a 100644 --- a/tests/testthat/test-pln.R +++ b/tests/testthat/test-pln.R @@ -142,7 +142,7 @@ test_that("PLN: Check that all univariate PLN models are equivalent with the mul univariate_full <- lapply(1:p, function(j) { Abundance <- trichoptera$Abundance[, j, drop = FALSE] - PLN(Abundance ~ 1 + offset(log(Offset)), control = PLN_param(trace = 0)) + PLN(Abundance ~ 1 + offset(log(Offset)), control = PLN_param(trace = 0, config_optim = list(algorithm = "LBFGS"))) }) univariate_diagonal <- lapply(1:p, function(j) { diff --git a/tests/testthat/test-ziplnfit.R b/tests/testthat/test-ziplnfit.R index dbc5bfd5..135a3032 100644 --- a/tests/testthat/test-ziplnfit.R +++ b/tests/testthat/test-ziplnfit.R @@ -127,7 +127,7 @@ test_that("ZIPLN fit: check sparse output and plot", { expect_is(myPLNfit, "ZIPLNfit") expect_equal(myPLNfit$vcov_model, "sparse") - expect_true(igraph::is.igraph(myPLNfit$plot_network(output = "igraph", plot = FALSE))) + expect_true(igraph::is_igraph(myPLNfit$plot_network(output = "igraph", plot = FALSE))) expect_true(inherits(myPLNfit$plot_network(output = "corrplot", plot = FALSE), "Matrix")) })