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