Louis BECQUEY / biorseo

Go to a project
Toggle navigation Toggle navigation pinning
Louis BECQUEY
4bfcaf37 1 parent ff314aea

Re-forbid every found solution (faster and more precise)

Inline Side-by-side
Showing 4 changed files with 35 additions and 101 deletions
......@@ -31,15 +31,15 @@ unsigned getNumConstraints(IloModel& m)
return count;
}
MOIP::MOIP() {}
MOIP::MOIP(const RNA& rna, const vector<Motif>& insertionSites, float pthreshold, bool verbose)
: verbose_{verbose}, rna_(rna), insertion_sites_(insertionSites), theta_{pthreshold}
MOIP::MOIP(const RNA& rna, const vector<Motif>& insertionSites, float theta, bool verbose)
: verbose_{verbose}, rna_(rna), insertion_sites_(insertionSites)
{
if (verbose_) rna_.print_basepair_p_matrix(theta_);
if (verbose_) rna_.print_basepair_p_matrix(theta);
if (verbose_) cout << "defining problem decision variables..." << endl;
basepair_dv_ = IloNumVarArray(env_);
......@@ -51,7 +51,7 @@ MOIP::MOIP(const RNA& rna, const vector<Motif>& insertionSites, float pthreshold
index_of_yuv_ = vector<vector<size_t>>(rna_.get_RNA_length() - 6, vector<size_t>(0));
for (u = 0; u < rna_.get_RNA_length() - 6; u++)
for (v = u + 4; v < rna_.get_RNA_length(); v++) // A basepair is possible iff v > u+3
if (rna_.get_pij(u, v) > theta_) {
if (rna_.get_pij(u, v) > theta) {
if (verbose_) cout << u << '-' << v << " ";
index_of_yuv_[u].push_back(c);
c++;
......@@ -130,7 +130,7 @@ bool MOIP::is_undominated_yet(const SecondaryStructure& s)
return true;
}
SecondaryStructure MOIP::solve_objective(int o, double min, double max, const IloConstraintArray& F)
SecondaryStructure MOIP::solve_objective(int o, double min, double max)
{
// Solves one of the objectives, under constraint that the other should be in [min, max]
......@@ -141,10 +141,9 @@ SecondaryStructure MOIP::solve_objective(int o, double min, double max, const Il
max = max - min;
}
// impose the bounds and the objective
IloObjective obj;
IloRange bounds;
// impose the bounds and the objective
switch (o) {
case 1:
obj = IloMaximize(env_, obj1);
......@@ -157,8 +156,6 @@ SecondaryStructure MOIP::solve_objective(int o, double min, double max, const Il
}
model_.add(obj);
model_.add(bounds);
// cout << "\t>adding bounds: " << bounds << endl;
model_.add(F);
IloCplex cplex_ = IloCplex(model_);
cplex_.setOut(env_.getNullStream());
......@@ -169,7 +166,6 @@ SecondaryStructure MOIP::solve_objective(int o, double min, double max, const Il
// Removing the objective from the model_
model_.remove(obj);
model_.remove(bounds);
model_.remove(F);
return SecondaryStructure(true);
}
......@@ -192,12 +188,11 @@ SecondaryStructure MOIP::solve_objective(int o, double min, double max, const Il
if (cplex_.getValue(y(u, v)) > 0.5) best_ss.set_basepair(u, v);
best_ss.sort(); // order the basepairs in the vector
// truncate the result of cplex's computation which is full of numerical unprecisions
best_ss.set_objective_score(2, cplex_.getValue(obj2));
best_ss.set_objective_score(1, cplex_.getValue(obj1));
// if (verbose_) cout << "\t\t>building the IP forbidding condition..." << endl;
// Forbidding to find best_ss later : save a constraint
// Forbidding to find best_ss later
IloExpr c(env_);
for (uint d = 0; d < insertion_dv_.getSize(); d++)
if (cplex_.getValue(insertion_dv_[d]) > 0.5)
......@@ -209,15 +204,11 @@ SecondaryStructure MOIP::solve_objective(int o, double min, double max, const Il
c += IloNum(1) - basepair_dv_[d];
else
c += basepair_dv_[d];
best_ss.forbid_this_ = (c >= IloNum(1));
// Removing the objective from the model_
// if (verbose_) cout << "\t\t>removing temporary constraints..." << endl;
model_.remove(obj);
model_.remove(bounds);
model_.remove(F);
model_.add(c >= IloNum(1));
// exit
model_.remove(bounds);
model_.remove(obj);
return best_ss;
}
......@@ -379,9 +370,9 @@ void MOIP::define_problem_constraints(void)
}
}
void MOIP::search_between(double lambdaMin, double lambdaMax, const IloConstraintArray& F_)
void MOIP::search_between(double lambdaMin, double lambdaMax)
{
SecondaryStructure s = solve_objective(obj_to_solve_, lambdaMin, lambdaMax, F_);
SecondaryStructure s = solve_objective(obj_to_solve_, lambdaMin, lambdaMax);
if (!s.is_empty_structure) { // A solution has been found
// if the solution is dominated, ignore it
......@@ -407,19 +398,14 @@ void MOIP::search_between(double lambdaMin, double lambdaMax, const IloConstrain
}
// search on top
double min = s.get_objective_score(3 - obj_to_solve_) + precision_;
double truemin = s.get_objective_score(3 - obj_to_solve_);
double max = lambdaMax;
IloConstraintArray F = forbid_solutions_between(min, max);
double min = s.get_objective_score(3 - obj_to_solve_) + precision_;
double max = lambdaMax;
if (verbose_)
cout << std::setprecision(-log10(precision_) + 4) << "\nSolving objective function " << obj_to_solve_
<< ", on top of " << s.get_objective_score(3 - obj_to_solve_) << ": Obj" << 3 - obj_to_solve_
<< " being in [" << std::setprecision(-log10(precision_) + 4) << min << ", "
<< std::setprecision(-log10(precision_) + 4) << max << "]..." << endl
<< "\t>forbidding " << F.getSize() << " solutions found in [" << std::setprecision(-log10(precision_) + 4)
<< std::setprecision(-log10(precision_) + 4) << min - precision_ << ", "
<< std::setprecision(-log10(precision_) + 4) << max + precision_ << ']' << endl;
search_between(min, max, F);
<< std::setprecision(-log10(precision_) + 4) << max << "]..." << endl;
search_between(min, max);
if (std::abs(max - min) - precision_ > precision_) {
......@@ -427,17 +413,12 @@ void MOIP::search_between(double lambdaMin, double lambdaMax, const IloConstrain
// search below
min = lambdaMin;
max = s.get_objective_score(3 - obj_to_solve_);
F.clear();
F = forbid_solutions_between(min, max);
if (verbose_)
cout << std::setprecision(-log10(precision_) + 4) << "\nSolving objective function " << obj_to_solve_
<< ", below (or eq. to) " << max << ": Obj" << 3 - obj_to_solve_ << " being in ["
<< std::setprecision(-log10(precision_) + 4) << min << ", "
<< std::setprecision(-log10(precision_) + 4) << max << "]..." << endl
<< "\t>forbidding " << F.getSize() << " solutions found in ["
<< std::setprecision(-log10(precision_) + 4) << min - precision_ << ", "
<< std::setprecision(-log10(precision_) + 4) << max + precision_ << ']' << endl;
search_between(min, max, F);
<< std::setprecision(-log10(precision_) + 4) << max << "]..." << endl;
search_between(min, max);
}
} else {
......@@ -476,18 +457,6 @@ bool MOIP::exists_horizontal_outdated_labels(const SecondaryStructure& s) const
return result;
}
IloConstraintArray MOIP::forbid_solutions_between(double min, double max)
{
IloConstraintArray F(env_);
// 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.add(prev.forbid_this_);
return F;
}
void MOIP::add_solution(const SecondaryStructure& s)
{
if (verbose_) cout << "\t>adding structure to Pareto set :\t" << s.to_string() << endl;
......@@ -515,9 +484,8 @@ bool MOIP::allowed_basepair(size_t u, size_t v) const
if (a >= rna_.get_RNA_length() - 6) return false;
if (b >= rna_.get_RNA_length()) return false;
if (get_yuv_index(a, b) == rna_.get_RNA_length() * rna_.get_RNA_length() + 1)
return false; // not allowed because proba < theta_
return false; // not allowed because proba < theta
return true;
}
void MOIP::remove_solution(uint i) { pareto_.erase(pareto_.begin() + i); }
\ No newline at end of file
......
......@@ -14,21 +14,20 @@ class MOIP
{
public:
MOIP(void);
MOIP(const RNA& rna, const vector<Motif>& motifSites, float pthreshold, bool verbose);
MOIP(const RNA& rna, const vector<Motif>& motifSites, float theta, bool verbose);
~MOIP(void);
SecondaryStructure solve_objective(int o, double min, double max, const IloConstraintArray& F);
SecondaryStructure solve_objective(int o, double min, double max);
SecondaryStructure solve_objective(int o);
uint get_n_solutions(void) const;
const SecondaryStructure& solution(uint i) const;
void search_between(double lambdaMin, double lambdaMax, const IloConstraintArray& F_);
void search_between(double lambdaMin, double lambdaMax);
bool allowed_basepair(size_t u, size_t v) const;
void add_solution(const SecondaryStructure& s);
void remove_solution(uint i);
IloConstraintArray forbid_solutions_between(double min, double max);
void forbid_solutions_between(double min, double max);
IloEnv& get_env(void);
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_;
private:
bool is_undominated_yet(const SecondaryStructure& s);
......@@ -47,10 +46,6 @@ class MOIP
RNA rna_; // RNA object
vector<Motif> insertion_sites_; // Potential Motif insertion sites
vector<SecondaryStructure> pareto_; // Vector of results
// uint n_sets_; // number of Pareto sets to return
// Objectives related
float theta_; // theta parameter for the probability function
// CPLEX objects
IloEnv env_; // environment CPLEX object
......@@ -68,9 +63,6 @@ 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(), IloConstraintArray(env_));
}
inline SecondaryStructure MOIP::solve_objective(int o) { return solve_objective(o, 0, rna_.get_RNA_length()); }
inline IloEnv& MOIP::get_env(void) { return env_; }
#endif // MOIP_H_
\ No newline at end of file
......
......@@ -4,8 +4,6 @@
#define IL_STD
#include "rna.h"
#include <ilconcert/ilomodel.h>
#include <ilcplex/ilocplex.h>
#include <iostream>
#include <string>
#include <vector>
......@@ -49,8 +47,6 @@ class SecondaryStructure
void insert_motif(const Motif& m);
double get_objective_score(int i) const;
void set_objective_score(int i, double s);
// void set_pareto_set(uint k);
// uint get_pareto_set(void) const;
uint get_n_motifs(void) const;
uint get_n_bp(void) const;
void print(void) const;
......@@ -60,13 +56,11 @@ class SecondaryStructure
vector<double> objective_scores_; // values of the different objective functions for that SecondaryStructure
vector<pair<uint, uint>> basepairs_; // values of the decision variable of the integer program
vector<Motif> motif_info_; // information about known motives in this secondary structure and their positions
size_t n_; // length of the RNA
size_t nBP_; // number of basepairs
// uint k_; // Secondary Structure belongs to the kth Pareto set
RNA rna_; // RNA object which is folded
bool is_empty_structure; // Empty structure, returned when the solver does not find solutions anymore
IloConstraint forbid_this_; // Add it to a cplex model to forbid that solution
vector<Motif> motif_info_; // information about known motives in this secondary structure and their positions
size_t n_; // length of the RNA
size_t nBP_; // number of basepairs
RNA rna_; // RNA object which is folded
bool is_empty_structure; // Empty structure, returned when the solver does not find solutions anymore
};
// return if this SecondaryStructure s1 dominates s2
......@@ -90,7 +84,5 @@ inline double SecondaryStructure::get_objective_score(int i) const { return obje
inline void SecondaryStructure::set_objective_score(int i, double s) { objective_scores_[i - 1] = s; }
inline uint SecondaryStructure::get_n_motifs(void) const { return motif_info_.size(); }
inline uint SecondaryStructure::get_n_bp(void) const { return nBP_; }
// inline uint SecondaryStructure::get_pareto_set(void) const { return k_; }
// inline void SecondaryStructure::set_pareto_set(uint k) { k_ = k; }
#endif
......
......@@ -114,15 +114,15 @@ int main(int argc, char* argv[])
if (verbose)
cout << "\t>" << csvname << " successfuly loaded (" << posInsertionSites.size() << " insertion sites)" << endl;
/* FIND K-PARETO SETS */
/* FIND PARETO SET */
MOIP myMOIP = MOIP(myRNA, posInsertionSites, theta_p_threshold, verbose);
double min, max;
IloConstraintArray F(myMOIP.get_env());
try {
bestSSO1 = myMOIP.solve_objective(1, -__DBL_MAX__, __DBL_MAX__, F);
bestSSO2 = myMOIP.solve_objective(2, -__DBL_MAX__, __DBL_MAX__, F);
bestSSO1 = myMOIP.solve_objective(1, -__DBL_MAX__, __DBL_MAX__);
bestSSO2 = myMOIP.solve_objective(2, -__DBL_MAX__, __DBL_MAX__);
if (verbose) {
cout << endl << "Best solution according to objective 1 :" << bestSSO1.to_string() << endl;
cout << "Best solution according to objective 2 :" << bestSSO2.to_string() << endl;
......@@ -145,7 +145,7 @@ int main(int argc, char* argv[])
cout << std::setprecision(-log10(MOIP::precision_) + 4) << "\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
myMOIP.search_between(min, max);
// extend the Pareto set below
......@@ -158,37 +158,19 @@ int main(int argc, char* argv[])
min = -__DBL_MAX__;
max = bestSSO2.get_objective_score(1);
}
F.clear();
F = myMOIP.forbid_solutions_between(min, max);
if (verbose)
cout << std::setprecision(-log10(MOIP::precision_) + 4) << "\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.getSize() << " solutions found in [" << std::setprecision(10)
<< min - MOIP::precision_ << ", " << max + MOIP::precision_ << ']' << endl;
myMOIP.search_between(min, max, F);
myMOIP.search_between(min, max);
} catch (IloCplex::Exception& e) {
cerr << "Cplex Exception: " << e.getMessage() << endl;
exit(EXIT_FAILURE);
}
/* REMOVE SOLUTIONS WITH TOO HIGH LABEL */
// vector<size_t> to_remove;
// if (verbose) cout << endl;
// for (uint i = 0; i < myMOIP.get_n_solutions(); i++)
// if (myMOIP.solution(i).get_pareto_set() > 1) { // Some solution is fromm a Pareto set of too high order
// if (verbose)
// cout << "Removing structure from Pareto set " << myMOIP.solution(i).get_pareto_set() << " : "
// << myMOIP.solution(i).to_string() << endl;
// to_remove.push_back(i);
// }
// if (to_remove.size()) {
// for (size_t i = to_remove.size() - 1; i != 0; i--) myMOIP.remove_solution(to_remove[i]);
// myMOIP.remove_solution(to_remove[0]);
// }
/* DISPLAY RESULTS */
// print the pareto set
......