add normal and uniform packaging strategies

CHANGELOG.rst

@@ -13,6 +13,7 @@ Since last release
 * Added macros in ``cmake/CyclusBuildSetup.cmake`` for common CMake boilerplate (#1793)
 * Added ``doxygen-awesome-css`` to doxygen docs for style (#1787)
 * Added installation of files for building docs to share/cyclus/doc (#1807)
+* New packaging strategies uniform and normal (#1813)
 
 **Changed:**
 
@@ -26,7 +27,7 @@ Since last release
 * Warning and limits on number of packages that can be created from a resource at once (#1771)
 * Use keep_packaging instead of unpackaged in ResBuf (#1778)
 * Temporarily pin Boost libraries to <1.86.0 (#1796)
-* Package GetFillMass returns fill mass and number of packages filled at that mass (#1790)
+* Package GetFillMass returns vector of masses (#1790, #1813)
 
 **Removed:**
 

share/cyclus-flat.rng.in

@@ -151,9 +151,9 @@ datatypeLibrary="http://www.w3.org/2001/XMLSchema-datatypes">
     <element name="package">
       <interleave>
         <element name="name"><text/></element>
-        <element name="fill_min"><data type="double"/></element>
-        <element name="fill_max"><data type="double"/></element>
-        <element name="strategy"><text/></element>
+        <optional><element name="fill_min"><data type="double"/></element></optional>
+        <optional><element name="fill_max"><data type="double"/></element></optional>
+        <optional><element name="strategy"><text/></element></optional>
       </interleave>
     </element>
   </zeroOrMore>
@@ -162,9 +162,9 @@ datatypeLibrary="http://www.w3.org/2001/XMLSchema-datatypes">
     <element name="transportunit">
       <interleave>
         <element name="name"><text/></element>
-        <element name="fill_min"><data type="nonNegativeInteger"/></element>
-        <element name="fill_max"><data type="nonNegativeInteger"/></element>
-        <element name="strategy"><text/></element>
+        <optional><element name="fill_min"><data type="nonNegativeInteger"/></element></optional>
+        <optional><element name="fill_max"><data type="nonNegativeInteger"/></element></optional>
+        <optional><element name="strategy"><text/></element></optional>
       </interleave>
     </element>
   </zeroOrMore>

share/cyclus.rng.in

@@ -183,9 +183,9 @@ datatypeLibrary="http://www.w3.org/2001/XMLSchema-datatypes">
     <element name="package">
       <interleave>
         <element name="name"><text/></element>
-        <element name="fill_min"><data type="double"/></element>
-        <element name="fill_max"><data type="double"/></element>
-        <element name="strategy"><text/></element>
+        <optional><element name="fill_min"><data type="double"/></element></optional>
+        <optional><element name="fill_max"><data type="double"/></element></optional>
+        <optional><element name="strategy"><text/></element></optional>
       </interleave>
     </element>
   </zeroOrMore>
@@ -194,9 +194,9 @@ datatypeLibrary="http://www.w3.org/2001/XMLSchema-datatypes">
     <element name="transportunit">
       <interleave>
         <element name="name"><text/></element>
-        <element name="fill_min"><data type="nonNegativeInteger"/></element>
-        <element name="fill_max"><data type="nonNegativeInteger"/></element>
-        <element name="strategy"><text/></element>
+        <optional><element name="fill_min"><data type="nonNegativeInteger"/></element></optional>
+        <optional><element name="fill_max"><data type="nonNegativeInteger"/></element></optional>
+        <optional><element name="strategy"><text/></element></optional>
       </interleave>
     </element>
   </zeroOrMore>

src/package.cc

@@ -30,12 +30,28 @@ Package::Ptr& Package::unpackaged() {
   return unpackaged_;
 }
 
-std::pair<double, int> Package::GetFillMass(double qty) {
+void Package::SetDistribution() {
+  if (strategy_ == "uniform") {
+    dist_ = UniformDoubleDist::Ptr (new UniformDoubleDist(fill_min_, fill_max_));
+  } else if (strategy_ == "normal") {
+    dist_ = NormalDoubleDist::Ptr (new NormalDoubleDist(
+      (fill_min_ + fill_max_) / 2,
+      (fill_max_ - fill_min_) / 6,
+      fill_min_,
+      fill_max_));
+  }
+}
+
+std::vector<double> Package::GetFillMass(double qty) {
+  std::vector<double> packages;
   if ((qty - fill_min_) < -eps_rsrc()) {
     // less than one pkg of material available
-    return std::pair<double, int> (0, 0);
+    return packages;
   }
 
+  // simple check for whether vector limits *might* be exceeded
+  ExceedsSplitLimits(qty / fill_max_);
+
   double fill_mass;
   int num_at_fill_mass;
 
@@ -52,9 +68,34 @@ std::pair<double, int> Package::GetFillMass(double qty) {
       fill_mass = fill_max_;
     }
   }
-  fill_mass = std::min(qty, fill_mass);
-  num_at_fill_mass = static_cast<int>(std::floor(qty / fill_mass));
-  return std::pair<double, int>(fill_mass, num_at_fill_mass);
+
+  if (strategy_ == "first" || strategy_ == "equal") {
+    fill_mass = std::min(qty, fill_mass);
+    num_at_fill_mass = static_cast<int>(std::floor(qty / fill_mass));
+    ExceedsSplitLimits(num_at_fill_mass);
+    packages.assign(num_at_fill_mass, fill_mass);
+
+    qty -= num_at_fill_mass * fill_mass;
+  }
+
+  if (strategy_ == "uniform" || strategy_ == "normal") {
+    // only use random if a full package amount is available. if less than one
+    // full amount is available, below will fill a partial package (no random).
+    while ((qty > 0) && (qty >= fill_max_)) {
+      fill_mass = dist_->sample();
+      packages.push_back(fill_mass);
+      qty -= fill_mass;
+    }
+  }
+
+  if (qty > eps_rsrc() && qty >= fill_min_) {
+    // leftover material is enough to fill one more partial package. 
+    packages.push_back(qty);
+  }
+
+  Package::ExceedsSplitLimits(packages.size());
+
+  return packages;
 }
 
 Package::Package(std::string name, double fill_min, double fill_max,
@@ -65,9 +106,10 @@ Package::Package(std::string name, double fill_min, double fill_max,
         throw ValueError("can't create a new package with name 'unpackaged'");
       }
   }
-    if (strategy != "first" && strategy != "equal") {
+    if (!IsValidStrategy(strategy_)) {
       throw ValueError("Invalid strategy for package: " + strategy_);
     }
+    SetDistribution();
 }
 
 TransportUnit::Ptr TransportUnit::unrestricted_ = NULL;

src/package.h

@@ -9,6 +9,7 @@
 
 #include "error.h"
 #include "logger.h"
+#include "random_number_generator.h"
 
 namespace cyclus {
 
@@ -38,13 +39,21 @@ class Package {
     /// This tries to find the optimal number and fill mass of packages given
     /// the packaging limitations. It does this by calculating bounding fills, 
     /// floor(quantity/fill_min) and ceiling(quantity/fill_max). 
+    /// Packaging strategy "uniform" fills packages with a random mass between
+    /// fill_min and fill_max, if at least fill_max is available. If less than
+    /// fill_max is available, a partial package is filled with the total mass.
+    /// Packaging strategy "normal" fills packages with a random mass between
+    /// fill_min and fill_max, with a normal distribution. Mean is the middle
+    /// of fill_min and fill_max, standard deviation is 1/6 of the range such
+    /// that 3 sigma is the range. If less than fill_max is available, a 
+    /// partial package is filled with the total mass (no dist sampling).
     /// There might be a scenario where there is no solution, i.e. an integer
     /// number of packages cannot be filled with no remainder. In this case,
     /// the most effective fill strategy is to fill to the max. Numeric example:
     /// quantity = 5, fill_min = 3, fill_max = 4. num_min_fill = floor(5/3) = 1,
     /// num_max_fill = ceil(5/4) = 2. num_min_fill < num_max_fill, so fill to
     /// the max.
-    std::pair<double, int> GetFillMass(double qty);
+    std::vector<double> GetFillMass(double qty);
 
     // returns package name
     std::string name() const { return name_; }
@@ -86,6 +95,15 @@ class Package {
       return (std::numeric_limits<int>::max() / 10); 
     }
 
+    bool IsValidStrategy(std::string strategy) {
+      return strategy == "first"
+          || strategy == "equal"
+          || strategy == "uniform"
+          || strategy == "normal";
+    }
+
+    void SetDistribution();
+
   private:
     Package(std::string name, 
             double fill_min = 0, 
@@ -99,6 +117,7 @@ class Package {
     double fill_min_;
     double fill_max_;
     std::string strategy_;
+    DoubleDistribution::Ptr dist_ = NULL;
 };
 
 /// TransportUnit is a class that can be used in conjunction with packages to

src/resource.h

@@ -141,22 +141,17 @@ std::vector<typename T::Ptr> Resource::Package(Package::Ptr pkg) {
   std::vector<typename T::Ptr> ts_pkgd;
   typename T::Ptr t_pkgd;
 
-  std::pair<double, int> fill = pkg->GetFillMass(quantity());
-  double fill_mass = fill.first;
-  if (fill_mass == 0) {
+  std::vector<double> packages = pkg->GetFillMass(quantity());
+  if (packages.size() == 0) {
     return ts_pkgd;
   }
 
-  // Check if the number of packages is within the limits, including if
-  // int overflow is reached
-  int approx_num_pkgs = fill.second;
-  Package::ExceedsSplitLimits(approx_num_pkgs);
-
-  while (quantity() > 0 && (quantity() - pkg->fill_min()) >= -eps_rsrc()) {
-    double pkg_fill = std::min(quantity(), fill_mass);
+  for (int i = 0; i < packages.size(); ++i) {
+    double pkg_fill = packages[i];
     t_pkgd = boost::dynamic_pointer_cast<T>(PackageExtract(pkg_fill, pkg->name()));
-    ts_pkgd.push_back(t_pkgd); 
+    ts_pkgd.push_back(t_pkgd);
   }
+
   return ts_pkgd;
 }
 

src/toolkit/matl_sell_policy.cc

@@ -47,14 +47,22 @@ void MatlSellPolicy::set_package(std::string x) {
   // if no real context, only unpackaged can be used (keep default)
   if (manager() != NULL) {
     Package::Ptr pkg = manager()->context()->GetPackage(x);
-    std::pair<double, int> fill = pkg->GetFillMass(quantize_);
-    double pkg_fill = fill.first;
-    if ((pkg->name() != Package::unpackaged_name()) && (quantize_ > 0) &&
-        (std::fmod(quantize_, pkg_fill) > 0)) { 
-      std::stringstream ss;
-      ss << "Quantize " << quantize_ << " is not fully packagable based on fill min/max values (" 
-         << pkg->fill_min() << ", " << pkg->fill_max() << ")";
-      throw ValueError(ss.str());
+    std::stringstream ss;
+
+    if (quantize_ > 0 && pkg->name() != Package::unpackaged_name()) {
+      if (pkg->strategy() == "first" || pkg->strategy() == "equal") {
+        std::vector<double> fill = pkg->GetFillMass(quantize_);
+        if (fill.size() > 1 ||  std::fmod(fill.front(), quantize_) > 0) {
+          ss << "Quantize " << quantize_ 
+             << " is not fully packagable based on fill min/max values (" 
+             << pkg->fill_min() << ", " << pkg->fill_max() << ")";
+          throw ValueError(ss.str());
+        }
+      } else {
+        ss << "Package strategy " << pkg->strategy() 
+            << " is not allowed for sell policies with quantize.";
+        throw ValueError(ss.str());
+      }
     }
     package_ = pkg;
   }
@@ -64,16 +72,19 @@ void MatlSellPolicy::set_transport_unit(std::string x) {
   if (manager() != NULL) {
     TransportUnit::Ptr tu = manager()->context()->GetTransportUnit(x);
 
-    std::pair<double, int> fill = package_->GetFillMass(quantize_);
-    int num_pkgs = fill.second;
-    int max_shippable = tu->MaxShippablePackages(num_pkgs);
-
-    if ((tu->name() != TransportUnit::unrestricted_name()) && quantize_ > 0 &&
-        (max_shippable != num_pkgs))  {
-      std::stringstream ss;
-      ss << "Quantize " << quantize_ << " packages cannot be shipped according to transport unit fill min/max values (" << tu->fill_min() << ", "
-       << tu->fill_max() << ")";
-      throw ValueError(ss.str());
+    if ((quantize_ > 0) && (tu->name() != TransportUnit::unrestricted_name())) {
+      std::vector<double> fill = package_->GetFillMass(quantize_);
+      int num_pkgs = fill.size();
+      int max_shippable = tu->MaxShippablePackages(num_pkgs);
+
+      if (max_shippable != num_pkgs)  {
+        std::stringstream ss;
+        ss << "Quantize " << quantize_ 
+           << " packages cannot be shipped according to transport unit fill "
+           << "min/max values (" 
+           << tu->fill_min() << ", " << tu->fill_max() << ")";
+        throw ValueError(ss.str());
+      }
     }
     transport_unit_ = tu;
   }
@@ -199,34 +210,15 @@ std::set<BidPortfolio<Material>::Ptr> MatlSellPolicy::GetMatlBids(
     for (rit = requests.begin(); rit != requests.end(); ++rit) {
       req = *rit;
       qty = std::min(req->target()->quantity(), limit);
-      std::pair<double, int> fill = package_->GetFillMass(qty);
-      bid_qty = excl ? quantize_ : fill.first;
-      if (bid_qty != 0) {
-        n_full_bids = excl ? std::floor(qty / quantize_) : fill.second;
-
-        // Throw if number of bids above limit or if casting to int caused
-        // overflow to negative int limit 
-        std::string s;
-        if (manager() != NULL)
-          s + " Agent: "
-            + Trader::manager()->prototype() + "-"
-            + std::to_string(Trader::manager()->id()) + ". ";
-        s + "This is likely due to too much material (did you forget a "
-          + "throughput?) or small package limits relative to the "
-          + "quantity available. qty: " + std::to_string(qty) 
-          + ", and each bid would be: " + std::to_string(bid_qty);
-        Package::ExceedsSplitLimits(n_full_bids, s);
-
-        bids.assign(n_full_bids, bid_qty);
-
-        remaining_qty = qty - (n_full_bids * bid_qty);
-        if ((!excl) && (remaining_qty > 0) && 
-            (remaining_qty >= package_->fill_min())) {
-          // leftover material is enough to fill one more partial package. Add
-          // to bids
-          bids.push_back(remaining_qty);
-        }
+
+      std::vector<double> bids;
+      if (excl) {
+        int n_full_bids = static_cast<int>(std::floor(qty / quantize_));
+        bids.assign(n_full_bids, quantize_);
+      } else {
+        bids = package_->GetFillMass(qty);
       }
+
       // check transportability
       int shippable_pkgs = transport_unit_->MaxShippablePackages(bids.size());
       if (shippable_pkgs < bids.size()) {