CYCLUS
greedy_solver.cc
Go to the documentation of this file.
1 #include "greedy_solver.h"
2 
3 #include <algorithm>
4 #include <cassert>
5 #include <functional>
6 #include <vector>
7 #include <boost/math/special_functions/next.hpp>
8 
9 #include "cyc_limits.h"
10 #include "error.h"
11 #include "logger.h"
12 
13 namespace cyclus {
14 
15 void Capacity(cyclus::Arc const&, double, double) {};
16 void Capacity(boost::shared_ptr<cyclus::ExchangeNode>, cyclus::Arc const&,
17  double) {};
18 
20  : conditioner_(c),
21  ExchangeSolver(exclusive_orders) {}
22 
23 GreedySolver::GreedySolver(bool exclusive_orders)
24  : ExchangeSolver(exclusive_orders) {
25  conditioner_ = new cyclus::GreedyPreconditioner();
26 }
27 
29  : conditioner_(c),
30  ExchangeSolver(true) {}
31 
33  conditioner_ = new cyclus::GreedyPreconditioner();
34 }
35 
37  if (conditioner_ != NULL)
38  delete conditioner_;
39 }
40 
42  if (conditioner_ != NULL)
43  conditioner_->Condition(graph_);
44 }
45 
47  std::for_each(graph_->request_groups().begin(),
48  graph_->request_groups().end(),
49  std::bind1st(
50  std::mem_fun(&GreedySolver::GetCaps),
51  this));
52 
53  std::for_each(graph_->supply_groups().begin(),
54  graph_->supply_groups().end(),
55  std::bind1st(
56  std::mem_fun(&GreedySolver::GetCaps),
57  this));
58 }
59 
61  double pseudo_cost = PseudoCost(); // from ExchangeSolver API
62  Condition();
63  obj_ = 0;
64  unmatched_ = 0;
65  n_qty_.clear();
66 
67  Init();
68 
69  std::for_each(graph_->request_groups().begin(),
70  graph_->request_groups().end(),
71  std::bind1st(
72  std::mem_fun(&GreedySolver::GreedilySatisfySet),
73  this));
74 
75  obj_ += unmatched_ * pseudo_cost;
76  return obj_;
77 }
78 
79 double GreedySolver::Capacity(const Arc& a, double u_curr_qty,
80  double v_curr_qty) {
81  bool min = true;
82  double ucap = Capacity(a.unode(), a, !min, u_curr_qty);
83  double vcap = Capacity(a.vnode(), a, min, v_curr_qty);
84 
85  CLOG(cyclus::LEV_DEBUG1) << "Capacity for unode of arc: " << ucap;
86  CLOG(cyclus::LEV_DEBUG1) << "Capacity for vnode of arc: " << vcap;
87  CLOG(cyclus::LEV_DEBUG1) << "Capacity for arc : "
88  << std::min(ucap, vcap);
89 
90  return std::min(ucap, vcap);
91 }
92 
93 double GreedySolver::Capacity(ExchangeNode::Ptr n, const Arc& a, bool min_cap,
94  double curr_qty) {
95  if (n->group == NULL) {
96  throw cyclus::StateError("An notion of node capacity requires a nodegroup.");
97  }
98 
99  if (n->unit_capacities[a].size() == 0) {
100  return n->qty - curr_qty;
101  }
102 
103  std::vector<double>& unit_caps = n->unit_capacities[a];
104  const std::vector<double>& group_caps = grp_caps_[n->group];
105  std::vector<double> caps;
106  double grp_cap, u_cap, cap;
107 
108  for (int i = 0; i < unit_caps.size(); i++) {
109  grp_cap = group_caps[i];
110  u_cap = unit_caps[i];
111  cap = grp_cap / u_cap;
112  CLOG(cyclus::LEV_DEBUG1) << "Capacity for node: ";
113  CLOG(cyclus::LEV_DEBUG1) << " group capacity: " << grp_cap;
114  CLOG(cyclus::LEV_DEBUG1) << " unit capacity: " << u_cap;
115  CLOG(cyclus::LEV_DEBUG1) << " capacity: " << cap;
116 
117  // special case for unlimited capacities
118  if (grp_cap == std::numeric_limits<double>::max()) {
119  caps.push_back(std::numeric_limits<double>::max());
120  } else {
121  caps.push_back(cap);
122  }
123  }
124 
125  if (min_cap) { // the smallest value is constraining (for bids)
126  cap = *std::min_element(caps.begin(), caps.end());
127  } else { // the largest value must be met (for requests)
128  cap = *std::max_element(caps.begin(), caps.end());
129  }
130  return std::min(cap, n->qty - curr_qty);
131 }
132 
133 void GreedySolver::GetCaps(ExchangeNodeGroup::Ptr g) {
134  for (int i = 0; i != g->nodes().size(); i++) {
135  n_qty_[g->nodes()[i]] = 0;
136  }
137  grp_caps_[g.get()] = g->capacities();
138 }
139 
140 void GreedySolver::GreedilySatisfySet(RequestGroup::Ptr prs) {
141  std::vector<ExchangeNode::Ptr>& nodes = prs->nodes();
142  std::stable_sort(nodes.begin(), nodes.end(), AvgPrefComp);
143 
144  std::vector<ExchangeNode::Ptr>::iterator req_it = nodes.begin();
145  double target = prs->qty();
146  double match = 0;
147 
148  ExchangeNode::Ptr u, v;
149  std::vector<Arc>::const_iterator arc_it;
150  std::vector<Arc> sorted;
151  double remain, tomatch, excl_val;
152 
153  CLOG(LEV_DEBUG1) << "Greedy Solving for " << target
154  << " amount of a resource.";
155 
156  while ((match <= target) && (req_it != nodes.end())) {
157  // this if statement is needed because map.at() will throw if the key does
158  // not exist, which is a corner case for when there is a request with no bid
159  // arcs associated with it
160  if (graph_->node_arc_map().count(*req_it) > 0) {
161  const std::vector<Arc>& arcs = graph_->node_arc_map().at(*req_it);
162  sorted = std::vector<Arc>(arcs); // make a copy for now
163  std::stable_sort(sorted.begin(), sorted.end(), ReqPrefComp);
164  arc_it = sorted.begin();
165 
166  while ((match <= target) && (arc_it != sorted.end())) {
167  remain = target - match;
168  const Arc& a = *arc_it;
169  u = a.unode();
170  v = a.vnode();
171  // capacity adjustment
172  tomatch = std::min(remain, Capacity(a, n_qty_[u], n_qty_[v]));
173 
174  // exclusivity adjustment
175  if (arc_it->exclusive()) {
176  excl_val = a.excl_val();
177 
178  // this careful float comparison is vital for preventing false positive
179  // constraint violations w.r.t. exclusivity-related capacity.
180  double dist = boost::math::float_distance(tomatch, excl_val);
181  if (dist >= float_ulp_eq ) {
182  tomatch = 0;
183  } else {
184  tomatch = excl_val;
185  }
186  }
187 
188  if (tomatch > eps()) {
189  CLOG(LEV_DEBUG1) << "Greedy Solver is matching " << tomatch
190  << " amount of a resource.";
191  UpdateCapacity(u, a, tomatch);
192  UpdateCapacity(v, a, tomatch);
193  n_qty_[u] += tomatch;
194  n_qty_[v] += tomatch;
195  graph_->AddMatch(a, tomatch);
196 
197  match += tomatch;
198  UpdateObj(tomatch, u->prefs[a]);
199  }
200  ++arc_it;
201  } // while( (match =< target) && (arc_it != arcs.end()) )
202  } // if(graph_->node_arc_map().count(*req_it) > 0)
203  ++req_it;
204  } // while( (match =< target) && (req_it != nodes.end()) )
205 
206  unmatched_ += target - match;
207 }
208 
209 void GreedySolver::UpdateObj(double qty, double pref) {
210  // updates minimizing object (i.e., 1/pref is a cost and the objective is cost
211  // * flow)
212  obj_ += qty / pref;
213 }
214 
215 void GreedySolver::UpdateCapacity(ExchangeNode::Ptr n, const Arc& a,
216  double qty) {
217  using cyclus::IsNegative;
218  using cyclus::ValueError;
219 
220  std::vector<double>& unit_caps = n->unit_capacities[a];
221  std::vector<double>& caps = grp_caps_[n->group];
222  assert(unit_caps.size() == caps.size());
223  for (int i = 0; i < caps.size(); i++) {
224  double prev = caps[i];
225  // special case for unlimited capacities
226  CLOG(cyclus::LEV_DEBUG1) << "Updating capacity value from: "
227  << prev;
228  caps[i] = (prev == std::numeric_limits<double>::max()) ?
229  std::numeric_limits<double>::max() :
230  prev - qty * unit_caps[i];
231  CLOG(cyclus::LEV_DEBUG1) << " to: "
232  << caps[i];
233  }
234 
235  if (IsNegative(n->qty - qty)) {
236  std::stringstream ss;
237  ss << "A bid for " << n->commod << " was set at " << n->qty
238  << " but has been matched to a higher value " << qty
239  << ". This could be due to a problem with your "
240  << "bid portfolio constraints.";
241  throw ValueError(ss.str());
242  }
243 }
244 
245 } // namespace cyclus
bool AvgPrefComp(ExchangeNode::Ptr l, ExchangeNode::Ptr r)
A comparison function for sorting a container of Nodes by the nodes preference in decensing order (i...
Definition: greedy_solver.h:39
double Capacity(const Arc &a, double u_curr_qty, double v_curr_qty)
the capacity of the arc
void Condition()
Uses the provided (or a default) GreedyPreconditioner to condition the solver&#39;s ExchangeGraph so that...
An arc represents a possible connection between two nodes in the bipartite resource exchange graph...
For failed object state expectations.
Definition: error.h:53
For values that are too big, too small, etc.
Definition: error.h:41
boost::shared_ptr< RequestGroup > Ptr
bool ReqPrefComp(const Arc &l, const Arc &r)
double Capacity(const Arc& a, double u_curr_qty, double v_curr_qty) { return 0; } ...
Definition: greedy_solver.h:26
boost::shared_ptr< ExchangeNode > Ptr
void Init()
Initialize member values based on the given graph.
#define CLOG(level)
Definition: logger.h:39
boost::shared_ptr< ExchangeNodeGroup > Ptr
ExchangeGraph * graph_
boost::shared_ptr< ExchangeNode > vnode() const
static const double float_ulp_eq
distance in ULP within which floating point numbers should be considered equal.
Definition: cyc_limits.h:35
virtual double SolveGraph()
the GreedySolver solves an ExchangeGraph by iterating over each RequestGroup and matching requests wi...
const double g
Definition: material.h:18
bool IsNegative(double d)
returns true if a double is less than 0 - eps()
Definition: cyc_limits.h:24
A GreedyPreconditioner conditions an ExchangeGraph for a GreedySolver by ordering the RequestGroups a...
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
void Capacity(cyclus::Arc const &, double, double)
GreedySolver()
GreedySolver constructor.
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
debugging information - least verbose
Definition: logger.h:58
double excl_val() const
a very simple interface for solving translated resource exchanges
const std::vector< RequestGroup::Ptr > & request_groups() const
boost::shared_ptr< ExchangeNode > unode() const
void Condition(ExchangeGraph *graph)
conditions the graph as described above
double PseudoCost()
Calculates the ratio of the maximum objective coefficient to minimum unit capacity plus an added cost...
double eps()
a generic epsilon value
Definition: cyc_limits.h:12