1
// Copyright 2012 The Closure Library Authors. All Rights Reserved.
2
//
3
// Licensed under the Apache License, Version 2.0 (the "License");
4
// you may not use this file except in compliance with the License.
5
// You may obtain a copy of the License at
6
//
7
//      http://www.apache.org/licenses/LICENSE-2.0
8
//
9
// Unless required by applicable law or agreed to in writing, software
10
// distributed under the License is distributed on an "AS-IS" BASIS,
11
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12
// See the License for the specific language governing permissions and
13
// limitations under the License.
14

15
/**
16
 * @fileoverview A one dimensional cubic spline interpolator with not-a-knot
17
 * boundary conditions.
18
 *
19
 * See http://en.wikipedia.org/wiki/Spline_interpolation.
20
 *
21
 */
22

23
goog.provide('goog.math.interpolator.Spline1');
24

25
goog.require('goog.array');
26
goog.require('goog.asserts');
27
goog.require('goog.math');
28
goog.require('goog.math.interpolator.Interpolator1');
29
goog.require('goog.math.tdma');
30

31

32

33
/**
34
 * A one dimensional cubic spline interpolator with natural boundary conditions.
35
 * @implements {goog.math.interpolator.Interpolator1}
36
 * @constructor
37
 */
38
goog.math.interpolator.Spline1 = function() {
39
  /**
40
   * The abscissa of the data points.
41
   * @type {!Array<number>}
42
   * @private
43
   */
44
  this.x_ = [];
45

46
  /**
47
   * The spline interval coefficients.
48
   * Note that, in general, the length of coeffs and x is not the same.
49
   * @type {!Array<!Array<number>>}
50
   * @private
51
   */
52
  this.coeffs_ = [[0, 0, 0, Number.NaN]];
53
};
54

55

56
/** @override */
57
goog.math.interpolator.Spline1.prototype.setData = function(x, y) {
58
  goog.asserts.assert(
59
      x.length == y.length,
60
      'input arrays to setData should have the same length');
61
  if (x.length > 0) {
62
    this.coeffs_ = this.computeSplineCoeffs_(x, y);
63
    this.x_ = x.slice();
64
  } else {
65
    this.coeffs_ = [[0, 0, 0, Number.NaN]];
66
    this.x_ = [];
67
  }
68
};
69

70

71
/** @override */
72
goog.math.interpolator.Spline1.prototype.interpolate = function(x) {
73
  var pos = goog.array.binarySearch(this.x_, x);
74
  if (pos < 0) {
75
    pos = -pos - 2;
76
  }
77
  pos = goog.math.clamp(pos, 0, this.coeffs_.length - 1);
78

79
  var d = x - this.x_[pos];
80
  var d2 = d * d;
81
  var d3 = d2 * d;
82
  var coeffs = this.coeffs_[pos];
83
  return coeffs[0] * d3 + coeffs[1] * d2 + coeffs[2] * d + coeffs[3];
84
};
85

86

87
/**
88
 * Solve for the spline coefficients such that the spline precisely interpolates
89
 * the data points.
90
 * @param {Array<number>} x The abscissa of the spline data points.
91
 * @param {Array<number>} y The ordinate of the spline data points.
92
 * @return {!Array<!Array<number>>} The spline interval coefficients.
93
 * @private
94
 */
95
goog.math.interpolator.Spline1.prototype.computeSplineCoeffs_ = function(x, y) {
96
  var nIntervals = x.length - 1;
97
  var dx = new Array(nIntervals);
98
  var delta = new Array(nIntervals);
99
  for (var i = 0; i < nIntervals; ++i) {
100
    dx[i] = x[i + 1] - x[i];
101
    delta[i] = (y[i + 1] - y[i]) / dx[i];
102
  }
103

104
  // Compute the spline coefficients from the 1st order derivatives.
105
  var coeffs = [];
106
  if (nIntervals == 0) {
107
    // Nearest neighbor interpolation.
108
    coeffs[0] = [0, 0, 0, y[0]];
109
  } else if (nIntervals == 1) {
110
    // Straight line interpolation.
111
    coeffs[0] = [0, 0, delta[0], y[0]];
112
  } else if (nIntervals == 2) {
113
    // Parabola interpolation.
114
    var c3 = 0;
115
    var c2 = (delta[1] - delta[0]) / (dx[0] + dx[1]);
116
    var c1 = delta[0] - c2 * dx[0];
117
    var c0 = y[0];
118
    coeffs[0] = [c3, c2, c1, c0];
119
  } else {
120
    // General Spline interpolation. Compute the 1st order derivatives from
121
    // the Spline equations.
122
    var deriv = this.computeDerivatives(dx, delta);
123
    for (var i = 0; i < nIntervals; ++i) {
124
      var c3 = (deriv[i] - 2 * delta[i] + deriv[i + 1]) / (dx[i] * dx[i]);
125
      var c2 = (3 * delta[i] - 2 * deriv[i] - deriv[i + 1]) / dx[i];
126
      var c1 = deriv[i];
127
      var c0 = y[i];
128
      coeffs[i] = [c3, c2, c1, c0];
129
    }
130
  }
131
  return coeffs;
132
};
133

134

135
/**
136
 * Computes the derivative at each point of the spline such that
137
 * the curve is C2. It uses not-a-knot boundary conditions.
138
 * @param {Array<number>} dx The spacing between consecutive data points.
139
 * @param {Array<number>} slope The slopes between consecutive data points.
140
 * @return {!Array<number>} The Spline derivative at each data point.
141
 * @protected
142
 */
143
goog.math.interpolator.Spline1.prototype.computeDerivatives = function(
144
    dx, slope) {
145
  var nIntervals = dx.length;
146

147
  // Compute the main diagonal of the system of equations.
148
  var mainDiag = new Array(nIntervals + 1);
149
  mainDiag[0] = dx[1];
150
  for (var i = 1; i < nIntervals; ++i) {
151
    mainDiag[i] = 2 * (dx[i] + dx[i - 1]);
152
  }
153
  mainDiag[nIntervals] = dx[nIntervals - 2];
154

155
  // Compute the sub diagonal of the system of equations.
156
  var subDiag = new Array(nIntervals);
157
  for (var i = 0; i < nIntervals; ++i) {
158
    subDiag[i] = dx[i + 1];
159
  }
160
  subDiag[nIntervals - 1] = dx[nIntervals - 2] + dx[nIntervals - 1];
161

162
  // Compute the super diagonal of the system of equations.
163
  var supDiag = new Array(nIntervals);
164
  supDiag[0] = dx[0] + dx[1];
165
  for (var i = 1; i < nIntervals; ++i) {
166
    supDiag[i] = dx[i - 1];
167
  }
168

169
  // Compute the right vector of the system of equations.
170
  var vecRight = new Array(nIntervals + 1);
171
  vecRight[0] =
172
      ((dx[0] + 2 * supDiag[0]) * dx[1] * slope[0] + dx[0] * dx[0] * slope[1]) /
173
      supDiag[0];
174
  for (var i = 1; i < nIntervals; ++i) {
175
    vecRight[i] = 3 * (dx[i] * slope[i - 1] + dx[i - 1] * slope[i]);
176
  }
177
  vecRight[nIntervals] =
178
      (dx[nIntervals - 1] * dx[nIntervals - 1] * slope[nIntervals - 2] +
179
       (2 * subDiag[nIntervals - 1] + dx[nIntervals - 1]) * dx[nIntervals - 2] *
180
           slope[nIntervals - 1]) /
181
      subDiag[nIntervals - 1];
182

183
  // Solve the system of equations.
184
  var deriv = goog.math.tdma.solve(subDiag, mainDiag, supDiag, vecRight);
185

186
  return deriv;
187
};
188

189

190
/**
191
 * Note that the inverse of a cubic spline is not a cubic spline in general.
192
 * As a result the inverse implementation is only approximate. In
193
 * particular, it only guarantees the exact inverse at the original input data
194
 * points passed to setData.
195
 * @override
196
 */
197
goog.math.interpolator.Spline1.prototype.getInverse = function() {
198
  var interpolator = new goog.math.interpolator.Spline1();
199
  var y = [];
200
  for (var i = 0; i < this.x_.length; i++) {
201
    y[i] = this.interpolate(this.x_[i]);
202
  }
203
  interpolator.setData(y, this.x_);
204
  return interpolator;
205
};
206

207