From 30b8681731f6d49bb2ea12af2cd389e7b8bba3e4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Wed, 20 Feb 2019 11:56:56 +0100
Subject: [PATCH] Dynare++: make multinomial coeffs computation a method of
 IntSequence

Also improve on the comments.
---
 dynare++/tl/cc/int_sequence.cc  | 34 +++++++++++++++++++++++++++++++++
 dynare++/tl/cc/int_sequence.hh  |  4 ++++
 dynare++/tl/cc/sparse_tensor.cc |  3 +--
 dynare++/tl/cc/tensor.cc        | 20 -------------------
 dynare++/tl/cc/tensor.hh        |  9 ---------
 5 files changed, 39 insertions(+), 31 deletions(-)

diff --git a/dynare++/tl/cc/int_sequence.cc b/dynare++/tl/cc/int_sequence.cc
index f50d0f80e..221b45885 100644
--- a/dynare++/tl/cc/int_sequence.cc
+++ b/dynare++/tl/cc/int_sequence.cc
@@ -3,6 +3,7 @@
 #include "int_sequence.hh"
 #include "symmetry.hh"
 #include "tl_exception.hh"
+#include "pascal_triangle.hh"
 
 #include <iostream>
 #include <limits>
@@ -262,3 +263,36 @@ IntSequence::print() const
     std::cout << operator[](i) << ' ';
   std::cout << ']' << std::endl;
 }
+
+/* Here we calculate the multinomial coefficients
+    $\left(\matrix{a\cr b_1,\ldots,b_n}\right)$, where $a=b_1+\ldots+b_n$.
+
+   See:
+    https://en.wikipedia.org/wiki/Binomial_coefficient#Generalization_to_multinomials
+    https://en.wikipedia.org/wiki/Multinomial_theorem
+
+   For n=1, the coefficient is equal to 1.
+   For n=2, the multinomial coeffs correspond to the binomial coeffs, i.e. the binomial
+    (a; b) is equal to the multinomial (a; b,a-b).
+   For n>=3, we have the identity
+   $$\left(\matrix{a\cr b_1,\ldots,b_n}\right)=\left(\matrix{b_1+b_2\cr b_1}\right)\cdot
+   \left(\matrix{a\cr b_1+b_2,b_3,\ldots,b_n}\right)$$ (where the first factor
+   on the right hand side is to be interpreted as a binomial coefficient)
+
+   This number is exactly a number of unfolded indices corresponding to one
+   folded index, where the sequence $b_1,\ldots,b_n$ is the symmetry of the
+   index. This can be easily seen if the multinomial coefficient is interpreted
+   as the number of unique permutations of a word, where $a$ is the length of
+   the word, $n$ is the number of distinct letters, and the $b_i$ are the
+   number of repetitions of each letter. For example, for a symmetry of the
+   form $y^4 u^2 v^3$, we want to compute the number of permutations of the word
+   $yyyyuuvvv$. This is equal to the multinomial coefficient (9; 4,2,3). */
+
+int
+IntSequence::noverseq()
+{
+  if (size() == 0 || size() == 1)
+    return 1;
+  data[1] += data[0];
+  return PascalTriangle::noverk(data[1], data[0]) * IntSequence(*this, 1, size()).noverseq();
+}
diff --git a/dynare++/tl/cc/int_sequence.hh b/dynare++/tl/cc/int_sequence.hh
index 44bf0b3c9..2495f914c 100644
--- a/dynare++/tl/cc/int_sequence.hh
+++ b/dynare++/tl/cc/int_sequence.hh
@@ -158,6 +158,10 @@ public:
   bool isConstant() const;
   bool isSorted() const;
   void print() const;
+  // Compute multinomial coefficient (sum(this); this)
+  /* Warning: this operation is destructive; make a copy if you want to keep
+     the original sequence */
+  int noverseq();
 };
 
 #endif
diff --git a/dynare++/tl/cc/sparse_tensor.cc b/dynare++/tl/cc/sparse_tensor.cc
index 5c3d59b91..308c88c59 100644
--- a/dynare++/tl/cc/sparse_tensor.cc
+++ b/dynare++/tl/cc/sparse_tensor.cc
@@ -82,8 +82,7 @@ SparseTensor::getUnfoldIndexFillFactor() const
   while (start_col != m.end())
     {
       const IntSequence &key = (*start_col).first;
-      Symmetry s(key);
-      cnt += Tensor::noverseq(s);
+      cnt += Symmetry(key).noverseq();
       start_col = m.upper_bound(key);
     }
 
diff --git a/dynare++/tl/cc/tensor.cc b/dynare++/tl/cc/tensor.cc
index 30a7222f2..2c59caea7 100644
--- a/dynare++/tl/cc/tensor.cc
+++ b/dynare++/tl/cc/tensor.cc
@@ -5,26 +5,6 @@
 #include "tl_static.hh"
 #include "pascal_triangle.hh"
 
-// |Tensor::noverseq_ip| static method
-/* Here we calculate a generalized combination number
-   $\left(\matrix{a\cr b_1,\ldots,b_n}\right)$, where $a=b_1+\ldots+
-   b_n$. We use the identity
-   $$\left(\matrix{a\cr b_1,\ldots,b_n}\right)=\left(\matrix{b_1+b_2\cr b_1}\right)\cdot
-   \left(\matrix{a\cr b_1+b_2,b_3,\ldots,b_n}\right)$$
-
-   This number is exactly a number of unfolded indices corresponding to
-   one folded index, where the sequence $b_1,\ldots,b_n$ is the symmetry
-   of the index. */
-
-int
-Tensor::noverseq_ip(IntSequence &s)
-{
-  if (s.size() == 0 || s.size() == 1)
-    return 1;
-  s[1] += s[0];
-  return PascalTriangle::noverk(s[1], s[0]) * noverseq(IntSequence(s, 1, s.size()));
-}
-
 /* Here we increment a given sequence within full symmetry given by
    |nv|, which is number of variables in each dimension. The underlying
    tensor is unfolded, so we increase the rightmost by one, and if it is
diff --git a/dynare++/tl/cc/tensor.hh b/dynare++/tl/cc/tensor.hh
index 24f24476a..0c90a26b0 100644
--- a/dynare++/tl/cc/tensor.hh
+++ b/dynare++/tl/cc/tensor.hh
@@ -217,15 +217,6 @@ public:
   {
     return in_end;
   }
-
-  static int
-  noverseq(const IntSequence &s)
-  {
-    IntSequence seq(s);
-    return noverseq_ip((IntSequence &) s);
-  }
-private:
-  static int noverseq_ip(IntSequence &s);
 };
 
 /* Here is an abstraction for unfolded tensor. We provide a pure