22921176376

Committed 10 Mar 2026 07:44PM UTC coverage: 36.502%. First build

Build # 22921176376

Build Type

Pull #1937

github

Committed by

web-flow

Commit Message

Merge c90a08f4d into 32dbe028c

Pull Request Pull Request #1937: Changing the DRE's Objective function to Surplus Maximization

Coverage Stats

249 of 283 new or added lines in 26 files covered. (87.99%)

52424 of 143618 relevant lines covered (36.5%)

74129.23 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

86.78

/src/greedy_solver.cc

#include "greedy_solver.h"

#include <algorithm>
#include <cassert>
#include <functional>
#include <vector>
#include <boost/math/special_functions/next.hpp>

#include "cyc_limits.h"
#include "error.h"
#include "logger.h"

namespace cyclus {

void Capacity(cyclus::Arc const&, double, double) {};
void Capacity(boost::shared_ptr<cyclus::ExchangeNode>, cyclus::Arc const&,
              double) {};

GreedySolver::GreedySolver(bool exclusive_orders, GreedyPreconditioner* c)
    : conditioner_(c), ExchangeSolver(exclusive_orders), shift_(0.0) {}

GreedySolver::GreedySolver(bool exclusive_orders)
    : ExchangeSolver(exclusive_orders), shift_(0.0) {
  conditioner_ = new cyclus::GreedyPreconditioner();
}

GreedySolver::GreedySolver(GreedyPreconditioner* c)
    : conditioner_(c), ExchangeSolver(true), shift_(0.0) {}

GreedySolver::GreedySolver() : ExchangeSolver(true), shift_(0.0) {
  conditioner_ = new cyclus::GreedyPreconditioner();
}

GreedySolver::~GreedySolver() {
  if (conditioner_ != NULL) delete conditioner_;
}

void GreedySolver::Condition() {
  if (conditioner_ != NULL) conditioner_->Condition(graph_);
}

void GreedySolver::Init() {
  std::for_each(graph_->request_groups().begin(),
                graph_->request_groups().end(),
                std::bind(&GreedySolver::GetCaps, this, std::placeholders::_1));

  std::for_each(graph_->supply_groups().begin(),
                graph_->supply_groups().end(),
                std::bind(&GreedySolver::GetCaps, this, std::placeholders::_1));
}

double GreedySolver::SolveGraph() {
  double pseudo_cost = PseudoCost();  // from ExchangeSolver API
  Condition();
  obj_ = 0;
  unmatched_ = 0;
  n_qty_.clear();

  Init();
  
  std::for_each(graph_->request_groups().begin(),
                graph_->request_groups().end(),
                std::bind(&GreedySolver::GreedilySatisfySet, this,
                          std::placeholders::_1));

  obj_ += unmatched_ * pseudo_cost;
  return obj_;
}

double GreedySolver::Capacity(const Arc& a, double u_curr_qty,
                              double v_curr_qty) {
  bool min = true;
  double ucap = Capacity(a.unode(), a, !min, u_curr_qty);
  double vcap = Capacity(a.vnode(), a, min, v_curr_qty);

  CLOG(cyclus::LEV_DEBUG1) << "Capacity for unode of arc: " << ucap;
  CLOG(cyclus::LEV_DEBUG1) << "Capacity for vnode of arc: " << vcap;
  CLOG(cyclus::LEV_DEBUG1) << "Capacity for arc         : "
                           << std::min(ucap, vcap);

  return std::min(ucap, vcap);
}

double GreedySolver::Capacity(ExchangeNode::Ptr n, const Arc& a, bool min_cap,
                              double curr_qty) {
  if (n->group == NULL) {
    throw cyclus::StateError(
        "An notion of node capacity requires a nodegroup.");
  }

  if (n->unit_capacities[a].size() == 0) {
    return n->qty - curr_qty;
  }

  std::vector<double>& unit_caps = n->unit_capacities[a];
  const std::vector<double>& group_caps = grp_caps_[n->group];
  std::vector<double> caps;
  double grp_cap, u_cap, cap;

  for (int i = 0; i < unit_caps.size(); i++) {
    grp_cap = group_caps[i];
    u_cap = unit_caps[i];
    cap = grp_cap / u_cap;
    CLOG(cyclus::LEV_DEBUG1) << "Capacity for node: ";
    CLOG(cyclus::LEV_DEBUG1) << "   group capacity: " << grp_cap;
    CLOG(cyclus::LEV_DEBUG1) << "    unit capacity: " << u_cap;
    CLOG(cyclus::LEV_DEBUG1) << "         capacity: " << cap;

    // special case for unlimited capacities
    if (grp_cap == std::numeric_limits<double>::max()) {
      caps.push_back(std::numeric_limits<double>::max());
    } else {
      caps.push_back(cap);
    }
  }

  if (min_cap) {  // the smallest value is constraining (for bids)
    cap = *std::min_element(caps.begin(), caps.end());
  } else {  // the largest value must be met (for requests)
    cap = *std::max_element(caps.begin(), caps.end());
  }
  return std::min(cap, n->qty - curr_qty);
}

void GreedySolver::GetCaps(ExchangeNodeGroup::Ptr g) {
  for (int i = 0; i != g->nodes().size(); i++) {
    n_qty_[g->nodes()[i]] = 0;
  }
  grp_caps_[g.get()] = g->capacities();
}

void GreedySolver::GreedilySatisfySet(RequestGroup::Ptr prs) {
  std::vector<ExchangeNode::Ptr>& nodes = prs->nodes();
  std::stable_sort(nodes.begin(), nodes.end(), AvgPrefComp(graph_));

  std::vector<ExchangeNode::Ptr>::iterator req_it = nodes.begin();
  double target = prs->qty();
  double match = 0;

  ExchangeNode::Ptr u, v;
  std::vector<Arc>::const_iterator arc_it;
  std::vector<Arc> sorted;
  double remain, tomatch, excl_val;

  CLOG(LEV_DEBUG1) << "Greedy Solving for " << target
                   << " amount of a resource.";

  while ((match <= target) && (req_it != nodes.end())) {
    // this if statement is needed because map.at() will throw if the key does
    // not exist, which is a corner case for when there is a request with no bid
    // arcs associated with it
    if (graph_->node_arc_map().count(*req_it) > 0) {
      const std::vector<Arc>& arcs = graph_->node_arc_map().at(*req_it);
      sorted = std::vector<Arc>(arcs);  // make a copy for now
      // ReqPrefComp no longer needs shift since it uses arc.pref() directly
      std::stable_sort(sorted.begin(), sorted.end(), ReqPrefComp(0.0));
      arc_it = sorted.begin();

      while ((match <= target) && (arc_it != sorted.end())) {
        remain = target - match;
        const Arc& a = *arc_it;
        u = a.unode();
        v = a.vnode();
        // capacity adjustment
        tomatch = std::min(remain, Capacity(a, n_qty_[u], n_qty_[v]));

        // exclusivity adjustment
        if (arc_it->exclusive()) {
          excl_val = a.excl_val();

          // this careful float comparison is vital for preventing false
          // positive constraint violations w.r.t. exclusivity-related capacity.
          double dist = boost::math::float_distance(tomatch, excl_val);
          if (dist >= float_ulp_eq) {
            tomatch = 0;
          } else {
            tomatch = excl_val;
          }
        }

        if (tomatch > eps()) {
          CLOG(LEV_DEBUG1) << "Greedy Solver is matching " << tomatch
                           << " amount of a resource.";
          UpdateCapacity(u, a, tomatch);
          UpdateCapacity(v, a, tomatch);
          n_qty_[u] += tomatch;
          n_qty_[v] += tomatch;
          graph_->AddMatch(a, tomatch);

          match += tomatch;
          // Use stored arc weight from pref() to avoid recalculation and ensure consistency
          UpdateObj(tomatch, a.pref());
        }
        ++arc_it;
      }  // while( (match =< target) && (arc_it != arcs.end()) )
    }  // if(graph_->node_arc_map().count(*req_it) > 0)
    ++req_it;
  }  // while( (match =< target) && (req_it != nodes.end()) )

  unmatched_ += target - match;
}

void GreedySolver::UpdateObj(double qty, double arc_weight) {
  // updates minimizing objective (arc_weight is the cost and the objective is cost * flow)
  obj_ += qty * arc_weight;
}

void GreedySolver::UpdateCapacity(ExchangeNode::Ptr n, const Arc& a,
                                  double qty) {
  using cyclus::IsNegative;
  using cyclus::ValueError;

  std::vector<double>& unit_caps = n->unit_capacities[a];
  std::vector<double>& caps = grp_caps_[n->group];
  assert(unit_caps.size() == caps.size());
  for (int i = 0; i < caps.size(); i++) {
    double prev = caps[i];
    // special case for unlimited capacities
    CLOG(cyclus::LEV_DEBUG1) << "Updating capacity value from: " << prev;
    caps[i] = (prev == std::numeric_limits<double>::max())
                  ? std::numeric_limits<double>::max()
                  : prev - qty * unit_caps[i];
    CLOG(cyclus::LEV_DEBUG1) << "                          to: " << caps[i];
  }

  if (IsNegative(n->qty - qty)) {
    std::stringstream ss;
    ss << "A bid for " << n->commod << " was set at " << n->qty
       << " but has been matched to a higher value " << qty
       << ". This could be due to a problem with your "
       << "bid portfolio constraints.";
    throw ValueError(ss.str());
  }
}

}  // namespace cyclus

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
19	GreedySolver::GreedySolver(bool exclusive_orders, GreedyPreconditioner* c)	50✔
20	: conditioner_(c), ExchangeSolver(exclusive_orders), shift_(0.0) {}	50✔
21
22	GreedySolver::GreedySolver(bool exclusive_orders)	4,184✔
23	: ExchangeSolver(exclusive_orders), shift_(0.0) {	4,184✔
24	conditioner_ = new cyclus::GreedyPreconditioner();	4,184✔
25	}	4,184✔
26
27	GreedySolver::GreedySolver(GreedyPreconditioner* c)	×
NEW 28	: conditioner_(c), ExchangeSolver(true), shift_(0.0) {}	×
29
30	GreedySolver::GreedySolver() : ExchangeSolver(true), shift_(0.0) {	13✔
31	conditioner_ = new cyclus::GreedyPreconditioner();	13✔
32	}	13✔
33
34	GreedySolver::~GreedySolver() {	4,493✔
35	if (conditioner_ != NULL) delete conditioner_;	4,247✔
36	}	4,493✔
37
38	void GreedySolver::Condition() {	5,030✔
39	if (conditioner_ != NULL) conditioner_->Condition(graph_);	5,030✔
40	}	5,030✔
41
42	void GreedySolver::Init() {	5,030✔
43	std::for_each(graph_->request_groups().begin(),
44	graph_->request_groups().end(),	5,030✔
45	std::bind(&GreedySolver::GetCaps, this, std::placeholders::_1));	5,030✔
46
47	std::for_each(graph_->supply_groups().begin(),
48	graph_->supply_groups().end(),	5,030✔
49	std::bind(&GreedySolver::GetCaps, this, std::placeholders::_1));	5,030✔
50	}	5,030✔
51
52	double GreedySolver::SolveGraph() {	5,029✔
53	double pseudo_cost = PseudoCost(); // from ExchangeSolver API	5,029✔
54	Condition();	5,029✔
55	obj_ = 0;	5,029✔
56	unmatched_ = 0;	5,029✔
57	n_qty_.clear();
58
59	Init();	5,029✔
60
61	std::for_each(graph_->request_groups().begin(),
62	graph_->request_groups().end(),	5,029✔
63	std::bind(&GreedySolver::GreedilySatisfySet, this,	5,029✔
64	std::placeholders::_1));
65
66	obj_ += unmatched_ * pseudo_cost;	5,029✔
67	return obj_;	5,029✔
68	}
69
70	double GreedySolver::Capacity(const Arc& a, double u_curr_qty,	548,635✔
71	double v_curr_qty) {
72	bool min = true;
73	double ucap = Capacity(a.unode(), a, !min, u_curr_qty);	1,097,270✔
74	double vcap = Capacity(a.vnode(), a, min, v_curr_qty);	548,635✔
75
76	CLOG(cyclus::LEV_DEBUG1) << "Capacity for unode of arc: " << ucap;	548,635✔
77	CLOG(cyclus::LEV_DEBUG1) << "Capacity for vnode of arc: " << vcap;	548,635✔
78	CLOG(cyclus::LEV_DEBUG1) << "Capacity for arc : "	548,635✔
79	<< std::min(ucap, vcap);	×
80
81	return std::min(ucap, vcap);	548,635✔
82	}
83
84	double GreedySolver::Capacity(ExchangeNode::Ptr n, const Arc& a, bool min_cap,	1,097,270✔
85	double curr_qty) {
86	if (n->group == NULL) {	1,097,270✔
87	throw cyclus::StateError(	×
88	"An notion of node capacity requires a nodegroup.");	×
89	}
90
91	if (n->unit_capacities[a].size() == 0) {	1,097,270✔
92	return n->qty - curr_qty;	12✔
93	}
94
95	std::vector<double>& unit_caps = n->unit_capacities[a];	1,097,258✔
96	const std::vector<double>& group_caps = grp_caps_[n->group];	1,097,258✔
97	std::vector<double> caps;
98	double grp_cap, u_cap, cap;
99
100	for (int i = 0; i < unit_caps.size(); i++) {	2,201,304✔
101	grp_cap = group_caps[i];	1,104,046✔
102	u_cap = unit_caps[i];	1,104,046✔
103	cap = grp_cap / u_cap;	1,104,046✔
104	CLOG(cyclus::LEV_DEBUG1) << "Capacity for node: ";	1,104,046✔
105	CLOG(cyclus::LEV_DEBUG1) << " group capacity: " << grp_cap;	1,104,046✔
106	CLOG(cyclus::LEV_DEBUG1) << " unit capacity: " << u_cap;	1,104,046✔
107	CLOG(cyclus::LEV_DEBUG1) << " capacity: " << cap;	1,104,046✔
108
109	// special case for unlimited capacities
110	if (grp_cap == std::numeric_limits<double>::max()) {	1,104,046✔
111	caps.push_back(std::numeric_limits<double>::max());	×
112	} else {
113	caps.push_back(cap);	1,104,046✔
114	}
115	}
116
117	if (min_cap) { // the smallest value is constraining (for bids)	1,097,258✔
118	cap = *std::min_element(caps.begin(), caps.end());	548,623✔
119	} else { // the largest value must be met (for requests)
120	cap = *std::max_element(caps.begin(), caps.end());	548,635✔
121	}
122	return std::min(cap, n->qty - curr_qty);	1,097,421✔
123	}
124
125	void GreedySolver::GetCaps(ExchangeNodeGroup::Ptr g) {	26,344✔
126	for (int i = 0; i != g->nodes().size(); i++) {	584,865✔
127	n_qty_[g->nodes()[i]] = 0;	558,521✔
128	}
129	grp_caps_[g.get()] = g->capacities();	26,344✔
130	}	26,344✔
131
132	void GreedySolver::GreedilySatisfySet(RequestGroup::Ptr prs) {	9,865✔
133	std::vector<ExchangeNode::Ptr>& nodes = prs->nodes();
134	std::stable_sort(nodes.begin(), nodes.end(), AvgPrefComp(graph_));	9,865✔
135
136	std::vector<ExchangeNode::Ptr>::iterator req_it = nodes.begin();
137	double target = prs->qty();
138	double match = 0;
139
140	ExchangeNode::Ptr u, v;	9,865✔
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	9,865✔
146	<< " amount of a resource.";	×
147
148	while ((match <= target) && (req_it != nodes.end())) {	19,748✔
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) {	9,883✔
153	const std::vector<Arc>& arcs = graph_->node_arc_map().at(*req_it);	9,881✔
154	sorted = std::vector<Arc>(arcs); // make a copy for now	9,881✔
155	// ReqPrefComp no longer needs shift since it uses arc.pref() directly
156	std::stable_sort(sorted.begin(), sorted.end(), ReqPrefComp(0.0));	9,881✔
157	arc_it = sorted.begin();
158
159	while ((match <= target) && (arc_it != sorted.end())) {	558,514✔
160	remain = target - match;	548,633✔
161	const Arc& a = *arc_it;
162	u = a.unode();	548,633✔
163	v = a.vnode();	548,633✔
164	// capacity adjustment
165	tomatch = std::min(remain, Capacity(a, n_qty_[u], n_qty_[v]));	753,511✔
166
167	// exclusivity adjustment
168	if (arc_it->exclusive()) {	548,633✔
169	excl_val = a.excl_val();
170
171	// this careful float comparison is vital for preventing false
172	// positive constraint violations w.r.t. exclusivity-related capacity.
173	double dist = boost::math::float_distance(tomatch, excl_val);
174	if (dist >= float_ulp_eq) {	58✔
175	tomatch = 0;
176	} else {
177	tomatch = excl_val;
178	}
179	}
180
181	if (tomatch > eps()) {	548,633✔
182	CLOG(LEV_DEBUG1) << "Greedy Solver is matching " << tomatch	12,105✔
183	<< " amount of a resource.";	×
184	UpdateCapacity(u, a, tomatch);	24,210✔
185	UpdateCapacity(v, a, tomatch);	12,105✔
186	n_qty_[u] += tomatch;	12,105✔
187	n_qty_[v] += tomatch;	12,105✔
188	graph_->AddMatch(a, tomatch);	12,105✔
189
190	match += tomatch;	12,105✔
191	// Use stored arc weight from pref() to avoid recalculation and ensure consistency
192	UpdateObj(tomatch, a.pref());	12,105✔
193	}
194	++arc_it;
195	} // while( (match =< target) && (arc_it != arcs.end()) )
196	} // if(graph_->node_arc_map().count(*req_it) > 0)
197	++req_it;
198	} // while( (match =< target) && (req_it != nodes.end()) )
199
200	unmatched_ += target - match;	9,865✔
201	}	19,730✔
202
203	void GreedySolver::UpdateObj(double qty, double arc_weight) {	12,105✔
204	// updates minimizing objective (arc_weight is the cost and the objective is cost * flow)
205	obj_ += qty * arc_weight;	12,105✔
206	}	12,105✔
207
208	void GreedySolver::UpdateCapacity(ExchangeNode::Ptr n, const Arc& a,	24,210✔
209	double qty) {
210	using cyclus::IsNegative;
211	using cyclus::ValueError;
212
213	std::vector<double>& unit_caps = n->unit_capacities[a];	24,210✔
214	std::vector<double>& caps = grp_caps_[n->group];	24,210✔
215	assert(unit_caps.size() == caps.size());
216	for (int i = 0; i < caps.size(); i++) {	55,196✔
217	double prev = caps[i];	30,986✔
218	// special case for unlimited capacities
219	CLOG(cyclus::LEV_DEBUG1) << "Updating capacity value from: " << prev;	30,986✔
220	caps[i] = (prev == std::numeric_limits<double>::max())	30,986✔
221	? std::numeric_limits<double>::max()	30,986✔
222	: prev - qty * unit_caps[i];	30,986✔
223	CLOG(cyclus::LEV_DEBUG1) << " to: " << caps[i];	30,986✔
224	}
225
226	if (IsNegative(n->qty - qty)) {	24,210✔
227	std::stringstream ss;	×
228	ss << "A bid for " << n->commod << " was set at " << n->qty	×
229	<< " but has been matched to a higher value " << qty
230	<< ". This could be due to a problem with your "
231	<< "bid portfolio constraints.";	×
232	throw ValueError(ss.str());	×
233	}	×
234	}	24,210✔
235
236	} // namespace cyclus

cyclus / cyclus / 22921176376

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous