MezzanineEngine 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cubicspline.h
Go to the documentation of this file.
1 // © Copyright 2010 - 2014 BlackTopp Studios Inc.
2 /* This file is part of The Mezzanine Engine.
3 
4  The Mezzanine Engine is free software: you can redistribute it and/or modify
5  it under the terms of the GNU General Public License as published by
6  the Free Software Foundation, either version 3 of the License, or
7  (at your option) any later version.
8 
9  The Mezzanine Engine is distributed in the hope that it will be useful,
10  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  GNU General Public License for more details.
13 
14  You should have received a copy of the GNU General Public License
15  along with The Mezzanine Engine. If not, see <http://www.gnu.org/licenses/>.
16 */
17 /* The original authors have included a copy of the license specified above in the
18  'Docs' folder. See 'gpl.txt'
19 */
20 /* We welcome the use of the Mezzanine engine to anyone, including companies who wish to
21  Build professional software and charge for their product.
22 
23  However there are some practical restrictions, so if your project involves
24  any of the following you should contact us and we will try to work something
25  out:
26  - DRM or Copy Protection of any kind(except Copyrights)
27  - Software Patents You Do Not Wish to Freely License
28  - Any Kind of Linking to Non-GPL licensed Works
29  - Are Currently In Violation of Another Copyright Holder's GPL License
30  - If You want to change our code and not add a few hundred MB of stuff to
31  your distribution
32 
33  These and other limitations could cause serious legal problems if you ignore
34  them, so it is best to simply contact us or the Free Software Foundation, if
35  you have any questions.
36 
37  Joseph Toppi - toppij@gmail.com
38  John Blackwood - makoenergy02@gmail.com
39 */
40 #ifndef _cubicspline_h
41 #define _cubicspline_h
42 
43 /// @file
44 /// @brief
45 
46 #include "datatypes.h"
47 
48 namespace Mezzanine
49 {
50  // The following license applies only to the Element and Spline class included below.
51  /* "THE BEER-WARE LICENSE" (Revision 42): Devin Lane wrote this file. As long as you retain
52  * this notice you can do whatever you want with this stuff. If we meet some day, and you
53  * think this stuff is worth it, you can buy me a beer in return. */
54  /// @brief Meta data nodes required to generate a cubic spline
55  template <typename X, typename Y>
57  {
58  public:
59 
60  /// @brief Time series meta-data
61  X x;
62 
63  /// @brief Meta-data corresponding to "a" on wikipedia description of a cubic splines
64  Y a;
65  /// @brief Meta-data corresponding to "b" on wikipedia description of a cubic splines
66  Y b;
67  /// @brief Meta-data corresponding to "c" on wikipedia description of a cubic splines
68  Y c;
69  /// @brief Meta-data corresponding to "d" on wikipedia description of a cubic splines
70  Y d;
71 
72  /// @brief Simple constructor
73  /// @param _x
74  CubicSplineElement(X _x) : x(_x) {}
75 
76  /// @brief Complete meta-data constructor
77  /// _x Time series location
78  /// _a Meta-data corresponding to "a" on wikipedia description of a cubic splines
79  /// _b Meta-data corresponding to "b" on wikipedia description of a cubic splines
80  /// _c Meta-data corresponding to "c" on wikipedia description of a cubic splines
81  /// _d Meta-data corresponding to "d" on wikipedia description of a cubic splines
82  CubicSplineElement(X _x, Y _a, Y _b, Y _c, Y _d)
83  : x(_x), a(_a), b(_b), c(_c), d(_d)
84  {}
85 
86  /// @brief The actual interpolation of the meta-data into a result
87  /// @param xx Location in the time series toget the corresponding Data series location on a cubic curve.
88  Y eval(const X& xx) const
89  {
90  X xix(xx - x);
91  return a + b * xix + c * (xix * xix) + d * (xix * xix * xix);
92  }
93 
94  /// @brief Sort Meta-data elements by time series locations in other Meta-Data
95  /// @param e The CubicSplineElement to compare this one too.
96  bool operator<(const CubicSplineElement& e) const
97  { return x < e.x; }
98 
99  /// @brief Sort Meta-data elements by time series locations
100  /// @param xx The CubicSplineElement to compare this one too.
101  bool operator<(const X& xx) const
102  { return x < xx; }
103  };
104 
105  /// @brief A class for interpolating data with arbitrary
106  /// @details Templated on type of X, Y. X and Y must have operator +, -, *, /. Y must have defined
107  /// a constructor that takes a scalar.
108  template <typename TimeSeriesType, typename DataSeriesType>
110  {
111  public:
112  /// @brief The meta data type stored by this container.
114  /// @brief The underlying type of the container that actually stores the instances of ContainerType.
115  typedef std::vector<element_type> ContainerType;
116  /// @brief An iterator suitable for the metadata this spline uses.
117  typedef typename ContainerType::iterator iterator;
118  /// @brief An iterator suitable for the metadata this spline uses that will not allow changes to the underlying data.
119  typedef typename ContainerType::const_iterator const_iterator;
120 
121  /// @brief The actual container of Metadata.
123 
124  /// @brief Create an empty, invalid spline
126 
127  /// @brief Create a Spline with a time series spread evenly between 0 and 1.
128  /// @param y The data series.
129  CubicSpline(const std::vector<DataSeriesType>& y)
130  {
131  std::vector<TimeSeriesType> x;
132 
133 
134  CalculateElements(x,y);
135  }
136 
137  /// @brief A spline with custom time and data series values
138  CubicSpline(const std::vector<TimeSeriesType>& x, const std::vector<DataSeriesType>& y)
139  { CalculateElements(x,y); }
140 
141  virtual ~CubicSpline() {}
142 
143  DataSeriesType operator[](const TimeSeriesType& x) const
144  { return interpolate(x); }
145 
146  DataSeriesType interpolate(const TimeSeriesType&x) const
147  {
148  if (DataElements.size() == 0)
149  { return DataSeriesType(); }
150  typename std::vector<element_type>::const_iterator it;
151  it = std::lower_bound(DataElements.begin(), DataElements.end(), element_type(x));
152  if (it != DataElements.begin())
153  { it--; }
154  return it->eval(x);
155  }
156 
157  std::vector<DataSeriesType> operator[](const std::vector<TimeSeriesType>& xx) const
158  { return interpolate(xx); }
159 
160  iterator begin()
161  { return DataElements.begin(); }
162  iterator end()
163  { return DataElements.end(); }
164 
165  const_iterator begin() const
166  { return DataElements.begin(); }
167  const_iterator end() const
168  { return DataElements.end(); }
169 
170  /// @brief Add another entry to the spline
171  /// @details Adjust the time series to evenly distribute it between 0 and 1 for
172  /// each entry, including the new one DataToAdd
173  void push_back(const DataSeriesType& DataToAdd)
174  {
175  // This needs to be implemented still.
176  std::vector<TimeSeriesType> x;
177  std::vector<DataSeriesType> y;
178 
179 
180  CalculateElements(x,y);
181  }
182 
183  /* Evaluate at multiple locations, assuming xx is sorted ascending */
184  std::vector<DataSeriesType> interpolate(const std::vector<TimeSeriesType>& xx) const
185  {
186  if (DataElements.size() == 0)
187  { return std::vector<DataSeriesType>(xx.size()); }
188 
189  typename std::vector<TimeSeriesType>::const_iterator it;
190  typename std::vector<element_type>::const_iterator it2;
191  it2 = DataElements.begin();
192  std::vector<DataSeriesType> ys;
193  for (it = xx.begin(); it != xx.end(); it++)
194  {
195  it2 = std::lower_bound(it2, DataElements.end(), element_type(*it));
196  if (it2 != DataElements.begin())
197  { it2--; }
198 
199  ys.push_back(it2->eval(*it));
200  }
201  return ys;
202  }
203 
204  protected:
205 
206  void CalculateElements(const std::vector<TimeSeriesType>& x, const std::vector<DataSeriesType>& y)
207  {
208  if (x.size() != y.size())
209  { MEZZ_EXCEPTION(Exception::PARAMETERS_RANGE_EXCEPTION,"Data series and time series must have the same count of elements."); }
210 
211  if (x.size() < 3)
212  { MEZZ_EXCEPTION(Exception::PARAMETERS_RANGE_EXCEPTION,"Must have at least three points for interpolation."); }
213 
214  DataElements.clear();
215 
216  typedef typename std::vector<TimeSeriesType>::difference_type size_type;
217 
218  size_type n = y.size() - 1;
219 
220  std::vector<DataSeriesType> b(n), d(n), a(n), c(n+1), l(n+1), u(n+1), z(n+1);
221  std::vector<TimeSeriesType> h(n+1);
222 
223  l[0] = DataSeriesType(1);
224  u[0] = DataSeriesType(0);
225  z[0] = DataSeriesType(0);
226  h[0] = x[1] - x[0];
227 
228  for (size_type i = 1; i < n; i++)
229  {
230  h[i] = x[i+1] - x[i];
231  l[i] = DataSeriesType(2 * (x[i+1] - x[i-1])) - DataSeriesType(h[i-1]) * u[i-1];
232  u[i] = DataSeriesType(h[i]) / l[i];
233  a[i] = (DataSeriesType(3) / DataSeriesType(h[i])) * (y[i+1] - y[i]) - (DataSeriesType(3) / DataSeriesType(h[i-1])) * (y[i] - y[i-1]);
234  z[i] = (a[i] - DataSeriesType(h[i-1]) * z[i-1]) / l[i];
235  }
236 
237  l[n] = DataSeriesType(1);
238  z[n] = c[n] = DataSeriesType(0);
239 
240  for (size_type j = n-1; j >= 0; j--)
241  {
242  c[j] = z[j] - u[j] * c[j+1];
243  b[j] = (y[j+1] - y[j]) / DataSeriesType(h[j]) - (DataSeriesType(h[j]) * (c[j+1] + DataSeriesType(2) * c[j])) / DataSeriesType(3);
244  d[j] = (c[j+1] - c[j]) / DataSeriesType(3 * h[j]);
245  }
246 
247  for (size_type i = 0; i < n; i++)
248  {
249  DataElements.push_back(element_type(x[i], y[i], b[i], c[i], d[i]));
250  }
251  }
252  };
253 
254 } // /namespace Mezzanine
255 
256 #endif // Include guard