DOCloudSolver.cc 8.09 KB
Newer Older
1 2
// (C) Copyright 2015 by Autodesk, Inc.

3 4 5 6 7 8 9 10 11 12
//=============================================================================
//
//  CLASS DOCloudSolver - IMPLEMENTATION
//
//=============================================================================

//== INCLUDES =================================================================

//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include "DOCloudSolver.hh"
13 14
//=============================================================================
#if COMISO_DOCLOUD_AVAILABLE
15
#include "DOCloudCache.hh"
16 17
#include "DOCloudJob.hh"
#include "cURLpp.hh"
18
#include "CoMISo/Utils/CoMISoError.hh"
19

20
#include <Base/Debug/DebUtils.hh>
21 22

#include <stdexcept>
23
#include <algorithm>
24 25
#include <fstream>
#include <iomanip>
26 27 28 29 30 31 32 33

DEB_module("DOCloudSolver")

//== NAMESPACES ===============================================================

namespace COMISO {

//== IMPLEMENTATION ==========================================================
34 35
namespace DOcloud {

36
namespace {
37

38
#define P(X) ((X).data())
39 40 41 42 43 44 45 46 47 48
#define XVAR(IDX) "x" << IDX

class WriteExpression
{
  // lp format allows a max line length = 560.
  // For simplicity we put end line when the number of written characters on the
  // same line exceeds LINE_TRESHOLD_LEN.
  enum { LINE_TRESHOLD_LEN = 100 };

public:
49
  WriteExpression(std::ostringstream& _out_str) : out_str_stream(_out_str)
50 51 52 53 54 55
  {
    start();
  }

  void start()
  {
56
    f_size_ = out_str_stream.tellp();
57 58 59 60 61 62 63 64
    at_start_ = true;
  }
  
  // Writes a monomial.
  void add_monomial(const double _coeff, const size_t _i_var)
  {
    if (_coeff == 0)
      return;
65 66 67 68 69 70 71 72 73 74 75
    add_monomial_internal(_coeff, _i_var);
    wrap_long_line();
  }

  // Writes a binomial.
  void add_binomial(const double _coeff, const size_t _i_var, const size_t _j_var)
  {
    if (_coeff == 0)
      return;
    add_monomial_internal(_coeff, _i_var);
    if (_j_var == _i_var)
76
      out_str_stream << "^2";
77
    else
78
      out_str_stream << " * " << XVAR(_j_var);
79 80 81 82 83 84 85
    wrap_long_line();
  }

private:

  void wrap_long_line()
  {
86
    const auto new_f_size = out_str_stream.tellp();
87 88
    if (new_f_size - f_size_ > LINE_TRESHOLD_LEN)
    {
89
      out_str_stream << std::endl;
90 91 92 93 94 95
      f_size_ = new_f_size;
    }
  }

  void add_monomial_internal(const double _coeff, const size_t _i_var)
  {
96 97 98
    if (_coeff == 1)
    {
      if (!at_start_)
99
        out_str_stream << " + ";
100 101
    }
    else if (_coeff == -1)
102
      out_str_stream << " - ";
103 104 105 106 107
    else
    {
      if (!at_start_)
      {
        if (_coeff > 0)
108
          out_str_stream << " + ";
109
        else
110
          out_str_stream << ' ';
111
      }
112
      out_str_stream << _coeff << ' ';
113
    }
114
    out_str_stream << XVAR(_i_var);
115 116 117 118
    at_start_ = false;
  }

private:
119
  std::ostringstream& out_str_stream;
120 121 122
  std::fstream::pos_type f_size_;
  bool at_start_;
};
123

124 125 126
// Create a lp file for the given constraints and object function.
// Here is the lp format specifications:
// http://www-01.ibm.com/support/knowledgecenter/SSSA5P_12.6.1/ilog.odms.cplex.help/CPLEX/FileFormats/topics/LP.html
127
std::string create_lp_string(
128 129 130 131 132 133
  NProblemInterface* _problem,
  const std::vector<NConstraintInterface*>& _constraints,
  const std::vector<PairIndexVtype>& _discrete_constraints,
  const std::vector<double>& _x
  )
{
134 135
  const int n_cols = _problem->n_unknowns(); // Unknowns #

136
  std::ostringstream lp_str_stream;
137

138
  // Set the ofstream options.
139
  lp_str_stream << std::setprecision(std::numeric_limits<double>::digits10 + 2);
140

141 142
  lp_str_stream << "\\Problem name: " << std::endl << std::endl;
  lp_str_stream << "Minimize" << std::endl;
143

144
  // Writes objective function.
145
  lp_str_stream << "obj: ";
146

147
  WriteExpression wrte_expr(lp_str_stream);
148
  // 1. Linear part.
149 150 151 152 153
  std::vector<double> objective(n_cols);
  _problem->eval_gradient(P(_x), P(objective));
  for (size_t i = 0; i < objective.size(); ++i)
    wrte_expr.add_monomial(objective[i], i);

154 155 156 157 158
  // 2. Quadratic part (if requested).
  if (!_problem->constant_gradient())
  {
    NProblemInterface::SMatrixNP H;
    _problem->eval_hessian(P(_x), H);
159
    lp_str_stream << " + [ ";
160 161 162 163 164
    for (int i = 0; i < H.outerSize(); ++i)
    {
      for (NProblemInterface::SMatrixNP::InnerIterator it(H, i); it; ++it)
        wrte_expr.add_binomial(it.value(), it.row(), it.col());
    }
165
    lp_str_stream << " ] / 2";
166 167 168
  }


169
  // Writes constraints.
170
  lp_str_stream << std::endl << "Subject To" << std::endl;
171
  for (const auto& cstr : _constraints)
172 173
  {
    NConstraintInterface::SVectorNC gc;
174
    cstr->eval_gradient(P(_x), gc);
175

176
    wrte_expr.start();
177 178
    for (NConstraintInterface::SVectorNC::InnerIterator v_it(gc); v_it; ++v_it)
    {
179 180 181 182 183 184
      auto coeff = v_it.value();
      wrte_expr.add_monomial(coeff, v_it.index());
    }
    switch (cstr->constraint_type())
    {
    case NConstraintInterface::NC_EQUAL:
185
      lp_str_stream << " = ";
186
      break;
187
    case NConstraintInterface::NC_GREATER_EQUAL:
188
      lp_str_stream << " >= ";
189
      break;
190
    case NConstraintInterface::NC_LESS_EQUAL:
191
      lp_str_stream << " <= ";
192 193
      break;
    }
194
    lp_str_stream << -cstr->eval_constraint(P(_x)) << std::endl;
195 196
  }

197
  // Writes the variables.
198
  lp_str_stream << "Bounds" << std::endl;
199
  for (size_t i = 0; i < n_cols; ++i)
200
    lp_str_stream << XVAR(i) << " Free" << std::endl;
201

202 203 204
  // Integer and binary variables.
  std::vector<unsigned int> int_var, bin_var;
  for (const auto& dc : _discrete_constraints)
205
  {
206 207 208 209 210
    if (dc.second == Integer)
      int_var.push_back(dc.first);
    else if (dc.second == Binary)
      bin_var.push_back(dc.first);
  }
211
  auto write_var_set = [&lp_str_stream](const std::vector<unsigned int>& _vars,
212 213 214 215 216
    const char* _type)
  {
    if (_vars.empty())
      return;
    // Writes integer variables.
217
    lp_str_stream << _type << std::endl;
218
    auto var_it = _vars.begin();
219
    lp_str_stream << XVAR(*var_it);
220 221
    size_t n_wrt_var = 1;
    while (++var_it != _vars.end())
222
    {
223
      if (n_wrt_var++ % 16) // 16 variables per line. Lines length must be < 560.
224
        lp_str_stream << ' ';
225
      else
226 227
        lp_str_stream << std::endl;
      lp_str_stream << XVAR(*var_it);
228
    }
229
    lp_str_stream << std::endl;
230 231 232 233 234 235 236 237

  };
  // Writes integer variables.
  write_var_set(int_var, "Integers");

  // Writes Binary variables.
  write_var_set(bin_var, "Binary");

238
  lp_str_stream << "End";
239

240
  return std::string(lp_str_stream.str());
241 242
}

243
#undef XVAR
244

245
} // namespace
246

247
} // namespace DOcloud
248

249
void DOCloudSolver::solve(
250 251 252
  NProblemInterface*                        _problem,
  const std::vector<NConstraintInterface*>& _constraints,
  const std::vector<PairIndexVtype>&        _discrete_constraints,
253
  const double                              _time_limit
254 255 256 257 258 259 260 261 262
)
{
  DEB_enter_func;
  DEB_warning_if(!_problem->constant_hessian(), 1,
    "DOCloudSolver received a problem with non-constant hessian!");
  DEB_warning_if(!_problem->constant_gradient(), 1,
    "DOCloudSolver received a problem with non-constant gradient!");

  std::vector<double> x(_problem->n_unknowns(), 0.0); // solution
263 264
  const std::string mip_lp = DOcloud::create_lp_string(_problem, _constraints,  
    _discrete_constraints, x);
265

266
  double obj_val;
267 268
  DOcloud::Cache cache(mip_lp);
  if (cache.restore_result(x, obj_val))
269
    DEB_line(3, "MIP cached.")
270 271
  else
  {
272
    DEB_line(1, "MIP not cached, computing optimization.");
273 274
    const std::string lp_hash = cache.hash() + ".lp";
    DOcloud::Job job(lp_hash, mip_lp);
275 276 277
    job.setup();
    job.wait();
    obj_val = job.solution(x);
278
    cache.store_result(x, obj_val);
279
  }
280
  COMISO_THROW_if(x.empty(), MIPS_NO_SOLUTION);
281 282 283 284 285

  DEB_only(
  // The lp problem ignores the constant term in the objective function.
  const std::vector<double> x_zero(_problem->n_unknowns(), 0.0);
  const auto zero_val = _problem->eval_f(P(x_zero));
286 287
  const auto test_obj_val = _problem->eval_f(P(x));
  DEB_error_if(fabs(obj_val + zero_val - test_obj_val) > (1e-8 + zero_val * 0.01),
288 289
               "DOCloudSolver solved a wrong problem.");
  )
290 291
  
  _problem->store_result(P(x));
292 293 294 295 296 297 298 299 300 301 302
}

#undef P

//=============================================================================
} // namespace COMISO
//=============================================================================

#endif // COMISO_DOCLOUD_AVAILABLE
//=============================================================================