smooth_linear_spline.cpp
Go to the documentation of this file.
1 
11 /*****************************************************************************
12 ** Includes
13 *****************************************************************************/
14 
15 #include <iostream>
16 #include <cmath>
17 #include <ecl/exceptions/standard_exception.hpp>
18 #include "../../include/ecl/geometry/smooth_linear_spline.hpp"
19 
20 /*****************************************************************************
21 ** Namespaces
22 *****************************************************************************/
23 
24 namespace ecl {
25 
26 /*****************************************************************************
27 ** Implementation
28 *****************************************************************************/
29 
30 SmoothLinearSpline::SmoothLinearSpline(const Array<double>& x_data, const Array<double>& y_data, double a_max) throw (DataException<int>) {
31 
32  if ( x_data.size() != y_data.size() ) {
33  throw DataException<int>(LOC,OutOfRangeError,"Input domain and range sets were not the same size.", 0);
34  }
35 
36  /*********************
37  ** Basic Parameters
38  **********************/
39  // n = number of linear segments
40  // 2*n = number of points (2 end points + 2*(n-1) corners) remember we're including the polynomial corners in this calculation
41  unsigned int n = x_data.size()-1;
42 
43  /******************************************
44  ** Storage
45  *******************************************/
46  segments.resize(n);
47  corners.resize(n-1);
48  discretised_domain.resize(2*n);
49 
50  /******************************************
51  ** Build Segments
52  *******************************************/
53  for ( unsigned int i = 0; i < n; ++i ) {
54  segments[i] = LinearFunction::Interpolation(x_data[i],y_data[i],x_data[i+1],y_data[i+1]);
55  }
56  discretised_domain[0] = x_data[0]; // Initial point
57  discretised_domain[2*n-1] = x_data[n]; // Last point
58 
59  /******************************************
60  ** Build Corners (this might throw)
61  *******************************************/
62  for ( unsigned int i = 1; i < n; ++i ) {
63  if ( segments[i-1].coefficients()[1] == segments[i].coefficients()[1] ) { // slope identical
64  // Just set the quintic to be a linear function.
65  corners[i-1].coefficients() << segments[i-1].coefficients()[0], segments[i-1].coefficients()[1], 0.0, 0.0, 0.0, 0.0;
66  discretised_domain[1+2*(i-1)] = x_data[i];
67  discretised_domain[1+2*(i-1)+1] = x_data[i];
68  } else {
69  for ( unsigned int j = 1; j <= 5; ++j ) {
70  double x_l, x_r;
71  if ( i == 1 ) { // First one, we can use the whole segment
72  x_l = x_data[i] - j*(x_data[i]-x_data[i-1])/5;
73  } else {
74  x_l = x_data[i] - j*(x_data[i]-x_data[i-1])/10;
75  }
76  if ( i == (n-1) ) { // Last one, we can use the whole segment
77  x_r = x_data[i] + j*(x_data[i+1]-x_data[i])/5;
78  } else {
79  x_r = x_data[i] + j*(x_data[i+1]-x_data[i])/10;
80  }
81  double y_l = segments[i-1](x_l);
82  double y_r = segments[i](x_r);
83  double ydot_l = segments[i-1].derivative(x_l);
84  double ydot_r = segments[i].derivative(x_r);
85  corners[i-1] = QuinticPolynomial::Interpolation(x_l,y_l,ydot_l,0.0,
86  x_r,y_r,ydot_r,0.0);
87  if ( ( fabs( Maximum<CubicPolynomial>()(x_l,x_r,corners[i-1].derivative().derivative())) < fabs(a_max) ) &&
88  ( fabs( Minimum<CubicPolynomial>()(x_l,x_r,corners[i-1].derivative().derivative())) < fabs(a_max) ) ) {
89  discretised_domain[1+2*(i-1)] = x_l;
90  discretised_domain[1+2*(i-1)+1] = x_r;
91  break;
92  }
93  /*********************
94  ** Debugging
95  **********************/
96 // if ( i == 1 ) {
97 // double max = fabs( Maximum<CubicPolynomial>()(x_l,x_r,corners[i-1].derivative().derivative()));
98 // if ( fabs( Minimum<CubicPolynomial>()(x_l,x_r,corners[i-1].derivative().derivative())) > max ) {
99 // max = fabs( Minimum<CubicPolynomial>()(x_l,x_r,corners[i-1].derivative().derivative()));
100 // }
101 // std::cout << "Max: " << max << std::endl;
102 // }
103 
104  if ( j == 5 ) { // Cannot exceed more than half the length of the segment.
105  throw DataException<int>(LOC,ConstructorError,"Max acceleration bound could not be satisfied at corner specified by the data element.",i);
106  }
107  }
108  }
109  }
110 }
111 
112 double SmoothLinearSpline::operator()(const double &x) const ecl_assert_throw_decl(StandardException) {
113  ecl_assert_throw( ( ( x >= discretised_domain.front() ) && ( x <= discretised_domain.back() ) ), StandardException(LOC,OutOfRangeError) );
114  int index = 0;
115  while ( x > discretised_domain[index+1] ) {
116  ++index;
117  }
118  if ( index % 2 == 0 ) { // linear
119  return segments[index/2](x);
120  } else { // quintic
121  return corners[(index-1)/2](x);
122  }
123 }
124 
125 double SmoothLinearSpline::derivative(const double &x) const ecl_assert_throw_decl(StandardException) {
126  ecl_assert_throw( ( ( x >= discretised_domain.front() ) && ( x <= discretised_domain.back() ) ), StandardException(LOC,OutOfRangeError) );
127  int index = 0;
128  while ( x > discretised_domain[index+1] ) {
129  ++index;
130  }
131  if ( index % 2 == 0 ) { // linear
132  return segments[index/2].derivative(x);
133  } else { // quintic
134  return corners[(index-1)/2].derivative(x);
135  }
136 }
137 
138 double SmoothLinearSpline::dderivative(const double &x) const ecl_assert_throw_decl(StandardException) {
139  ecl_assert_throw( ( ( x >= discretised_domain.front() ) && ( x <= discretised_domain.back() ) ), StandardException(LOC,OutOfRangeError) );
140  int index = 0;
141  while ( x > discretised_domain[index+1] ) {
142  ++index;
143  }
144  if ( index % 2 == 0 ) { // linear
145  return segments[index/2].dderivative(x);
146  } else { // quintic
147  return corners[(index-1)/2].dderivative(x);
148  }
149 }
150 
151 } // namespace ecl
152 
153 
154 
155 
ecl::ConstructorError
ConstructorError
ecl
Embedded control libraries.
ecl::SmoothLinearSpline::dderivative
double dderivative(const double &x) const ecl_assert_throw_decl(StandardException)
Spline second derivative.
Definition: smooth_linear_spline.cpp:138
ecl::Minimum< CubicPolynomial >
Mathematical minimum on a compact interval for cubic polynomials.
Definition: polynomial.hpp:1063
LOC
#define LOC
ecl_assert_throw
#define ecl_assert_throw(expression, exception)
ecl::OutOfRangeError
OutOfRangeError
ecl::SmoothLinearSpline::operator()
double operator()(const double &x) const ecl_assert_throw_decl(StandardException)
Spline function.
Definition: smooth_linear_spline.cpp:112
ecl::Maximum< CubicPolynomial >
Mathematical maximum on a compact interval for cubic polynomials.
Definition: polynomial.hpp:1014
ecl_assert_throw_decl
#define ecl_assert_throw_decl(exception)
ecl::SmoothLinearSpline::SmoothLinearSpline
SmoothLinearSpline()
Default constructor.
Definition: smooth_linear_spline.hpp:53
standard_exception.hpp
ecl::SmoothLinearSpline::derivative
double derivative(const double &x) const ecl_assert_throw_decl(StandardException)
Spline derivative.
Definition: smooth_linear_spline.cpp:125

ecl_geometry
Author(s): Daniel Stonier
autogenerated on Mon Jun 10 2019 13:08:37