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