skip Cholesky decomposition in is>>n_mv_dist
authorAlexandre Oliva <aoliva@gcc.gnu.org>
Fri, 9 Aug 2019 09:20:58 +0000 (09:20 +0000)
committerAlexandre Oliva <aoliva@gcc.gnu.org>
Fri, 9 Aug 2019 09:20:58 +0000 (09:20 +0000)
normal_mv_distribution maintains the variance-covariance matrix param
in Cholesky-decomposed form.  Existing param_type constructors, when
taking a full or lower-triangle varcov matrix, perform Cholesky
decomposition to convert it to the internal representation.  This
internal representation is visible both in the varcov() result, and in
the streamed-out representation of a normal_mv_distribution object.

The problem is that when that representation is streamed back in, the
read-back decomposed varcov matrix is used as a lower-triangle
non-decomposed varcov matrix, and it undergoes Cholesky decomposition
again.  So, each cycle of stream-out/stream-in changes the varcov
matrix to its "square root", instead of restoring the original
params.

This patch includes Corentin's changes that introduce verification in
testsuite/ext/random/normal_mv_distribution/operators/serialize.cc and
other similar tests that the object read back in compares equal to the
written-out object: the modified tests pass only if (u == v).

This patch also fixes the error exposed by his change, introducing an
alternate private constructor for param_type, used only by operator>>.

for  libstdc++-v3/ChangeLog

* include/ext/random
(normal_mv_distribution::param_type::param_type): New private
ctor taking a decomposed varcov matrix, for use by...
(operator>>): ... this, befriended.
* include/ext/random.tcc (operator>>): Use it.
(normal_mv_distribution::param_type::_M_init_lower): Adjust
member function name in exception message.

for  libstdc++-v3/ChangeLog
from  Corentin Gay  <gay@adacore.com>

* testsuite/ext/random/beta_distribution/operators/serialize.cc,
testsuite/ext/random/hypergeometric_distribution/operators/serialize.cc,
testsuite/ext/random/normal_mv_distribution/operators/serialize.cc,
testsuite/ext/random/triangular_distribution/operators/serialize.cc,
testsuite/ext/random/von_mises_distribution/operators/serialize.cc:
Add call to `VERIFY`.

From-SVN: r274233

libstdc++-v3/ChangeLog
libstdc++-v3/include/ext/random
libstdc++-v3/include/ext/random.tcc
libstdc++-v3/testsuite/ext/random/beta_distribution/operators/serialize.cc
libstdc++-v3/testsuite/ext/random/hypergeometric_distribution/operators/serialize.cc
libstdc++-v3/testsuite/ext/random/normal_mv_distribution/operators/serialize.cc
libstdc++-v3/testsuite/ext/random/triangular_distribution/operators/serialize.cc
libstdc++-v3/testsuite/ext/random/von_mises_distribution/operators/serialize.cc

index 29418eb9ad96d07c67b73918b7ac4508402faa66..5c02cb36488ad65b4a5d7b0acdee0f5abe98543c 100644 (file)
@@ -1,3 +1,22 @@
+2019-08-09  Corentin Gay  <gay@adacore.com>
+
+       * testsuite/ext/random/beta_distribution/operators/serialize.cc,
+       testsuite/ext/random/hypergeometric_distribution/operators/serialize.cc,
+       testsuite/ext/random/normal_mv_distribution/operators/serialize.cc,
+       testsuite/ext/random/triangular_distribution/operators/serialize.cc,
+       testsuite/ext/random/von_mises_distribution/operators/serialize.cc:
+       Add call to `VERIFY`.
+
+2019-08-09  Alexandre Oliva <oliva@adacore.com>
+
+       * include/ext/random
+       (normal_mv_distribution::param_type::param_type): New private
+       ctor taking a decomposed varcov matrix, for use by...
+       (operator>>): ... this, befriended.
+       * include/ext/random.tcc (operator>>): Use it.
+       (normal_mv_distribution::param_type::_M_init_lower): Adjust
+       member function name in exception message.
+
 2019-08-08  Jonathan Wakely  <jwakely@redhat.com>
 
        P0325R4 to_array from LFTS with updates
index 41a2962c8f6e5de9bd5a90e960ba2c8822c364ff..d5574e02ba02cc432da82de06f912b586c86f9c1 100644 (file)
@@ -752,6 +752,21 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
                                _InputIterator2 __varbegin,
                                _InputIterator2 __varend);
 
+       // param_type constructors apply Cholesky decomposition to the
+       // varcov matrix in _M_init_full and _M_init_lower, but the
+       // varcov matrix output ot a stream is already decomposed, so
+       // we need means to restore it as-is when reading it back in.
+       template<size_t _Dimen1, typename _RealType1,
+                typename _CharT, typename _Traits>
+       friend std::basic_istream<_CharT, _Traits>&
+       operator>>(std::basic_istream<_CharT, _Traits>& __is,
+                  __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>&
+                  __x);
+       param_type(std::array<_RealType, _Dimen> const &__mean,
+                  std::array<_RealType, _M_t_size> const &__varcov)
+         : _M_mean (__mean), _M_t (__varcov)
+       {}
+
        std::array<_RealType, _Dimen> _M_mean;
        std::array<_RealType, _M_t_size> _M_t;
       };
index 31dc33a2555ed423545eb738ece65108cde9c57b..a8a49a3a9fa6a48a2911da122565aeded03d3f03 100644 (file)
@@ -581,7 +581,7 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
            __sum = *__varcovbegin++ - __sum;
            if (__builtin_expect(__sum <= _RealType(0), 0))
              std::__throw_runtime_error(__N("normal_mv_distribution::"
-                                            "param_type::_M_init_full"));
+                                            "param_type::_M_init_lower"));
            *__w++ = std::sqrt(__sum);
          }
       }
@@ -709,9 +709,11 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
 
       __is >> __x._M_nd;
 
+      // The param_type temporary is built with a private constructor,
+      // to skip the Cholesky decomposition that would be performed
+      // otherwise.
       __x.param(typename normal_mv_distribution<_Dimen, _RealType>::
-               param_type(__mean.begin(), __mean.end(),
-                          __varcov.begin(), __varcov.end()));
+               param_type(__mean, __varcov));
 
       __is.flags(__flags);
       return __is;
index b05417156d191adff40dffad99bef6352859bcf4..a4925fc1c41be34d88f71d578145db6d4c7e950b 100644 (file)
@@ -23,6 +23,7 @@
 
 #include <ext/random>
 #include <sstream>
+#include <testsuite_hooks.h>
 
 void
 test01()
@@ -35,6 +36,7 @@ test01()
   str << u;
 
   str >> v;
+  VERIFY( u == v );
 }
 
 int main()
index 8d83f9e6966d2ad49d0e5a64e73e4a181261c94d..f5fbc42a686f09f66833f8d06be05c8b3b215b2d 100644 (file)
@@ -23,6 +23,7 @@
 
 #include <ext/random>
 #include <sstream>
+#include <testsuite_hooks.h>
 
 void
 test01()
@@ -35,6 +36,7 @@ test01()
   str << u;
 
   str >> v;
+  VERIFY( u == v );
 }
 
 int main()
index cf17fea8b03ffd7a065b26ef4370b1aa02152c33..75e16cf0437cbd8e4b77c0c505e93825f55cfc2e 100644 (file)
@@ -23,6 +23,7 @@
 
 #include <ext/random>
 #include <sstream>
+#include <testsuite_hooks.h>
 
 void
 test01()
@@ -35,6 +36,7 @@ test01()
   str << u;
 
   str >> v;
+  VERIFY( u == v );
 }
 
 int main()
index f3d7912e314bac49355d762f2a3f2df15fb50037..b32a31dee64212ba342e4e005c14dd4642675144 100644 (file)
@@ -23,6 +23,7 @@
 
 #include <ext/random>
 #include <sstream>
+#include <testsuite_hooks.h>
 
 void
 test01()
@@ -35,6 +36,7 @@ test01()
   str << u;
 
   str >> v;
+  VERIFY( u == v );
 }
 
 int main()