opt.h

Go to the documentation of this file.

// Copyright (C) 2006 Ivo Nowak and Stefan Vigerske
// All Rights Reserved.
// This code is published under the Common Public License.
//
// Author: Ivo Nowak, Stefan Vigerske

#ifndef OPT_H
#define OPT_H

#include "standard.h"
#include "usermatrix.h"
#include "func.h"
#include "param.h"
#include "dualqqp.h"

class Solver {
      friend ostream& operator << (ostream& out, const Solver& a)
      {  a.print(out); return out; };

   protected:
      int dim_;

      double time_;

      double opt_val_;

      int iter_;
   public:
      Pointer<ostream> out_solver_p;
      Pointer<ostream> out_solver_log_p;
#define out_solver if (out_solver_p) (*out_solver_p)

#define out_solver_log if (out_solver_log_p) (*out_solver_log_p)

      double tol;

      int iter_max;

      dvector sol_point;

      Solver(int n, Pointer<ostream> out_solver_p_=out_out_p, Pointer<ostream> out_solver_log_p_=out_log_p)
      : dim_(n), sol_point(n), iter_(0), time_(-1), tol(-1.), opt_val_(INFINITY), iter_max(INF),
        out_solver_p(out_solver_p_), out_solver_log_p(out_solver_log_p_)
      {  };

      virtual ~Solver() { };

      int dim() const { return dim_;  };

      double opt_val() { return opt_val_;  };

      double time() { return time_;  };

      int iter() { return iter_;  };

      virtual int solve()=0;

      virtual int solve(dvector &x) { sol_point=x; return solve(); };

      virtual void print(ostream &out) const {
                                out << "dim: " << dim() << " tol: " << tol << endl;
                                out << "solution time: " << time_ << endl;
                                out << "iterations: " << iter_ << endl;
                        };
};


class DualSolver : public Solver {
  private:
    bool threshold_cntrl;

    bool conv_rate_cntrl;
    double stopping_rho;
    int minor_iter;
    double last_major_val;
    double last_val;
    double first_major_val;
    double max_rel_improvement;
    int improve_iter;

  protected:
    Param& param;

    DualFunc& obj;

    void do_log();

    int check(double val);

  public:
    double last_rel_improvement;

    double threshold;

    dvector lower_bound;

    Pointer<vector<double> > dual_vals;
    Pointer<vector<dvector> > dual_points;
    Pointer<vector<vector<dvector> > > orig_points;
    int log_frequency;
    bool store_primals;

    DualSolver(DualFunc& f, vector<double> lower_bound_, Param& param_, Pointer<ostream> out_solver_p_=out_out_p, Pointer<ostream> out_solver_log_p_=out_log_p)
    : Solver(f.dim(), out_solver_p_, out_solver_log_p_), threshold(INFINITY), lower_bound(f.dim()), obj(f),
      log_frequency(0), store_primals(true), param(param_), last_major_val(-INFINITY), last_val(-INFINITY), improve_iter(-1)
    { for (int i=0; i<dim(); i++) lower_bound[i]=lower_bound_[i];
      param.read();
      threshold_cntrl=param.get_i("control threshold", 1);
      conv_rate_cntrl=param.get_i("control convergence rate", 1);
                        stopping_rho=param.get_d("stopping rho", 0.1);
                        minor_iter=param_.get_i("minor iterations", 5);
                }

                virtual ~DualSolver() { }

                double rel_improvement1() {
                  return (last_major_val-first_major_val)/fabs(first_major_val);
                }

                double rel_improvement2() {
                  return last_rel_improvement;
                }

                double rel_improvement_max() {
                  return max_rel_improvement;
                }

                int serious_steps() {
                  return improve_iter;
                }

};

class LocOpt : public Solver {
        protected:
                Timer timer;
        public:
                static bool nlp_solver_available();
                static Pointer<LocOpt> get_solver(const Pointer<MinlpProblem> prob, Pointer<Param> param, char* param_prefix=NULL, Pointer<ostream> out_solver_p_=out_out_p, Pointer<ostream> out_solver_log_p_=out_log_p);

                static Pointer<LocOpt> get_lp_solver(const Pointer<MinlpProblem> prob, Pointer<Param> param, char* param_prefix=NULL, Pointer<ostream> out_solver_p_=out_out_p, Pointer<ostream> out_solver_log_p_=out_log_p);

                static Pointer<LocOpt> get_solver_origprob(const Pointer<MinlpProblem> prob, Pointer<Param> param, char* param_prefix=NULL, Pointer<ostream> out_solver_p_=out_out_p, Pointer<ostream> out_solver_log_p_=out_log_p);

                LocOpt(int n, Pointer<ostream> out_solver_p_=out_out_p, Pointer<ostream> out_solver_log_p_=out_log_p)
                : Solver(n, out_solver_p_, out_solver_log_p_)
                { }

                virtual ~LocOpt() { }

                virtual double time() { return timer; }

                virtual void reinit() { };

                virtual dvector get_lag_multipliers()=0;
};

class SimpleCut;


class MIPSolver {
        public:
                typedef enum { SOLVED, FEASIBLE, UNBOUNDED, INFEASIBLE, ITERATIONLIMITEXCEEDED, ABORTED, UNKNOWN } SolutionStatus;

                friend ostream& operator<<(ostream& out, const SolutionStatus status) {
                        switch(status) {
                                case SOLVED: out << "Solved"; break;
                                case FEASIBLE: out << "Feasible"; break;
                                case UNBOUNDED: out << "Unbounded"; break;
                                case INFEASIBLE: out << "Infeasible"; break;
                                case ITERATIONLIMITEXCEEDED: out << "Iterationlimit exceeded"; break;
                                case ABORTED: out << "Aborted"; break;
                                case UNKNOWN: out << "Unknown"; break;
                                default: out << "SolutionStatus not known!!";
                        }
                        return out;
                }

                class RowItem {
                        public:
                                RowItem() { }
                                virtual ~RowItem() { }
                };

                class ColItem {
                        public:
                                ColItem() { }
                                virtual ~ColItem() { }
                };

                static Pointer<MIPSolver> get_solver(const MipProblem& mip, Pointer<Param> param);

                virtual ~MIPSolver() { }

                virtual int nr_col()=0;
                virtual int nr_row()=0;

                virtual void set_tol(double tol)=0;
                virtual void set_maxiter(int maxiter)=0;

                virtual void reset()=0;

                virtual SolutionStatus solve()=0;

                virtual SolutionStatus solve(const UserVector<double>& x) { return solve(); };

                virtual SolutionStatus solveMIP()=0;

                virtual SolutionStatus feasible()=0;

                virtual void get_primal(UserVector<double>& x)=0;
                virtual dvector get_primal() { dvector x(nr_col()); get_primal(x); return x; }
                virtual double get_primal(const ColItem& colitem)=0;

                virtual int get_colindex(const ColItem& colitem)=0;

                virtual void get_dual(UserVector<double>& mu)=0;
                virtual double get_dual(const RowItem& rowitem)=0;

                virtual void get_reducedcosts(UserVector<double>& rc)=0;
                virtual double get_reducedcosts(const ColItem& colitem)=0;
                
                virtual void get_rowactivity(UserVector<double>& rowact)=0;
                virtual double get_rowactivity(const RowItem& rowitem)=0;
                

                virtual double get_optval()=0;

                virtual int get_iter()=0;

                virtual void set_obj(const UserVector<double>& obj, double obj_const=0.)=0;
                virtual void modify_obj(int i, double coeff)=0;

                virtual const RowItem* add_row(const UserVector<double>& row, double low, double up)=0;
                virtual const RowItem* add_row(const UserVector<double>& row, const ivector& indices, double low, double up) {
                        dvector bigrow(nr_col()); bigrow.set_block(row, indices);
                        return add_row(bigrow, low, up);
                };
                virtual void add_rows(const vector<Pointer<UserVector<double> > >& rows, const dvector& low, const dvector& up) {
                        for (int i=0; i<rows.size(); i++) add_row(*rows[i], low[i], up[i]);
                }
                virtual void add_rows(const vector<pair<dvector, ivector> >& rows, const dvector& low, const dvector& up) {
                        for (int i=0; i<rows.size(); i++) add_row(rows[i].first, rows[i].second, low[i], up[i]);
                } 
                virtual void add_rows(list<const RowItem*>& rowitems, const vector<pair<dvector, ivector> >& rows, const dvector& low, const dvector& up) {
                        for (int i=0; i<rows.size(); i++) rowitems.push_back(add_row(rows[i].first, rows[i].second, low[i], up[i]));
                }
                virtual void delete_row(const RowItem& rowitem)=0;
                virtual void delete_rows(const list<const RowItem*>& rowitems) {
                        for (list<const RowItem*>::const_iterator it(rowitems.begin()); it!=rowitems.end(); ++it) delete_row(**it);
                }
                virtual void modify_row(const RowItem& rowitem, double low, double up)=0;
                virtual void modify_row(int index, double low, double up)=0;

                virtual const ColItem* add_col(double low, double up, MipProblem::VarType=MipProblem::CONTINUOUS)=0;
                virtual const ColItem* add_col(const UserVector<double>& col, double obj_coeff, double low, double up, MipProblem::VarType=MipProblem::CONTINUOUS)=0;
                virtual void add_cols(list<const ColItem*>& colitems, const dvector& low, const dvector& up)=0;
                virtual void add_cols(list<const ColItem*>& colitems, vector<Pointer<UserVector<double> > >& cols, const vector<double>& obj_coeff, const dvector& low, const dvector& up)=0;
                virtual void delete_col(const ColItem& colitem)=0;
                virtual void delete_cols(const list<const ColItem*>& colitems) {
                        for (list<const ColItem*>::const_iterator it(colitems.begin()); it!=colitems.end(); ++it) delete_col(**it);
                }
                virtual void modify_col(const ColItem& colitem, double low, double up, MipProblem::VarType vartype)=0;
                virtual void modify_col(const ColItem& colitem, const UserVector<double>& col, double obj_coeff, double low, double up, MipProblem::VarType vartype)=0;
                virtual void modify_col(int index, double low, double up, MipProblem::VarType type)=0;
                
                virtual double get_collow(int index)=0;
                virtual double get_colup(int index)=0;
                
                virtual int generate_cuts(list<Pointer<SimpleCut> >& rowcuts)=0;
};

class LPSolver : public LocOpt {
        private:
                Pointer<MIPSolver> mipsolver;
                const Pointer<MinlpProblem> prob;
//      Pointer<Param> param;
//              char* param_prefix;

                void initmip();

                // called from solve or solve(x)
                int solved(MIPSolver::SolutionStatus);

        public:
                LPSolver(const Pointer<MinlpProblem> prob_/*, Pointer<Param> param_, char* param_prefix_=NULL*/);

                dvector get_lag_multipliers();

                void reinit() { initmip(); }

                int solve() { return solved(mipsolver->solve()); }
                int solve(dvector& start) { return solved(mipsolver->solve(start)); }
};

class BoxMinimizer : public Solver {
        protected:
          const Func& f;

          Pointer<dvector> lower, upper;

        public:
          BoxMinimizer(const Func& f_, Pointer<dvector> lower_, Pointer<dvector> upper_)
          : Solver(f_.dim()), f(f_), lower(lower_), upper(upper_)
          { }

          virtual void set_box(Pointer<dvector> lower_, Pointer<dvector> upper_) {
                lower=lower_;
                upper=upper_;
          }


};

class BoxLocOpt : public BoxMinimizer {
        private:
                Pointer<LocOpt> locopt;
                Pointer<MinlpProblem> prob;
                Pointer<Param> param;

        public:
                BoxLocOpt(Func &f_, Pointer<dvector> lower_, Pointer<dvector> upper_, Pointer<Param> param_=NULL, char* param_prefix=NULL)
                : BoxMinimizer(f_, lower_, upper_), param(param_),
                  prob(new MinlpProblem(new SepQcFunc(NULL, NULL, Pointer<Func>(&f_, false)), lower_ ? *lower_ : dvector(f_.dim()), upper_ ? *upper_ : dvector(f_.dim())))
                { locopt=LocOpt::get_solver(prob, param, param_prefix, NULL, NULL);
                }

                BoxLocOpt(const SepQcFunc& f_, Pointer<dvector> lower_, Pointer<dvector> upper_, Pointer<Param> param_=NULL, char* param_prefix=NULL)
                : BoxMinimizer(f_, lower_, upper_), param(param_),
                        prob(new MinlpProblem(Pointer<SepQcFunc>((SepQcFunc*)&f_, false), lower_ ? *lower_ : dvector(f_.dim()), upper_ ? *upper_ : dvector(f_.dim())))
                { locopt=LocOpt::get_solver(prob, param, param_prefix, NULL, NULL);
                }

                BoxLocOpt(SepQcFunc& f_, Pointer<dvector> lower_, Pointer<dvector> upper_, Pointer<Param> param_=NULL, char* param_prefix=NULL)
                : BoxMinimizer(f_, lower_, upper_), param(param_),
                        prob(new MinlpProblem(Pointer<SepQcFunc>(&f_, false), lower_ ? *lower_ : dvector(f_.dim()), upper_ ? *upper_ : dvector(f_.dim())))
                { locopt=LocOpt::get_solver(prob, param, param_prefix, NULL, NULL);
                }

                void set_box(Pointer<dvector> lower_, Pointer<dvector> upper_) {
                        if (lower==lower_ && upper==upper_) return; // same bounds, so no change needed
      BoxMinimizer::set_box(lower_, upper_);
      prob->lower=*lower;
      prob->upper=*upper;
      locopt->iter_max=iter_max;
      locopt->reinit();
                }

                int solve() {
                        locopt->iter_max=iter_max;
                        int ret=locopt->solve();
                        sol_point=locopt->sol_point;
                        opt_val_=locopt->opt_val();
                        return ret;
                }

                int solve(dvector& x) {
                        locopt->iter_max=iter_max;
                        int ret=locopt->solve(x);
                        sol_point=locopt->sol_point;
                        opt_val_=locopt->opt_val();
                        return ret;
                }

};

#endif // OPT_H

---