From 50e1643e4effc0df73cbcf17416a68139a72586a Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Thu, 9 May 2019 07:18:18 +0200 Subject: [PATCH] add Utilities::MPI::mean_and_standard_deviation() --- doc/news/changes/minor/20190509DenisDavydov | 3 + include/deal.II/base/mpi.h | 55 +++++++ tests/mpi/collective_mean_std.cc | 148 ++++++++++++++++++ tests/mpi/collective_mean_std.mpirun=2.output | 5 + 4 files changed, 211 insertions(+) create mode 100644 doc/news/changes/minor/20190509DenisDavydov create mode 100644 tests/mpi/collective_mean_std.cc create mode 100644 tests/mpi/collective_mean_std.mpirun=2.output diff --git a/doc/news/changes/minor/20190509DenisDavydov b/doc/news/changes/minor/20190509DenisDavydov new file mode 100644 index 0000000000..3ff2b2ceb6 --- /dev/null +++ b/doc/news/changes/minor/20190509DenisDavydov @@ -0,0 +1,3 @@ +New: Add Utilities::MPI::mean_and_standard_deviation() to calculate mean and standard deviation. +
+(Denis Davydov, 2019/05/09) diff --git a/include/deal.II/base/mpi.h b/include/deal.II/base/mpi.h index 74bef20899..99f70cce59 100644 --- a/include/deal.II/base/mpi.h +++ b/include/deal.II/base/mpi.h @@ -19,8 +19,10 @@ #include #include +#include #include +#include #include #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 + std::pair::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 + std::pair::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::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(size); + Std sq_sum = 0.; + std::for_each(begin, end, [&mean, &sq_sum](const Number &v) { + sq_sum += numbers::NumberTraits::abs_square(v - mean); + }); + sq_sum = Utilities::MPI::sum(sq_sum, comm); + return std::make_pair(mean, + std::sqrt(sq_sum / static_cast(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 index 0000000000..2cb17d827a --- /dev/null +++ b/tests/mpi/collective_mean_std.cc @@ -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 + +#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 values; + std::vector> val_c; + std::vector same; + std::vector 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>( + 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 index 0000000000..323753c8c5 --- /dev/null +++ b/tests/mpi/collective_mean_std.mpirun=2.output @@ -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 -- 2.39.5