Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : #pragma once
      10             : 
      11             : #include "MooseError.h"
      12             : #include <numeric>
      13             : #include <vector>
      14             : #include <deque>
      15             : #include "libmesh/libmesh_common.h"
      16             : 
      17             : using namespace libMesh;
      18             : 
      19             : namespace StochasticTools
      20             : {
      21             : /* Compute the Cartesian Product of the supplied vectors.
      22             :  * http://phrogz.net/lazy-cartesian-product
      23             :  * https://github.com/iamtheburd/lazy-cartesian-product-python/blob/master/LazyCartesianProduct.py
      24             :  */
      25             : template <class T>
      26             : class CartesianProduct
      27             : {
      28             : public:
      29             :   CartesianProduct(const std::vector<std::vector<T>> & items);
      30             : 
      31             :   /// Compute the complete Cartesian product matrix
      32             :   std::vector<std::vector<T>> computeMatrix() const;
      33             : 
      34             :   /// Compute specified row of Cartesian product matrix
      35             :   std::vector<T> computeRow(std::size_t row) const;
      36             : 
      37             :   /// Compute specific value, given row and column, of the Cartesian product matrix
      38             :   T computeValue(std::size_t row, std::size_t col) const;
      39             : 
      40             :   /// Total number of rows in the complete matrix
      41      728720 :   std::size_t numRows() const { return _n_rows; }
      42             : 
      43             :   /// Total number of columns in the complete matrix
      44        1102 :   std::size_t numCols() const { return _n_cols; }
      45             : 
      46             : protected:
      47             :   /// Number of rows/columns
      48             :   const std::size_t _n_rows;
      49             :   const std::size_t _n_cols;
      50             : 
      51             : private:
      52             :   /// Data used to create Cartesian product
      53             :   /// use a copy because a temporary can be supplied, as is the case in the CartesianProductSampler
      54             :   const std::vector<std::vector<T>> _items;
      55             : 
      56             :   /// Containers for lazy Cartesian product calculation
      57             :   std::deque<unsigned int> _denomenators;
      58             :   std::deque<unsigned int> _moduli;
      59             : 
      60             :   /// Helper to compute the rows in initialization list to allow _n_rows to be const
      61             :   static std::size_t computeRowCount(const std::vector<std::vector<T>> & items);
      62             : };
      63             : 
      64             : template <typename T>
      65        4315 : CartesianProduct<T>::CartesianProduct(const std::vector<std::vector<T>> & items)
      66        4315 :   : _n_rows(computeRowCount(items)), _n_cols(items.size()), _items(items)
      67             : {
      68             :   dof_id_type d = 1;
      69             :   for (typename std::vector<std::vector<T>>::const_reverse_iterator iter = _items.rbegin();
      70       24373 :        iter != _items.rend();
      71             :        ++iter)
      72             :   {
      73             :     std::size_t n = iter->size();
      74       20058 :     _denomenators.push_front(d);
      75       20058 :     _moduli.push_front(n);
      76       20058 :     d *= n;
      77             :   }
      78        4315 : }
      79             : 
      80             : template <typename T>
      81             : std::vector<std::vector<T>>
      82           1 : CartesianProduct<T>::computeMatrix() const
      83             : {
      84           1 :   std::vector<std::vector<T>> output(_n_rows, std::vector<T>(_n_cols));
      85          25 :   for (std::size_t row = 0; row < _n_rows; ++row)
      86          96 :     for (std::size_t col = 0; col < _n_cols; ++col)
      87          72 :       output[row][col] = computeValue(row, col);
      88           1 :   return output;
      89             : }
      90             : 
      91             : template <typename T>
      92             : std::vector<T>
      93     1229888 : CartesianProduct<T>::computeRow(std::size_t row) const
      94             : {
      95     1229888 :   std::vector<T> output(_n_cols);
      96     6608886 :   for (std::size_t col = 0; col < _n_cols; ++col)
      97     5378998 :     output[col] = computeValue(row, col);
      98     1229888 :   return output;
      99             : }
     100             : 
     101             : template <typename T>
     102             : T
     103     7429885 : CartesianProduct<T>::computeValue(std::size_t row, std::size_t col) const
     104             : {
     105             :   mooseAssert(row < _n_rows, "Row index out of range.");
     106             :   mooseAssert(col < _n_cols, "Column index out of range.");
     107     7429885 :   return _items[col][(row / _denomenators[col]) % _moduli[col]];
     108             : }
     109             : 
     110             : template <typename T>
     111             : std::size_t
     112             : CartesianProduct<T>::computeRowCount(const std::vector<std::vector<T>> & items)
     113             : {
     114             :   std::size_t n_rows = 1;
     115       24373 :   for (const auto & inner : items)
     116       20058 :     n_rows *= inner.size();
     117             :   return n_rows;
     118             : }
     119             : 
     120             : /*
     121             :  * Add ability to compute weighting values with the Cartesian Product
     122             :  */
     123             : template <class T, class W>
     124        1595 : class WeightedCartesianProduct : public CartesianProduct<T>
     125             : {
     126             : public:
     127             :   WeightedCartesianProduct(const std::vector<std::vector<T>> & items,
     128             :                            const std::vector<std::vector<W>> & weights);
     129             : 
     130             :   /// Compute complete vector of weights
     131             :   std::vector<W> computeWeightVector() const;
     132             : 
     133             :   /// Compute specific weight value, given row
     134             :   W computeWeight(std::size_t row) const;
     135             : 
     136             : private:
     137             :   /// Data used to create Cartesian product; use a copy because a temporary can be supplied
     138             :   const CartesianProduct<W> _weight;
     139             : };
     140             : 
     141             : template <typename T, typename W>
     142        1539 : WeightedCartesianProduct<T, W>::WeightedCartesianProduct(
     143             :     const std::vector<std::vector<T>> & items, const std::vector<std::vector<W>> & weights)
     144        1539 :   : CartesianProduct<T>(items), _weight(weights)
     145             : {
     146             :   mooseAssert(items.size() == weights.size(),
     147             :               "The supplied items and weights must be the same size.");
     148             :   for (std::size_t i = 0; i < items.size(); ++i)
     149             :     mooseAssert(items[i].size() == weights[i].size(),
     150             :                 "Internal vector of the supplied items and weights must be the same size.");
     151        1539 : }
     152             : 
     153             : template <typename T, typename W>
     154             : std::vector<W>
     155           1 : WeightedCartesianProduct<T, W>::computeWeightVector() const
     156             : {
     157           1 :   std::vector<W> output(this->_n_rows);
     158          25 :   for (std::size_t i = 0; i < output.size(); ++i)
     159          24 :     output[i] = computeWeight(i);
     160           1 :   return output;
     161           0 : }
     162             : 
     163             : template <typename T, typename W>
     164             : W
     165      503804 : WeightedCartesianProduct<T, W>::computeWeight(std::size_t row) const
     166             : {
     167      503804 :   std::vector<W> vec = _weight.computeRow(row);
     168      503804 :   return std::accumulate(vec.begin(), vec.end(), static_cast<W>(1), std::multiplies<W>());
     169      503804 : }
     170             : } // namespace