function [ x, y, xs, ys ] = bsplineRunge(num_points, degree, omit_knots)
% bsplineRunge: Generates some bspline interpolants of Runge's function.
%
%   bsplineRunge(num_points, degree, omit_knots)
%
% This function constructs some spline interpolants (continuous derivatives
% up to one less than their degree) of Runge's famous example (taken from
% Heath section 7.3):
%
%        f(x) = 1 / (1 + 25x^2)
%
% A bspline basis (equivalent to piecewise polynomials) is used.  The
% results are plotted.  Also plotted is the error in the interpolant.
%
% Input Parameters:
%
%   num_points: Odd positive integer.  Number of sample points on Runge's
%   function.  Points are chosen equally spaced in the unit interval.
%   Default = 9.
%
%   degree: Positive integer.  Degree of the piecewise polynomials.  Must be
%   less than or equal to the number of points, but is typically much
%   smaller.  Default = 3 (cubic spline).
%
%   omit_knots: Boolean.  For degree ~= 1, the number of functions in the
%   basis set with a knot at each abscissae is not equal to the number of
%   abscissae, and hence the linear system we need to solve for the
%   coefficients will not be square.  There are two ways of getting a
%   correct number of basis functions.  One is to omit some abscissae from
%   the knot set (or if degree == 0, add an abscissa).  The alternative is
%   to omit some basis functions (does not work for degree == 0).  The
%   latter is appealing, except that it is not balanced when degree is
%   even, and does not work for degree == 0.  Default = 1.
%
% Output Parameters:
%
%   x: The abscissae.
%
%   y: Runge's function evaluated at the abscissae.
%
%   xs: The sample points for the plot.
%
%   ys: The interpolant's values at the sample points xs.
%
% Ian Mitchell for CPSC 303 2/9/08.

  % Assign default inputs.
  if(nargin < 1)
    num_points = 9;
  elseif(mod(num_points, 2) == 0)
    error('num_points argument must be odd.');
  end
  
  if(nargin < 2)
    degree = 3;
  end

  if(nargin < 3)
    omit_knots = 1;
  end
  
  % What is our function?  Use element-wise divide and multiply because we
  % want to compute f for multiple x points all in one function call.
  f = inline('1 ./ (1 + 25 * x.^2)', 'x');

  % Over what range will be choose abscissae?
  x_range = [ -1, +1 ];
  
  % Over what range do we want to display the results.  Choose something
  % a little larger than the interpolation range, and a little larger
  % than the function's range.
  plot_range_x = x_range + [ -0.25, +0.25 ];
  % We should really evaluate f to figure out this value, but I'll
  % hardcode it for simplicity.
  plot_range_y = [ -0.1, +2.1 ];
  
  % Get a new figure window.
  figure;
  
  % Interpolation abscissae.
  % Apply transpose to get column vector.
  x = linspace(x_range(1), x_range(2), num_points)';
    
  % Interpolation data values.
  y = feval(f, x);

  if omit_knots
    % Knots for the b-splines are the sample points; however, we need to repeat
    % the first and last knot an extra degree times (for a total of
    % degree+1).  For degrees greater than 1, we also need to omit some
    % points from the knot set in order to get the right number of bases.
    % For even degrees, we will omit the middle abscissa.  For degrees three
    % or more, we will omit the necessary number of abscissae from each end
    % starting at the second and second to last.  For degree zero, we need
    % to add a knot, which we do by choosing all the points in between
    % abscissae (num_points - 1) as well as two points just beyond the first
    % and last abscissae.  Finally, we need to move the last (duplicated)
    % knots outward slightly so that there is a basis with nonzero value at
    % the last abscissa.
    mid_point = 0.5 * (num_points - 1) + 1;
    switch(degree)
     case 0
      t = [ 1.5 * x(1) - 0.5 * x(2); ...
            0.5 * (x(1:end-1) + x(2:end));
            1.5 * x(end) - 0.5 * x(end-1) ];
     case 1
      t = [ x(1); x; x(end) + eps ];
     case 2
      t = [ x(1) * ones(degree,1); x(1:mid_point-1); ...
            x(mid_point+1:end); (x(end) + eps) * ones(degree, 1) ];
     case 3
      t = [ x(1) * ones(degree+1,1); x(3:end-2); ...
            (x(end) + eps) * ones(degree+1, 1) ];
     case 4
      t = [ x(1) * ones(degree+1,1); x(3:mid_point-1); ...
            x(mid_point+1:end-2); (x(end) + eps) * ones(degree+1, 1) ];
     otherwise
      error('Knot set has not been defined for degree %d', degree);
    end

    % We have the right number of knots to form a square linear system,
    % so we need not omit any bases.
    omit_left = 0;
    omit_right = 0;
    
  else

    if(degree == 0)
      error('Zero degree b-spline basis needs extra knots.');
    end
    
    % Create dummy knots outside the abscissae.
    t = [ 2 * x(1) - flipud(x(2:degree+1)); x; ...
          2 * x(end) - flipud(x(end-degree:end-1)) ];

    % We may need to omit some bases to get a square linear system.  Omit
    % bases from the left and right, trying to keep the omission
    % balanced.  Unfortunately, for even degree, we need to omit one more
    % (we choose to omit the extra one on the right).
    omit_left = floor(0.5 * degree);
    omit_right = degree - omit_left - 1;
  end
    
  % How many bspline bases are there?
  num_bases = length(t) - 1 - degree - 1;
  
  % Basis function matrix.  
  A = zeros(length(x), length(x));
  for j = omit_left : num_bases - omit_right
    A(:,j+1) = bsplineBasis(x, t, j, degree);
  end
    
  % Find coefficients.
  c = A \ y;
    
  % Range of interest for the plot.
  xs = linspace(plot_range_x(1), plot_range_x(2), 201);

  % Evaluate the function we are trying to match.
  yf = feval(f, xs);

  % Evaluate the interpolant.
  ys = zeros(size(xs));
  for j = omit_left : num_bases - omit_right
    ys = ys + c(j+1) * bsplineBasis(xs, t, j, degree);
  end
  
  % Clear the figure window.
  clf;

  subplot(2,1,1);
  % Plot the function we are trying to match.
  hf = plot(xs, yf, 'b-');
  hold on;
    
  % Also plot the interpolant between interpolation points, to compare
  % results.
  hi = plot(xs, ys, 'r.');
    
  % Plot the interpolation points.
  hp = plot(x, y, 'r*');

  % Label the plot.
  ylabel('f(x)');
  title('B-spline example - Ian Mitchell for CPSC 303');
  legend('Runge function', [ 'degree ' num2str(degree) ' spline' ]);
    
  % Easier to see lines.
  set([ hf; hi ], 'LineWidth', 2);
  set(hp, 'MarkerSize', 9);

  % Interpolant can get really crazy, so we'll tighten the plot bounds.
  axis([ plot_range_x, plot_range_y ]);
    
  % Plot the error.  We will only look at magnitude of error, and we will use
  % a separate subplot so that we can examine the error with logarithmic
  % scaling.  Note how error gets much smaller near the interpolation
  % points (it goes to exactly zero at the interpolation points).
  subplot(2,1,2);
  plot_error = abs(yf - ys);
  he = semilogy(xs, plot_error);

  % Matlab seems to use weird axis bounds for this plot, so we will
  % tighten them up.  First, get the current bounds so as to keep the
  % range in the y-axis.
  axis_error = axis;
  % Now adjust the range in the x-axis to agree with the range for the
  % plot above.
  axis_error(1:2) = plot_range_x;
  % And finally, reset the axis.
  axis(axis_error);
    
  % Label the plot.  Both axes share the same label in the x-axis.
  % Only one line in this plot, so y-axis label substitutes for a legend.
  ylabel('error');
  xlabel('x');

  set(he, 'LineWidth', 2);