]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add Utilities::MPI::mean_and_standard_deviation() 8048/head
authorDenis Davydov <davydden@gmail.com>
Thu, 9 May 2019 05:18:18 +0000 (07:18 +0200)
committerDenis Davydov <davydden@gmail.com>
Thu, 9 May 2019 15:02:28 +0000 (17:02 +0200)
doc/news/changes/minor/20190509DenisDavydov [new file with mode: 0644]
include/deal.II/base/mpi.h
tests/mpi/collective_mean_std.cc [new file with mode: 0644]
tests/mpi/collective_mean_std.mpirun=2.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20190509DenisDavydov b/doc/news/changes/minor/20190509DenisDavydov
new file mode 100644 (file)
index 0000000..3ff2b2c
--- /dev/null
@@ -0,0 +1,3 @@
+New: Add Utilities::MPI::mean_and_standard_deviation() to calculate mean and standard deviation.
+<br>
+(Denis Davydov, 2019/05/09)
index 74bef20899e633766b28b6ba7c4d41e1938a1d6f..99f70cce59d165199d2d1058ffcd693c248cbbdf 100644 (file)
 #include <deal.II/base/config.h>
 
 #include <deal.II/base/array_view.h>
+#include <deal.II/base/numbers.h>
 
 #include <map>
+#include <numeric>
 #include <vector>
 
 #if !defined(DEAL_II_WITH_MPI) && !defined(DEAL_II_WITH_PETSC)
@@ -238,6 +240,29 @@ namespace Utilities
     create_ascending_partitioning(const MPI_Comm &           comm,
                                   const IndexSet::size_type &local_size);
 
+#ifdef DEAL_II_WITH_MPI
+    /**
+     * Calculate mean and standard deviation across the MPI communicator @p comm
+     * for values provided as a range `[begin,end)`.
+     * The mean is computed as $\bar x=\frac 1N \sum x_k$ where the $x_k$ are
+     * the elements pointed to by the `begin` and `end` iterators on all
+     * processors (i.e., each processor's `[begin,end)` range points to a subset
+     * of the overall number of elements). The standard deviation is calculated
+     * as $\sigma=\sqrt{\frac {1}{N-1} \sum |x_k -\bar x|^2}$, which is known as
+     * unbiased sample variance.
+     *
+     * @tparam Number specifies the type to store the mean value.
+     * The standard deviation is stored as the corresponding real type.
+     * This allows, for example, to calculate statistics from integer input
+     * values.
+     */
+    template <class Iterator, typename Number = long double>
+    std::pair<Number, typename numbers::NumberTraits<Number>::real_type>
+    mean_and_standard_deviation(const Iterator  begin,
+                                const Iterator  end,
+                                const MPI_Comm &comm);
+#endif
+
     /**
      * Return the sum over all processors of the value @p t. This function is
      * collective over all processors given in the
@@ -952,6 +977,36 @@ namespace Utilities
 #  endif
     }
 
+
+#  ifdef DEAL_II_WITH_MPI
+    template <class Iterator, typename Number>
+    std::pair<Number, typename numbers::NumberTraits<Number>::real_type>
+    mean_and_standard_deviation(const Iterator  begin,
+                                const Iterator  end,
+                                const MPI_Comm &comm)
+    {
+      // below we do simple and straight-forward implementation. More elaborate
+      // options are:
+      // http://dx.doi.org/10.1145/2807591.2807644 section 3.1.2
+      // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm
+      // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online
+      using Std        = typename numbers::NumberTraits<Number>::real_type;
+      const Number sum = std::accumulate(begin, end, Number(0.));
+
+      const auto size = Utilities::MPI::sum(std::distance(begin, end), comm);
+      Assert(size > 0, ExcDivideByZero());
+      const Number mean =
+        Utilities::MPI::sum(sum, comm) / static_cast<Std>(size);
+      Std sq_sum = 0.;
+      std::for_each(begin, end, [&mean, &sq_sum](const Number &v) {
+        sq_sum += numbers::NumberTraits<Number>::abs_square(v - mean);
+      });
+      sq_sum = Utilities::MPI::sum(sq_sum, comm);
+      return std::make_pair(mean,
+                            std::sqrt(sq_sum / static_cast<Std>(size - 1)));
+    }
+#  endif
+
 #endif
   } // end of namespace MPI
 } // end of namespace Utilities
diff --git a/tests/mpi/collective_mean_std.cc b/tests/mpi/collective_mean_std.cc
new file mode 100644 (file)
index 0000000..2cb17d8
--- /dev/null
@@ -0,0 +1,148 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2009 - 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check Utilities::MPI::mean_and_standard_deviation()
+
+#include <deal.II/base/mpi.h>
+
+#include "../tests.h"
+
+/* MWE in Python
+>>> import numpy as np
+
+>>> a = np.array([-3, 55, 0, 1, 11, -12])
+>>> np.mean(a)
+8.6666666666666661
+>>> np.std(a, ddof=1)
+23.871880249923059
+
+>>> b = np.array([-1+2.j, 3+7.j, -5.j,6])
+>>> np.mean(b)
+(2+1j)
+>>> np.std(b,ddof=1)
+5.8878405775518976
+
+>>> c = np.array([1,1,1])
+>>> np.mean(c)
+1.0
+>>> np.std(c,ddof=1)
+0.0
+
+>>> d = np.array([1,2])
+>>> np.mean(d)
+1.5
+>>> np.std(d,ddof=1)
+0.70710678118654757
+*/
+
+void
+test()
+{
+  unsigned int     myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  std::vector<int> values;
+  std::vector<std::complex<double>> val_c;
+  std::vector<int>                  same;
+  std::vector<int>                  empty;
+  if (myid == 0)
+    {
+      values.resize(4);
+      values[0] = -3;
+      values[1] = 55;
+      values[2] = 0;
+      values[3] = 1;
+
+      val_c.resize(1);
+      val_c[0].real(-1);
+      val_c[0].imag(2.);
+
+      same.resize(2);
+      same[0] = 1;
+      same[1] = 1;
+
+      empty.resize(2);
+      empty[0] = 1;
+      empty[1] = 2;
+    }
+  else
+    {
+      values.resize(2);
+      values[0] = 11;
+      values[1] = -12;
+
+      val_c.resize(3);
+      val_c[0].real(3);
+      val_c[0].imag(7);
+      val_c[1].real(0);
+      val_c[1].imag(-5);
+      val_c[2].real(6);
+      val_c[2].imag(0);
+
+      same.resize(1);
+      same[0] = 1;
+    }
+
+  const auto pair = Utilities::MPI::mean_and_standard_deviation(values.begin(),
+                                                                values.end(),
+                                                                MPI_COMM_WORLD);
+
+  const auto pair_c =
+    Utilities::MPI::mean_and_standard_deviation<decltype(val_c.begin()),
+                                                std::complex<double>>(
+      val_c.begin(), val_c.end(), MPI_COMM_WORLD);
+
+  const auto pair_same =
+    Utilities::MPI::mean_and_standard_deviation(same.begin(),
+                                                same.end(),
+                                                MPI_COMM_WORLD);
+
+  const auto pair_empty =
+    Utilities::MPI::mean_and_standard_deviation(empty.begin(),
+                                                empty.end(),
+                                                MPI_COMM_WORLD);
+
+  if (myid == 0)
+    deallog << pair.first << ' ' << pair.second << std::endl
+            << pair_c.first << ' ' << pair_c.second << std::endl
+            << pair_same.first << ' ' << pair_same.second << std::endl
+            << pair_empty.first << ' ' << pair_empty.second << std::endl;
+}
+
+
+int
+main(int argc, char *argv[])
+{
+#ifdef DEAL_II_WITH_MPI
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+#else
+  (void)argc;
+  (void)argv;
+  compile_time_error;
+
+#endif
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    {
+      initlog();
+
+      deallog.push("mpi");
+      test();
+      deallog.pop();
+    }
+  else
+    test();
+}
diff --git a/tests/mpi/collective_mean_std.mpirun=2.output b/tests/mpi/collective_mean_std.mpirun=2.output
new file mode 100644 (file)
index 0000000..323753c
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL:mpi::8.66667 23.8719
+DEAL:mpi::(2.00000,1.00000) 5.88784
+DEAL:mpi::1.00000 0.00000
+DEAL:mpi::1.50000 0.707107

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.