CYCLUS
prog_translator.cc
Go to the documentation of this file.
1 #include "prog_translator.h"
2 
3 #include <algorithm>
4 
5 #include "CoinPackedVector.hpp"
6 #include "OsiSolverInterface.hpp"
7 
8 #include "cyc_limits.h"
9 #include "error.h"
10 #include "exchange_graph.h"
11 #include "exchange_solver.h"
12 #include "logger.h"
13 
14 namespace cyclus {
15 
16 ProgTranslator::ProgTranslator(ExchangeGraph* g, OsiSolverInterface* iface)
17  : g_(g),
18  iface_(iface),
19  excl_(false),
20  pseudo_cost_(std::numeric_limits<double>::max()) {
21  Init();
22 }
23 
24 ProgTranslator::ProgTranslator(ExchangeGraph* g, OsiSolverInterface* iface,
25  bool exclusive)
26  : g_(g),
27  iface_(iface),
28  excl_(exclusive),
29  pseudo_cost_(std::numeric_limits<double>::max()) {
30  Init();
31 }
32 
33 ProgTranslator::ProgTranslator(ExchangeGraph* g, OsiSolverInterface* iface,
34  double pseudo_cost)
35  : g_(g),
36  iface_(iface),
37  excl_(false),
38  pseudo_cost_(pseudo_cost) {
39  Init();
40 }
41 
42 ProgTranslator::ProgTranslator(ExchangeGraph* g, OsiSolverInterface* iface,
43  bool exclusive, double pseudo_cost)
44  : g_(g),
45  iface_(iface),
46  excl_(exclusive),
47  pseudo_cost_(pseudo_cost) {
48  Init();
49 }
50 
51 void ProgTranslator::Init() {
52  arc_offset_ = g_->arcs().size();
53  int n_cols = arc_offset_ + g_->request_groups().size();
54  ctx_.obj_coeffs.resize(n_cols);
55  ctx_.col_ubs.resize(n_cols);
56  ctx_.col_lbs.resize(n_cols);
57  ctx_.m = CoinPackedMatrix(false, 0, 0);
58 }
59 
60 void ProgTranslator::CheckPref(double pref) {
61  if (pref <= 0) {
62  std::stringstream ss;
63  ss << "Preference value found to be nonpositive (" << pref
64  << "). Preferences must be positive when using an optimization solver."
65  << " If using Cyclus in simulation mode (e.g., from the command line),"
66  << " this error is likely a bug in Cyclus. Please report it to the developer's "
67  << "list (https://groups.google.com/forum/#!forum/cyclus-dev).";
68  throw ValueError(ss.str());
69  }
70 }
71 
73  // number of variables = number of arcs + 1 faux arc per request group with arcs
74  int nfalse = 0;
75  std::vector<RequestGroup::Ptr>& rgs = g_->request_groups();
76  for (int i = 0; i != g_->request_groups().size(); ++i)
77  nfalse += rgs[i].get()->HasArcs() ? 1 : 0;
78  int n_cols = g_->arcs().size() + nfalse;
79  ctx_.m.setDimensions(0, n_cols);
80 
81  bool request;
82  std::vector<ExchangeNodeGroup::Ptr>& sgs = g_->supply_groups();
83  for (int i = 0; i != sgs.size(); i++) {
84  request = false;
85  XlateGrp_(sgs[i].get(), request);
86  }
87 
88  for (int i = 0; i != rgs.size(); i++) {
89  request = true;
90  XlateGrp_(rgs[i].get(), request);
91  }
92 
93  // add each false arc
94  CLOG(LEV_DEBUG1) << "Adding " << arc_offset_ - g_->arcs().size()
95  << " false arcs.";
96  double inf = iface_->getInfinity();
97  for (int i = g_->arcs().size(); i != arc_offset_; i++) {
98  ctx_.obj_coeffs[i] = pseudo_cost_;
99  ctx_.col_lbs[i] = 0;
100  ctx_.col_ubs[i] = inf;
101  }
102 }
103 
105  iface_->setObjSense(1.0); // minimize
106 
107  // load er up!
108  iface_->loadProblem(ctx_.m, &ctx_.col_lbs[0], &ctx_.col_ubs[0],
109  &ctx_.obj_coeffs[0], &ctx_.row_lbs[0], &ctx_.row_ubs[0]);
110 
111 
112  if (excl_) {
113  std::vector<Arc>& arcs = g_->arcs();
114  for (int i = 0; i != arcs.size(); i++) {
115  Arc& a = arcs[i];
116  if (a.exclusive()) {
117  iface_->setInteger(g_->arc_ids()[a]);
118  }
119  }
120  }
121 
122 }
123 
125  Translate();
126  Populate();
127 }
128 
129 void ProgTranslator::XlateGrp_(ExchangeNodeGroup* grp, bool request) {
130  double inf = iface_->getInfinity();
131  std::vector<double>& caps = grp->capacities();
132 
133  if (request && !grp->HasArcs())
134  return; // no arcs, no reason to add variables/constraints
135 
136  std::vector<CoinPackedVector> cap_rows;
137  std::vector<CoinPackedVector> excl_rows;
138  for (int i = 0; i != caps.size(); i++) {
139  cap_rows.push_back(CoinPackedVector());
140  }
141 
142  std::vector<ExchangeNode::Ptr>& nodes = grp->nodes();
143  for (int i = 0; i != nodes.size(); i++) {
144  std::map<Arc, std::vector<double> >& ucap_map = nodes[i]->unit_capacities;
145  std::map<Arc, std::vector<double> >::iterator cap_it;
146 
147  // add each arc
148  for (cap_it = ucap_map.begin(); cap_it != ucap_map.end(); ++cap_it) {
149  const Arc& a = cap_it->first;
150  std::vector<double>& ucaps = cap_it->second;
151  int arc_id = g_->arc_ids()[a];
152 
153  // add each unit capacity coefficient
154  for (int j = 0; j != ucaps.size(); j++) {
155  double coeff = ucaps[j];
156  if (excl_ && a.exclusive()) {
157  coeff *= a.excl_val();
158  }
159 
160  cap_rows[j].insert(arc_id, coeff);
161  }
162 
163  if (request) {
164  CheckPref(a.pref());
165  ctx_.obj_coeffs[arc_id] = ExchangeSolver::Cost(a, excl_);
166  ctx_.col_lbs[arc_id] = 0;
167  ctx_.col_ubs[arc_id] = (excl_ && a.exclusive()) ? 1 :
168  std::min(nodes[i]->qty, inf);
169  }
170  }
171  }
172 
173  int faux_id;
174  if (request) {
175  faux_id = arc_offset_++;
176  }
177 
178  // add all capacity rows
179  for (int i = 0; i != cap_rows.size(); i++) {
180  if (request) {
181  cap_rows[i].insert(faux_id, 1.0); // faux arc
182  }
183 
184  // 1e15 is the largest value that doesn't make the solver fall over
185  // (by emperical testing)
186  double rlb = std::min(caps[i], 1e15);
187  ctx_.row_lbs.push_back(request ? rlb : 0);
188  ctx_.row_ubs.push_back(request ? inf : caps[i]);
189  ctx_.m.appendRow(cap_rows[i]);
190  }
191 
192  if (excl_) {
193  // add exclusive arcs
194  std::vector< std::vector<ExchangeNode::Ptr> >& exngs =
195  grp->excl_node_groups();
196  for (int i = 0; i != exngs.size(); i++) {
197  CoinPackedVector excl_row;
198  std::vector<ExchangeNode::Ptr>& nodes = exngs[i];
199  for (int j = 0; j != nodes.size(); j++) {
200  std::vector<Arc>& arcs = g_->node_arc_map()[nodes[j]];
201  for (int k = 0; k != arcs.size(); k++) {
202  excl_row.insert(g_->arc_ids()[arcs[k]], 1.0);
203  }
204  }
205  if (excl_row.getNumElements() > 0) {
206  excl_rows.push_back(excl_row);
207  }
208  }
209 
210  // add all exclusive rows
211  for (int i = 0; i != excl_rows.size(); i++) {
212  ctx_.row_lbs.push_back(0.0);
213  ctx_.row_ubs.push_back(1.0);
214  ctx_.m.appendRow(excl_rows[i]);
215  }
216  }
217 }
218 
220  const double* sol = iface_->getColSolution();
221  std::vector<Arc>& arcs = g_->arcs();
222  double flow;
223  for (int i = 0; i < arcs.size(); i++) {
224  Arc& a = g_->arc_by_id().at(i);
225  flow = sol[i];
226  flow = (excl_ && a.exclusive()) ? flow * a.excl_val() : flow;
227  if (flow > cyclus::eps()) {
228  g_->AddMatch(a, flow);
229  }
230  }
231 }
232 
234  throw DepricationError("Class ProgTranslator::Context is now deprecated "
235  "in favor of ProgTranslatorContext.");
236 }
237 
239  throw DepricationError("Class ProgTranslator::Context is now deprecated "
240  "in favor of ProgTranslatorContext.");
241 }
242 
243 
244 } // namespace cyclus
const std::vector< std::vector< ExchangeNode::Ptr > > & excl_node_groups() const
exclusive node groups represent nodes over whose combined arcs flow can only exist on one arc ...
void Populate()
populates the solver interface with values from the translators Context
An arc represents a possible connection between two nodes in the bipartite resource exchange graph...
bool exclusive() const
For values that are too big, too small, etc.
Definition: error.h:41
const std::vector< ExchangeNode::Ptr > & nodes() const
const std::vector< Arc > & arcs() const
A ExchangeNodeGroup is a collection of ExchangeNodes, and is the ExchangeGraph representation of a Bi...
const std::vector< double > & capacities() const
the flow capacities assocaited with this group
#define CLOG(level)
Definition: logger.h:39
void Translate()
translates the graph, filling the translators Context
void ToProg()
translates graph into mathematic program via iface.
ProgTranslator(ExchangeGraph *g, OsiSolverInterface *iface)
constructor
const double g
Definition: material.h:18
const std::vector< ExchangeNodeGroup::Ptr > & supply_groups() const
void AddMatch(const Arc &a, double qty)
adds a match for a quanity of flow along an arc
For depricating API until a next major release.
Definition: error.h:77
const std::map< Arc, int > & arc_ids() const
const std::map< ExchangeNode::Ptr, std::vector< Arc > > & node_arc_map() const
Code providing rudimentary logging capability for the Cyclus core.
taken directly from OsiSolverInterface.cpp on 2/17/14 from https://projects.coin-or.org/Osi/browser/trunk.
Definition: agent.cc:14
const std::map< int, Arc > & arc_by_id() const
double pref() const
debugging information - least verbose
Definition: logger.h:58
std::vector< double > row_lbs
std::vector< double > col_ubs
double excl_val() const
void FromProg()
translates solution from iface back into graph matches
const std::vector< RequestGroup::Ptr > & request_groups() const
std::vector< double > obj_coeffs
An ExchangeGraph is a resource-neutral representation of a ResourceExchange.
std::vector< double > row_ubs
double eps()
a generic epsilon value
Definition: cyc_limits.h:12
std::vector< double > col_lbs
static double Cost(const Arc &a, bool exclusive_orders=kDefaultExclusive)
return the cost of an arc