Averager.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <boost/none.hpp>
8 #include <boost/optional.hpp>
9 #include <cstddef>
10 #include <deque>
11 
12 #include "DataStructures/DataVector.hpp"
13 
14 /// \ingroup ControlSystemGroup
15 /// A weighted exponential averager of \f$Q\f$ and its derivatives
16 /// implementing Appendix A in \cite Hemberger2012jz
17 ///
18 /// The purpose of the averager is to provide \f$Q\f$ and 'smoothed' numerical
19 /// derivatives of \f$Q\f$ up to the `DerivOrder`'th derivative. The 0th
20 /// derivative of \f$Q\f$ is not typically averaged, since no differencing is
21 /// necessary and no noise is introduced. The `average_0th_deriv_of_q` option
22 /// allows for specifying that the 0th derivative should be averaged in addition
23 /// to the derivatives. This may be desirable for control systems where
24 /// \f$Q\f$ contains some 'noise', e.g. size control typically uses an averaged
25 /// \f$Q\f$ since the raw computed \f$Q\f$ is a function of the minimum on a
26 /// surface and may jump around discontinuously.
27 ///
28 /// The averager is designed to support an arbitrary DerivOrder, however the
29 /// private get_derivs method currently only supports DerivOrder = 2.
30 /// If an additional DerivOrder is needed, finite differencing needs to be
31 /// implemented to specifically handle that order (as it seems like overkill to
32 /// generalize the differencing stencil at this time).
33 template <size_t DerivOrder>
34 class Averager {
35  public:
36  /// `avg_timescale_frac` determines the exponential averaging timescale
37  /// through \f$\tau_\mathrm{avg} = \f$`avg_timescale_frac`\f$\times \tau\f$,
38  /// where \f$\tau\f$ is the damping time.
39  /// `avg_timescale_frac` must be positive.
40  /// `average_0th_deriv_of_q` determines whether the call operator returns an
41  /// averaged or unaveraged quantity for the 0th derivative of \f$Q\f$.
42  /// `true` returns an averaged 0th derivative of \f$Q\f$.
43  /// `false` returns the raw 0th derivative of \f$Q\f$.
44  /// The derivatives are always averaged (to reduce noise due to numerical
45  /// differentiation), so the `average_0th_deriv_of_q` option only specifies
46  /// whether to return an averaged value for the 0th derivative piece of the
47  /// function.
48  Averager(double avg_timescale_frac, bool average_0th_deriv_of_q) noexcept;
49 
50  // Explicitly defined move constructor due to the fact that the std::deque
51  // move constructor is not marked noexcept
52  Averager(Averager&& rhs) noexcept;
53  Averager& operator=(Averager&& rhs) noexcept;
54  Averager(const Averager&) = delete;
55  Averager& operator=(const Averager&) = delete;
56  ~Averager() = default;
57 
58  /// Returns \f$Q\f$ and its derivatives at \f$t=\f$`time`, provided there is
59  /// sufficient data. The averager is limited by the need for at least
60  /// (`DerivOrder` + 1) data points in order to provide the `DerivOrder`'th
61  /// derivative. If sufficient data is available, it returns \f$Q\f$ and
62  /// averaged derivatives of \f$Q\f$ up to the `DerivOrder`'th derivative, at
63  /// \f$t=\f$`time`. If `using_average_0th_deriv_of_q()` is `true`, then the
64  /// returned 0th derivative of \f$Q\f$ is also averaged. In the case that
65  /// there is insufficient data, the operator returns an
66  /// unitialized `boost::optional` (`boost::none`).
67  const boost::optional<std::array<DataVector, DerivOrder + 1>>& operator()(
68  double time) const noexcept;
69  /// A function that allows for resetting the averager.
70  void clear() noexcept;
71  /// The function responsible for updating the averaged values
72  /// at \f$t=\f$`time`. Requires `raw_q` (the raw components of \f$Q(t)\f$)
73  /// and `timescales` (the associated damping times for each component).
74  void update(double time, const DataVector& raw_q,
75  const DataVector& timescales) noexcept;
76  /// Returns the latest time at which the averager has sufficient data to
77  /// return \f$Q\f$ and its derivatives.
78  double last_time_updated() const noexcept;
79  /// Returns the exponentially averaged time at \f$t=\f$`time`. The time is
80  /// averaged along side \f$Q\f$ to determine the effective time at which
81  /// the average value is computed. The effective time is retarded, due to the
82  /// weighting of past times.
83  double average_time(double time) const noexcept;
84  /// Returns a bool corresponding to whether `average_0th_deriv_of_q`
85  /// is `true`/`false`.
86  bool using_average_0th_deriv_of_q() const noexcept {
87  return average_0th_deriv_of_q_;
88  };
89 
90  private:
91  /// Returns the function and numerical derivatives up to the
92  /// `DerivOrder`'th derivative using backward finite differencing that
93  /// supports non-uniform time spacing. Since we use a stencil size of
94  /// (`DerivOrder` + 1), the order of the finite difference varies for
95  /// different derivatives. Assuming \f$M\lte\f$`DerivOrder`, the \f$M\f$'th
96  /// derivative is of order (`DerivOrder` + 1 - \f$M\f$).
97  ///
98  /// NOTE: the derivative function currently only supports DerivOrder = 2,
99  /// and contains a static_assert to guard against the possible instantiation
100  /// of another DerivOrder, for which finite differencing does not
101  /// currently support.
102  std::array<DataVector, DerivOrder + 1> get_derivs() const noexcept;
103 
104  double avg_tscale_frac_;
105  bool average_0th_deriv_of_q_;
106  boost::optional<std::array<DataVector, DerivOrder + 1>> averaged_values_{};
107  boost::optional<std::array<DataVector, DerivOrder + 1>> boost_none_ =
108  boost::none;
109  std::deque<double> times_;
110  std::deque<DataVector> raw_qs_;
111  double weight_k_ = 0.0;
112  double tau_k_ = 0.0;
113 };
A weighted exponential averager of and its derivatives implementing Appendix A in ...
Definition: Averager.hpp:34
double average_time(double time) const noexcept
Returns the exponentially averaged time at time. The time is averaged along side to determine the ef...
Definition: Averager.cpp:152
void clear() noexcept
A function that allows for resetting the averager.
Definition: Averager.cpp:62
double last_time_updated() const noexcept
Returns the latest time at which the averager has sufficient data to return and its derivatives...
Definition: Averager.cpp:144
Averager(double avg_timescale_frac, bool average_0th_deriv_of_q) noexcept
avg_timescale_frac determines the exponential averaging timescale through avg_timescale_frac ...
Definition: Averager.cpp:14
bool using_average_0th_deriv_of_q() const noexcept
Returns a bool corresponding to whether average_0th_deriv_of_q is true/false.
Definition: Averager.hpp:86
void update(double time, const DataVector &raw_q, const DataVector &timescales) noexcept
The function responsible for updating the averaged values at time. Requires raw_q (the raw components...
Definition: Averager.cpp:71
Stores a collection of function values.
Definition: DataVector.hpp:46
const boost::optional< std::array< DataVector, DerivOrder+1 > > & operator()(double time) const noexcept
Returns and its derivatives at time, provided there is sufficient data. The averager is limited by t...
Definition: Averager.cpp:54