#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)
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
# 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+}
projects / dealii.git / commitdiff