Skip to content

Commit

Permalink
Merge pull request #122 from PLN-team/diag-optim-PLN
Browse files Browse the repository at this point in the history 
simplication in optim and output objective in optim_par
  • Loading branch information
@jchiquet
jchiquet authored Mar 27, 2024
2 parents d6b0d1a + 461e9d6 commit 59868f2
Show file tree
Hide file tree
Showing 11 changed files with 94 additions and 195 deletions.
37 changes: 22 additions & 15 deletions src/nlopt_wrapper.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,11 +50,18 @@ std::unique_ptr<NloptStruct, NloptDeleter> 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<double>(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<double>(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<double>(config["ftol_abs"])) != NLOPT_SUCCESS) {
throw Rcpp::exception("nlopt_set_ftol_abs");
Expand Down Expand Up @@ -93,20 +100,20 @@ std::unique_ptr<NloptStruct, NloptDeleter> 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<double> & 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<double> & 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
Expand Down Expand Up @@ -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");
Expand Down
4 changes: 2 additions & 2 deletions src/nlopt_wrapper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,8 @@ std::unique_ptr<NloptStruct, NloptDeleter> 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<double> & xtol_abs);
// void set_uniform_xtol_abs(NloptStruct * opt, double value);
// void set_per_value_xtol_abs(NloptStruct * opt, const std::vector<double> & 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.
Expand Down
42 changes: 13 additions & 29 deletions src/optim_diag_cov.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,24 +34,11 @@ Rcpp::List nlopt_optimize_diagonal(
metadata.map<S_ID>(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<double>(value)) {
set_uniform_xtol_abs(optimizer.get(), Rcpp::as<double>(value));
} else {
auto per_param_list = Rcpp::as<Rcpp::List>(value);
auto packed = std::vector<double>(metadata.packed_size);
set_from_r_sexp(metadata.map<B_ID>(packed.data()), per_param_list["B"]);
set_from_r_sexp(metadata.map<M_ID>(packed.data()), per_param_list["M"]);
set_from_r_sexp(metadata.map<S_ID>(packed.data()), per_param_list["S"]);
set_per_value_xtol_abs(optimizer.get(), packed);
}
}

std::vector<double> 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<B_ID>(params);
const arma::mat M = metadata.map<M_ID>(params);
const arma::mat S = metadata.map<S_ID>(params);
Expand All @@ -64,7 +51,10 @@ Rcpp::List nlopt_optimize_diagonal(

metadata.map<B_ID>(grad) = (X.each_col() % w).t() * (A - Y);
metadata.map<M_ID>(grad) = diagmat(w) * ((M.each_row() / diag_sigma) + A - Y);
metadata.map<S_ID>(grad) = diagmat(w) * (S.each_row() % pow(diag_sigma, -1) + S % A - pow(S, -1));
metadata.map<S_ID>(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);
Expand Down Expand Up @@ -102,6 +92,7 @@ Rcpp::List nlopt_optimize_diagonal(
Rcpp::Named("monitoring", Rcpp::List::create(
Rcpp::Named("status", static_cast<int>(result.status)),
Rcpp::Named("backend", "nlopt"),
Rcpp::Named("objective", objective_vec),
Rcpp::Named("iterations", result.nb_iterations)
))
);
Expand Down Expand Up @@ -135,21 +126,10 @@ Rcpp::List nlopt_optimize_vestep_diagonal(
metadata.map<S_ID>(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<double>(value)) {
set_uniform_xtol_abs(optimizer.get(), Rcpp::as<double>(value));
} else {
auto per_param_list = Rcpp::as<Rcpp::List>(value);
auto packed = std::vector<double>(metadata.packed_size);
set_from_r_sexp(metadata.map<M_ID>(packed.data()), per_param_list["M"]);
set_from_r_sexp(metadata.map<S_ID>(packed.data()), per_param_list["S"]);
set_per_value_xtol_abs(optimizer.get(), packed);
}
}
std::vector<double> 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<M_ID>(params);
const arma::mat S = metadata.map<S_ID>(params);

Expand All @@ -162,6 +142,9 @@ Rcpp::List nlopt_optimize_vestep_diagonal(

metadata.map<M_ID>(grad) = diagmat(w) * (M * arma::diagmat(omega2) + A - Y);
metadata.map<S_ID>(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);
Expand All @@ -186,6 +169,7 @@ Rcpp::List nlopt_optimize_vestep_diagonal(
Rcpp::Named("monitoring", Rcpp::List::create(
Rcpp::Named("status", static_cast<int>(result.status)),
Rcpp::Named("backend", "nlopt"),
Rcpp::Named("objective", objective_vec),
Rcpp::Named("iterations", result.nb_iterations)
))
);
Expand Down
24 changes: 8 additions & 16 deletions src/optim_fixed_cov.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,23 +34,11 @@ Rcpp::List nlopt_optimize_fixed(
metadata.map<M_ID>(parameters.data()) = init_M;
metadata.map<S_ID>(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<double>(value)) {
set_uniform_xtol_abs(optimizer.get(), Rcpp::as<double>(value));
} else {
auto per_param_list = Rcpp::as<Rcpp::List>(value);
auto packed = std::vector<double>(metadata.packed_size);
set_from_r_sexp(metadata.map<B_ID>(packed.data()), per_param_list["B"]);
set_from_r_sexp(metadata.map<M_ID>(packed.data()), per_param_list["M"]);
set_from_r_sexp(metadata.map<S_ID>(packed.data()), per_param_list["S"]);
set_per_value_xtol_abs(optimizer.get(), packed);
}
}
std::vector<double> 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<B_ID>(params);
const arma::mat M = metadata.map<M_ID>(params);
const arma::mat S = metadata.map<S_ID>(params);
Expand All @@ -63,7 +51,10 @@ Rcpp::List nlopt_optimize_fixed(

metadata.map<B_ID>(grad) = (X.each_col() % w).t() * (A - Y);
metadata.map<M_ID>(grad) = diagmat(w) * (M * Omega + A - Y);
metadata.map<S_ID>(grad) = diagmat(w) * (S.each_row() % diagvec(Omega).t() + S % A - pow(S, -1));
metadata.map<S_ID>(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);
Expand Down Expand Up @@ -94,6 +85,7 @@ Rcpp::List nlopt_optimize_fixed(
Rcpp::Named("monitoring", Rcpp::List::create(
Rcpp::Named("status", static_cast<int>(result.status)),
Rcpp::Named("backend", "nlopt"),
Rcpp::Named("objective", objective_vec),
Rcpp::Named("iterations", result.nb_iterations)
))
);
Expand Down
68 changes: 18 additions & 50 deletions src/optim_full_cov.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,53 +33,28 @@ Rcpp::List nlopt_optimize(
metadata.map<M_ID>(parameters.data()) = init_M;
metadata.map<S_ID>(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<double>(value)) {
set_uniform_xtol_abs(optimizer.get(), Rcpp::as<double>(value));
} else {
auto per_param_list = Rcpp::as<Rcpp::List>(value);
auto packed = std::vector<double>(metadata.packed_size);
set_from_r_sexp(metadata.map<B_ID>(packed.data()), per_param_list["B"]);
set_from_r_sexp(metadata.map<M_ID>(packed.data()), per_param_list["M"]);
set_from_r_sexp(metadata.map<S_ID>(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<double>(value)) {
// set_uniform_x_weights(optimizer.get(), Rcpp::as<double>(value));
// } else {
// auto per_param_list = Rcpp::as<Rcpp::List>(value);
// auto packed = std::vector<double>(metadata.packed_size);
// set_from_r_sexp(metadata.map<B_ID>(packed.data()), per_param_list["B"]);
// set_from_r_sexp(metadata.map<M_ID>(packed.data()), per_param_list["M"]);
// set_from_r_sexp(metadata.map<S_ID>(packed.data()), per_param_list["S"]);
// set_per_value_x_weights(optimizer.get(), packed);
// }
// }

std::vector<double> 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<B_ID>(params);
const arma::mat M = metadata.map<M_ID>(params);
const arma::mat S = metadata.map<S_ID>(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<B_ID>(grad) = (X.each_col() % w).t() * (A - Y);
metadata.map<M_ID>(grad) = diagmat(w) * (M * Omega + A - Y);
metadata.map<S_ID>(grad) = diagmat(w) * (S.each_row() % diagvec(Omega).t() + S % A - pow(S, -1));
metadata.map<B_ID>(grad) = (X.each_col() % w).t() * (A - Y) ;
metadata.map<M_ID>(grad) = diagmat(w) * (M * Omega + A - Y) ;
metadata.map<S_ID>(grad) = diagmat(w) * (S.each_row() % diagvec(Omega).t() + S % A - pow(S, -1)) ;

objective_vec.push_back(objective) ;

return objective;
};
Expand Down Expand Up @@ -114,6 +89,7 @@ Rcpp::List nlopt_optimize(
Rcpp::Named("monitoring", Rcpp::List::create(
Rcpp::Named("status", static_cast<int>(result.status)),
Rcpp::Named("backend", "nlopt"),
Rcpp::Named("objective", objective_vec),
Rcpp::Named("iterations", result.nb_iterations)
))
);
Expand Down Expand Up @@ -145,22 +121,11 @@ Rcpp::List nlopt_optimize_vestep(
metadata.map<M_ID>(parameters.data()) = init_M;
metadata.map<S_ID>(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<double>(value)) {
set_uniform_xtol_abs(optimizer.get(), Rcpp::as<double>(value));
} else {
auto per_param_list = Rcpp::as<Rcpp::List>(value);
auto packed = std::vector<double>(metadata.packed_size);
set_from_r_sexp(metadata.map<M_ID>(packed.data()), per_param_list["M"]);
set_from_r_sexp(metadata.map<S_ID>(packed.data()), per_param_list["S"]);
set_per_value_xtol_abs(optimizer.get(), packed);
}
}
std::vector<double> 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<M_ID>(params);
const arma::mat S = metadata.map<S_ID>(params);

Expand All @@ -173,6 +138,8 @@ Rcpp::List nlopt_optimize_vestep(
metadata.map<M_ID>(grad) = diagmat(w) * (M * Omega + A - Y);
metadata.map<S_ID>(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);
Expand All @@ -196,6 +163,7 @@ Rcpp::List nlopt_optimize_vestep(
Rcpp::Named("monitoring", Rcpp::List::create(
Rcpp::Named("status", static_cast<int>(result.status)),
Rcpp::Named("backend", "nlopt"),
Rcpp::Named("objective", objective_vec),
Rcpp::Named("iterations", result.nb_iterations)
))
);
Expand Down
14 changes: 0 additions & 14 deletions src/optim_genet_cov.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,20 +36,6 @@ Rcpp::List nlopt_optimize_genetic_modeling(
metadata.map<RHO_ID>(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<double>(value)) {
set_uniform_xtol_abs(optimizer.get(), Rcpp::as<double>(value));
} else {
auto per_param_list = Rcpp::as<Rcpp::List>(value);
auto packed = std::vector<double>(metadata.packed_size);
set_from_r_sexp(metadata.map<THETA_ID>(packed.data()), per_param_list["Theta"]);
set_from_r_sexp(metadata.map<M_ID>(packed.data()), per_param_list["M"]);
set_from_r_sexp(metadata.map<S_ID>(packed.data()), per_param_list["S"]);
set_from_r_sexp(metadata.map<RHO_ID>(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);
Expand Down
Loading

0 comments on commit 59868f2

Please sign in to comment.