DOCloudSolver.cc 8.05 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
//=============================================================================
//
//  CLASS DOCloudSolver - IMPLEMENTATION
//
//=============================================================================

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

//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include "DOCloudSolver.hh"
11
12
//=============================================================================
#if COMISO_DOCLOUD_AVAILABLE
13
#include "DOCloudCache.hh"
14
15
#include "DOCloudJob.hh"
#include "cURLpp.hh"
16

17
#include <Base/Debug/DebUtils.hh>
18
19
20
#include <Base/Utils/OutcomeUtils.hh>

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

DEB_module("DOCloudSolver")

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

namespace COMISO {

//== IMPLEMENTATION ==========================================================
32
33
namespace DOcloud {

34
namespace {
35

36
#define P(X) ((X).data())
37
38
39
40
41
42
43
44
45
46
#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:
47
  WriteExpression(std::ostringstream& _out_str) : out_str_stream(_out_str)
48
49
50
51
52
53
  {
    start();
  }

  void start()
  {
54
    f_size_ = out_str_stream.tellp();
55
56
57
58
59
60
61
62
    at_start_ = true;
  }
  
  // Writes a monomial.
  void add_monomial(const double _coeff, const size_t _i_var)
  {
    if (_coeff == 0)
      return;
63
64
65
66
67
68
69
70
71
72
73
    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)
74
      out_str_stream << "^2";
75
    else
76
      out_str_stream << " * " << XVAR(_j_var);
77
78
79
80
81
82
83
    wrap_long_line();
  }

private:

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

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

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

122
123
124
// 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
125
std::string create_lp_string(
126
127
128
129
130
131
  NProblemInterface* _problem,
  const std::vector<NConstraintInterface*>& _constraints,
  const std::vector<PairIndexVtype>& _discrete_constraints,
  const std::vector<double>& _x
  )
{
132
133
  const int n_cols = _problem->n_unknowns(); // Unknowns #

134
  std::ostringstream lp_str_stream;
135

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

139
140
  lp_str_stream << "\\Problem name: " << std::endl << std::endl;
  lp_str_stream << "Minimize" << std::endl;
141

142
  // Writes objective function.
143
  lp_str_stream << "obj: ";
144

145
  WriteExpression wrte_expr(lp_str_stream);
146
  // 1. Linear part.
147
148
149
150
151
  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);

152
153
154
155
156
  // 2. Quadratic part (if requested).
  if (!_problem->constant_gradient())
  {
    NProblemInterface::SMatrixNP H;
    _problem->eval_hessian(P(_x), H);
157
    lp_str_stream << " + [ ";
158
159
160
161
162
    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());
    }
163
    lp_str_stream << " ] / 2";
164
165
166
  }


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

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

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

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

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

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

236
  lp_str_stream << "End";
237

238
  return std::string(lp_str_stream.str());
239
240
}

241
#undef XVAR
242

243
} // namespace
244

245
} // namespace DOcloud
246

247
void DOCloudSolver::solve(
248
249
250
  NProblemInterface*                        _problem,
  const std::vector<NConstraintInterface*>& _constraints,
  const std::vector<PairIndexVtype>&        _discrete_constraints,
251
  const double                              _time_limit
252
253
254
255
256
257
258
259
260
)
{
  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
261
262
  const std::string mip_lp = DOcloud::create_lp_string(_problem, _constraints,  
    _discrete_constraints, x);
263

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

  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));
284
285
  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),
286
287
               "DOCloudSolver solved a wrong problem.");
  )
288
289
  
  _problem->store_result(P(x));
290
291
292
293
294
295
296
297
298
299
300
}

#undef P

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

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