Louis BECQUEY / biorseo

Go to a project
Toggle navigation Toggle navigation pinning
Louis BECQUEY
16f5a95a 1 parent 990ba67f

Optimisations

Inline Side-by-side
Showing 3 changed files with 96 additions and 79 deletions
......@@ -18,17 +18,14 @@ using std::make_pair;
using std::vector;
uint MOIP::obj_to_solve_ = 1;
double MOIP::precision_ = 1e-7;
double MOIP::epsilon_ = 1e-4;
double MOIP::precision_ = 1e-5;
unsigned getNumConstraints(IloModel& m)
{
unsigned count = 0;
IloModel::Iterator iter(m);
while (iter.ok()) {
if ((*iter).asConstraint().getImpl()) {
++count;
}
if ((*iter).asConstraint().getImpl()) ++count;
++iter;
}
return count;
......@@ -133,7 +130,7 @@ bool MOIP::is_undominated_yet(const SecondaryStructure& s)
return true;
}
SecondaryStructure MOIP::solve_objective(int o, double min, double max, bool below)
SecondaryStructure MOIP::solve_objective(int o, double min, double max, const vector<IloConstraint>& F)
{
// Solves one of the objectives, under constraint that the other should be in [min, max]
......@@ -144,37 +141,18 @@ SecondaryStructure MOIP::solve_objective(int o, double min, double max, bool bel
max = max - min;
}
if (verbose_) {
cout << std::setprecision(8) << "\nSolving objective function " << o << ", ";
if (below)
cout << "below " << max;
else
cout << "on top of " << min;
cout << ": Obj" << 3 - o << " being in [" << min - precision_ << ", " << max + precision_ << "]..." << endl;
}
IloObjective obj;
IloRange bounds;
// gather known solutions in the search zone to forbid them
if (verbose_)
cout << "\t>forbidding solutions found in [" << std::setprecision(10) << (min) << ", " << (max) << ']' << endl;
vector<IloConstraint> F;
for (const SecondaryStructure& prev : pareto_)
if ((min) <= prev.get_objective_score(3 - o) and prev.get_objective_score(3 - o) <= (max)) {
F.push_back(prev.forbid_this_);
// if (verbose_) cout << "\t\t>forbidding " << prev.to_string() << endl;
}
if (verbose_) cout << "\t>forbidding " << F.size() << " solutions already found in that zone" << endl;
// impose the bounds and the objective
switch (o) {
case 1:
obj = IloMaximize(env_, obj1);
bounds = IloRange(env_, min - precision_, obj2, max + precision_);
bounds = IloRange(env_, min, obj2, max);
break;
case 2:
obj = IloMaximize(env_, obj2);
bounds = IloRange(env_, min - precision_, obj1, max + precision_);
bounds = IloRange(env_, min, obj1, max);
break;
}
model_.add(obj);
......@@ -398,9 +376,9 @@ void MOIP::define_problem_constraints(void)
}
}
void MOIP::search_between(double lambdaMin, double lambdaMax, bool below)
void MOIP::search_between(double lambdaMin, double lambdaMax, const vector<IloConstraint>& F_)
{
SecondaryStructure s = solve_objective(obj_to_solve_, lambdaMin, lambdaMax, below);
SecondaryStructure s = solve_objective(obj_to_solve_, lambdaMin, lambdaMax, F_);
if (!s.is_empty_structure) { // A solution has been found
// Attribute the correct pareto set label
......@@ -419,7 +397,7 @@ void MOIP::search_between(double lambdaMin, double lambdaMax, bool below)
s.set_pareto_set(max + 1);
// if (verbose_) cout << "\t>belongs to Pareto set " << s.get_pareto_set() << endl;
if (s.get_pareto_set() <= n_sets_ + 1) {
if (s.get_pareto_set() <= n_sets_) {
// adding the SecondaryStructure s to the set pareto_
add_solution(s);
......@@ -442,32 +420,40 @@ void MOIP::search_between(double lambdaMin, double lambdaMax, bool below)
}
}
// if (exists_horizontal_outdated_labels(s))
// for (vector<SecondaryStructure>::iterator x = pareto_.end() - 2; x >= pareto_.begin(); x--)
// if (
// abs(x->get_objective_score(3 - obj_to_solve_) - s.get_objective_score(3 - obj_to_solve_)) < precision_ and
// precision_ < s.get_objective_score(obj_to_solve_) - x->get_objective_score(obj_to_solve_)) {
// uint k = x->get_pareto_set();
// if (k <= n_sets_) {
// if (verbose_)
// cout << "\t>moving a structure from Pareto set " << k << " to " << k + 1 << endl;
// x->set_pareto_set(k + 1);
// } else {
// if (verbose_)
// cout << "\t>removing structure from Pareto set " << k << ":\tobj" << 3 - obj_to_solve_
// << '=' << x->get_objective_score(3 - obj_to_solve_) << endl;
// pareto_.erase(x);
// }
// }
// search below and on top of s
search_between(s.get_objective_score(3 - obj_to_solve_) + epsilon_, lambdaMax, false);
if (s.get_pareto_set() <= n_sets_)
search_between(lambdaMin, s.get_objective_score(3 - obj_to_solve_), true);
} else {
if (verbose_) cout << "\t>solution ignored." << endl;
double min = s.get_objective_score(3 - obj_to_solve_) + precision_;
double max = lambdaMax;
vector<IloConstraint> F = forbid_solutions_between(min, max);
// search on top
if (verbose_)
cout << std::setprecision(8) << "\nSolving objective function " << obj_to_solve_ << ", on top of "
<< s.get_objective_score(3 - obj_to_solve_) << ": Obj" << 3 - obj_to_solve_ << " being in ["
<< min << ", " << max << "]..." << endl
<< "\t>forbidding " << F.size() << " solutions found in [" << std::setprecision(10)
<< min - precision_ << ", " << max + precision_ << ']' << endl;
search_between(min, max, F);
if (s.get_pareto_set() <= n_sets_) {
min = lambdaMin;
max = s.get_objective_score(3 - obj_to_solve_);
F.clear();
F = forbid_solutions_between(min, max);
// search below
if (verbose_)
cout << std::setprecision(8) << "\nSolving objective function " << obj_to_solve_ << ", below (or eq. to) "
<< max << ": Obj" << 3 - obj_to_solve_ << " being in [" << min << ", " << max << "]..." << endl
<< "\t>forbidding " << F.size() << " solutions found in [" << std::setprecision(10)
<< min - precision_ << ", " << max + precision_ << ']' << endl;
search_between(min, max, F);
}
}
} else {
if (verbose_) cout << "\t>solution ignored." << endl;
}
}
......@@ -503,21 +489,21 @@ bool MOIP::exists_horizontal_outdated_labels(const SecondaryStructure& s) const
}
vector<IloConstraint> MOIP::forbid_solutions_between(double min, double max)
{
vector<IloConstraint> F;
// gather known solutions in the search zone to forbid them
for (const SecondaryStructure& prev : pareto_)
if (min - precision_ <= prev.get_objective_score(3 - obj_to_solve_) and prev.get_objective_score(3 - obj_to_solve_) <= max + precision_)
F.push_back(prev.forbid_this_);
return F;
}
void MOIP::add_solution(const SecondaryStructure& s)
{
// vector<size_t> to_remove;
// for (uint i = 0; i < pareto_.size(); i++)
// if (s > pareto_[i]) { // A solution from the set is now dominated
// // This should only happen in the case some structures have the same optimal Obj1 value.
// if (verbose_) cout << "\t>removing structure from Pareto set : " << pareto_[i].to_string() << endl;
// to_remove.push_back(i);
// }
// if (to_remove.size()) {
// for (size_t i = to_remove.size() - 1; i != 0; i--) pareto_.erase(pareto_.begin() + to_remove[i]);
// pareto_.erase(pareto_.begin() + to_remove[0]);
// }
if (verbose_)
cout << "\t>adding structure to Pareto set " << s.get_pareto_set() << " : " << s.to_string() << endl;
// if (verbose_)
cout << "\t>adding structure to Pareto set " << s.get_pareto_set() << " : " << s.to_string() << endl;
pareto_.push_back(s);
}
......
......@@ -16,14 +16,15 @@ class MOIP
MOIP(void);
MOIP(const RNA& rna, const vector<Motif>& motifSites, uint nsets, float pthreshold, bool verbose);
~MOIP(void);
SecondaryStructure solve_objective(int o, double min, double max, bool below);
SecondaryStructure solve_objective(int o, double min, double max, const vector<IloConstraint>& F);
SecondaryStructure solve_objective(int o);
uint get_n_solutions(void) const;
const SecondaryStructure& solution(uint i) const;
void search_between(double lambdaMin, double lambdaMax, bool below);
void search_between(double lambdaMin, double lambdaMax, const vector<IloConstraint>& F_);
bool allowed_basepair(size_t u, size_t v) const;
void add_solution(const SecondaryStructure& s);
void remove_solution(uint i);
vector<IloConstraint> forbid_solutions_between(double min, double max);
static uint obj_to_solve_; // What objective do you prefer to solve in mono-objective portions of the algorithm ?
static double precision_; // decimals to keep in objective values, to avoid numerical issues. otherwise, solution with objective 5.0000000009 dominates solution with 5.0 =(
static double epsilon_;
......@@ -66,5 +67,8 @@ inline uint MOIP::get_n_solutions(void) const { return pare
inline const SecondaryStructure& MOIP::solution(uint i) const { return pareto_[i]; }
inline IloNumExprArg& MOIP::y(size_t u, size_t v) { return basepair_dv_[get_yuv_index(u, v)]; }
inline IloNumExprArg& MOIP::C(size_t x, size_t i) { return insertion_dv_[get_Cpxi_index(x, i)]; }
inline SecondaryStructure MOIP::solve_objective(int o) { return solve_objective(o, 0, rna_.get_RNA_length(), false); }
inline SecondaryStructure MOIP::solve_objective(int o)
{
return solve_objective(o, 0, rna_.get_RNA_length(), vector<IloConstraint>());
}
#endif // MOIP_H_
\ No newline at end of file
......
......@@ -118,8 +118,8 @@ int main(int argc, char* argv[])
MOIP myMOIP = MOIP(myRNA, posInsertionSites, 1, theta_p_threshold, verbose);
try {
bestSSO1 = myMOIP.solve_objective(1, -__DBL_MAX__, __DBL_MAX__, false);
bestSSO2 = myMOIP.solve_objective(2, -__DBL_MAX__, __DBL_MAX__, false);
bestSSO1 = myMOIP.solve_objective(1, -__DBL_MAX__, __DBL_MAX__, vector<IloConstraint>());
bestSSO2 = myMOIP.solve_objective(2, -__DBL_MAX__, __DBL_MAX__, vector<IloConstraint>());
bestSSO1.set_pareto_set(1);
bestSSO2.set_pareto_set(1);
if (verbose) {
......@@ -127,20 +127,47 @@ int main(int argc, char* argv[])
cout << "Best solution according to objective 2 :" << bestSSO2.to_string() << endl;
}
// extend to the whole pareto set
double min, max;
vector<IloConstraint> F;
// extend the Pareto set on top
if (MOIP::obj_to_solve_ == 1) {
myMOIP.add_solution(bestSSO1);
min = bestSSO1.get_objective_score(2) + MOIP::precision_;
max = bestSSO2.get_objective_score(2);
if (verbose) cout << endl << "Solving obj1 on top of best solution 1." << endl;
myMOIP.search_between(bestSSO1.get_objective_score(2) + MOIP::epsilon_, bestSSO2.get_objective_score(2), false);
if (verbose) cout << endl << "Solving obj1 below best solution 1." << endl;
myMOIP.search_between(-__DBL_MAX__, bestSSO1.get_objective_score(2), true);
} else {
myMOIP.add_solution(bestSSO2);
min = bestSSO2.get_objective_score(1) + MOIP::precision_;
max = bestSSO1.get_objective_score(1);
if (verbose) cout << endl << "Solving obj2 on top of best solution 2." << endl;
myMOIP.search_between(bestSSO2.get_objective_score(1) + MOIP::epsilon_, bestSSO1.get_objective_score(1), false);
}
if (verbose)
cout << std::setprecision(8) << "\nSolving objective function " << MOIP::obj_to_solve_ << ", on top of "
<< min << ": Obj" << 3 - MOIP::obj_to_solve_ << " being in [" << min << ", " << max << "]..." << endl;
myMOIP.search_between(min, max, F); // F is empty
// extend the Pareto set below
if (MOIP::obj_to_solve_ == 1) {
if (verbose) cout << endl << "Solving obj1 below best solution 1." << endl;
min = -__DBL_MAX__;
max = bestSSO1.get_objective_score(2);
} else {
if (verbose) cout << endl << "Solving obj2 below best solution 2." << endl;
myMOIP.search_between(-__DBL_MAX__, bestSSO2.get_objective_score(1), true);
min = -__DBL_MAX__;
max = bestSSO2.get_objective_score(1);
}
F.clear();
F = myMOIP.forbid_solutions_between(min, max);
if (verbose)
cout << std::setprecision(8) << "\nSolving objective function " << MOIP::obj_to_solve_ << ", below (or eq. to) "
<< max << ": Obj" << 3 - MOIP::obj_to_solve_ << " being in [" << min << ", " << max << "]..." << endl
<< "\t>forbidding " << F.size() << " solutions found in [" << std::setprecision(10)
<< min - MOIP::precision_ << ", " << max + MOIP::precision_ << ']' << endl;
myMOIP.search_between(min, max, F);
} catch (IloAlgorithm::NotExtractedException& e) {
cerr << e << endl;
exit(EXIT_FAILURE);
......