Exotica
incremental_gaussian.h
Go to the documentation of this file.
1 //
2 // Copyright (c) 2018, University of Edinburgh
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without
6 // modification, are permitted provided that the following conditions are met:
7 //
8 // * Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright
11 // notice, this list of conditions and the following disclaimer in the
12 // documentation and/or other materials provided with the distribution.
13 // * Neither the name of nor the names of its contributors may be used to
14 // endorse or promote products derived from this software without specific
15 // prior written permission.
16 //
17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 // POSSIBILITY OF SUCH DAMAGE.
28 //
29 
30 // The algorithm is based on the paper :
31 // Chan, Tony F.; Golub, Gene H.; LeVeque, Randall J. (1979), “Updating Formulae and a Pairwise Algorithm for Computing Sample Variances.”, Technical Report STAN-CS-79-773, Department of Computer Science, Stanford University.
32 
33 #ifndef EXOTICA_AICO_SOLVER_INCREMENTAL_GAUSSIAN_H_
34 #define EXOTICA_AICO_SOLVER_INCREMENTAL_GAUSSIAN_H_
35 
36 #include <Eigen/Dense>
37 #include <vector>
38 
39 namespace exotica
40 {
42 {
43 private:
44  int D = 0;
45  double W = 0;
46  Eigen::VectorXd T = Eigen::VectorXd(0);
47  Eigen::VectorXd dX = Eigen::VectorXd(0);
48  Eigen::MatrixXd S = Eigen::MatrixXd(0, 0);
49 
50 public:
51  SinglePassMeanCovariance() = default;
52 
54  {
55  resize(D_);
56  }
57 
58  void resize(int D_)
59  {
60  D = D_;
61 
62  T.resize(D_);
63  dX.resize(D_);
64  S.resize(D_, D_);
65 
66  clear();
67  }
68 
69  void clear()
70  {
71  T.setZero();
72  dX.setZero();
73  S.setZero();
74  W = 0.;
75  }
76 
77  void add(const Eigen::Ref<const Eigen::VectorXd>& x)
78  {
79  if (W == 0.)
80  {
81  W = 1.;
82  T = x;
83  S.setZero();
84  return;
85  }
86  W += 1.;
87  T += x;
88  double f = 1. / W / (W - 1.);
89  dX = W * x - T;
90  for (int r = 0; r < D; ++r)
91  {
92  for (int c = 0; c < D; ++c)
93  {
94  S(r, c) += f * dX(r) * dX(c);
95  }
96  }
97  }
98 
99  inline void add(SinglePassMeanCovariance& M)
100  {
101  add(M.W, M.T, M.S);
102  }
103 
104  void add(double& W_, const Eigen::Ref<const Eigen::VectorXd>& T_,
105  const Eigen::Ref<const Eigen::VectorXd>& S_)
106  {
107  if (W == 0.)
108  {
109  W = W_;
110  T = T_;
111  S = S_;
112  return;
113  }
114  dX = T_ / W_ - T / W;
115 
116  double f = W * W_ / (W + W_);
117  for (int r = 0; r < D; ++r)
118  {
119  for (int c = 0; c < D; ++c)
120  {
121  S(r, c) += S_(r, c) + f * dX(r) * dX(c);
122  }
123  }
124  T += T_;
125  W += W_;
126  }
127 
128  inline void addw(double w, const Eigen::Ref<const Eigen::VectorXd>& x)
129  {
130  if (W == 0.)
131  {
132  W = w;
133  T = w * x;
134  S.setZero();
135  return;
136  }
137 
138  dX = x - T / W;
139 
140  double f = W * w / (W + w);
141  for (int r = 0; r < D; ++r)
142  {
143  for (int c = 0; c < D; ++c)
144  {
145  S(r, c) += f * dX(r) * dX(c);
146  }
147  }
148 
149  T += w * x;
150  W += w;
151  }
152 
153  void mean(Eigen::VectorXd& mu)
154  {
155  mu = T / W;
156  }
157  void cov(Eigen::MatrixXd& sig)
158  {
159  sig = S / W;
160  }
161  void covp(Eigen::MatrixXd& sig)
162  {
163  sig = S / (W - 1.);
164  }
165 };
166 } // namespace exotica
167 
168 #endif // EXOTICA_AICO_SOLVER_INCREMENTAL_GAUSSIAN_H_
exotica::SinglePassMeanCovariance::SinglePassMeanCovariance
SinglePassMeanCovariance()=default
exotica::SinglePassMeanCovariance::add
void add(SinglePassMeanCovariance &M)
Definition: incremental_gaussian.h:99
exotica::SinglePassMeanCovariance::W
double W
Definition: incremental_gaussian.h:45
exotica::SinglePassMeanCovariance::add
void add(double &W_, const Eigen::Ref< const Eigen::VectorXd > &T_, const Eigen::Ref< const Eigen::VectorXd > &S_)
Definition: incremental_gaussian.h:104
exotica::SinglePassMeanCovariance::clear
void clear()
Definition: incremental_gaussian.h:69
exotica::SinglePassMeanCovariance::T
Eigen::VectorXd T
Definition: incremental_gaussian.h:46
exotica
Definition: cartpole_dynamics_solver.h:38
exotica::SinglePassMeanCovariance::covp
void covp(Eigen::MatrixXd &sig)
Definition: incremental_gaussian.h:161
exotica::SinglePassMeanCovariance::cov
void cov(Eigen::MatrixXd &sig)
Definition: incremental_gaussian.h:157
exotica::SinglePassMeanCovariance::mean
void mean(Eigen::VectorXd &mu)
Definition: incremental_gaussian.h:153
exotica::SinglePassMeanCovariance::resize
void resize(int D_)
Definition: incremental_gaussian.h:58
exotica::SinglePassMeanCovariance::D
int D
Definition: incremental_gaussian.h:44
exotica::SinglePassMeanCovariance::addw
void addw(double w, const Eigen::Ref< const Eigen::VectorXd > &x)
Definition: incremental_gaussian.h:128
exotica::SinglePassMeanCovariance::dX
Eigen::VectorXd dX
Definition: incremental_gaussian.h:47
exotica::SinglePassMeanCovariance
Definition: incremental_gaussian.h:41
exotica::SinglePassMeanCovariance::S
Eigen::MatrixXd S
Definition: incremental_gaussian.h:48
exotica::SinglePassMeanCovariance::SinglePassMeanCovariance
SinglePassMeanCovariance(int D_)
Definition: incremental_gaussian.h:53
exotica::SinglePassMeanCovariance::add
void add(const Eigen::Ref< const Eigen::VectorXd > &x)
Definition: incremental_gaussian.h:77