modeling.hpp

#pragma once

/*
 * Model a non-convex optimization problem by defining Cost and Constraint objects
 * which know how to generate a convex approximation
 *
 *
 */

#include <vector>
#include <boost/shared_ptr.hpp>
#include "sco/sco_fwd.hpp"
#include "sco/solver_interface.hpp"

namespace sco {

using std::vector;

class ConvexObjective {
public:
  ConvexObjective(Model* model) : model_(model) {}
  void addAffExpr(const AffExpr&);
  void addQuadExpr(const QuadExpr&);
  void addHinge(const AffExpr&, double coeff);
  void addAbs(const AffExpr&, double coeff);
  void addHinges(const AffExprVector&);
  void addL1Norm(const AffExprVector&);
  void addL2Norm(const AffExprVector&);
  void addMax(const AffExprVector&);
  
  bool inModel() {
    return model_ != NULL;
  }
  void addConstraintsToModel();
  void removeFromModel();
  double value(const vector<double>& x);
  
  ~ConvexObjective();

  
  Model* model_;
  QuadExpr quad_;
  vector<Var> vars_;
  vector<AffExpr> eqs_;
  vector<AffExpr> ineqs_;
  vector<Cnt> cnts_;
private:
  ConvexObjective()  {}
  ConvexObjective(ConvexObjective&)  {}
};

class ConvexConstraints {
public:
  ConvexConstraints(Model* model) : model_(model) {}
  void addEqCnt(const AffExpr&);
  void addIneqCnt(const AffExpr&);
  void setModel(Model* model) {
    assert(!inModel());
    model_ = model;
  }
  bool inModel() {
    return model_ != NULL;
  }
  void addConstraintsToModel();
  void removeFromModel();

  vector<double> violations(const vector<double>& x);
  double violation(const vector<double>& x);

  ~ConvexConstraints();
  vector<AffExpr> eqs_;
  vector<AffExpr> ineqs_;
private:
   Model* model_;
   vector<Cnt> cnts_;
   ConvexConstraints() : model_(NULL) {}
   ConvexConstraints(ConvexConstraints&) {}
};

class Cost {
public:
  virtual double value(const vector<double>&) = 0;
  virtual ConvexObjectivePtr convex(const vector<double>& x, Model* model) = 0;

  string name() {return name_;}
  void setName(const string& name) {name_=name;}
  Cost() : name_("unnamed") {}
  Cost(const string& name) : name_(name) {}
  virtual ~Cost() {}
protected:
  string name_;
};

class Constraint {
public:

  virtual ConstraintType type() = 0;
  virtual vector<double> value(const vector<double>& x) = 0;
  virtual ConvexConstraintsPtr convex(const vector<double>& x, Model* model) = 0;
  vector<double> violations(const vector<double>& x);
  double violation(const vector<double>& x);

  string name() {return name_;}
  void setName(const string& name) {name_=name;}
  Constraint() : name_("unnamed") {}
  Constraint(const string& name) : name_(name) {}
  virtual ~Constraint() {}

protected:
  string name_;
};

class EqConstraint : public Constraint{
public:
  ConstraintType type() {return EQ;}
};

class IneqConstraint : public Constraint {
public:
  ConstraintType type() {return INEQ;}
};

class OptProb {
public:
  OptProb();
  VarVector createVariables(const vector<string>& names);
  VarVector createVariables(const vector<string>& names, const vector<double>& lb, const vector<double>& ub);
  void setLowerBounds(const vector<double>& lb);
  void setUpperBounds(const vector<double>& ub);
  void setLowerBounds(const vector<double>& lb, const vector<Var>& vars);
  void setUpperBounds(const vector<double>& ub, const vector<Var>& vars);
  void addLinearConstraint(const AffExpr&, ConstraintType type);
  void addCost(CostPtr);
  void addConstraint(ConstraintPtr);
  void addEqConstraint(ConstraintPtr);
  void addIneqConstraint(ConstraintPtr);
  virtual ~OptProb() {}
  vector<double> getCentralFeasiblePoint(const vector<double>& x);
  vector<double> getClosestFeasiblePoint(const vector<double>& x);

  vector<ConstraintPtr> getConstraints() const;
  vector<CostPtr>& getCosts() {return costs_;}
  vector<ConstraintPtr>& getIneqConstraints() {return ineqcnts_;}
  vector<ConstraintPtr>& getEqConstraints() {return eqcnts_;}
  DblVec& getLowerBounds() {return lower_bounds_;}
  DblVec& getUpperBounds() {return upper_bounds_;}
  ModelPtr getModel() {return model_;}
  vector<Var>& getVars() {return vars_;}
  int getNumCosts() {return costs_.size();}
  int getNumConstraints() {return eqcnts_.size() + ineqcnts_.size();}
  int getNumVars() {return vars_.size();}

protected:
  ModelPtr model_;
  vector<Var> vars_;
  vector<double> lower_bounds_;
  vector<double> upper_bounds_;
  vector<CostPtr> costs_;
  vector<ConstraintPtr> eqcnts_;
  vector<ConstraintPtr> ineqcnts_;

  OptProb(OptProb&);
};

template <typename VecType>
inline void setVec(DblVec& x, const VarVector& vars, const VecType& vals) {
  assert(vars.size() == vals.size());
  for (int i = 0; i < vars.size(); ++i) {
    x[vars[i].var_rep->index] = vals[i];
  }
}
template <typename OutVecType>
inline OutVecType getVec1(const vector<double>& x, const VarVector& vars) {
  OutVecType out(vars.size());
  for (unsigned i=0; i < vars.size(); ++i) out[i] = x[vars[i].var_rep->index];
  return out;
}


}