function varargout = Numderiv(varargin)
% NUMDERIV Application M-file for Numderiv.fig
%    FIG = NUMDERIV launch Numderiv GUI.
%    NUMDERIV('callback_name', ...) invoke the named callback.

% Last Modified by GUIDE v2.0 03-Mar-2002 18:59:11

if nargin == 0  % LAUNCH GUI

	fig = openfig(mfilename,'reuse');
    hold on;

	% Use system color scheme for figure:
	set(fig,'Color',get(0,'defaultUicontrolBackgroundColor'));

	% Generate a structure of handles to pass to callbacks, and store it. 
	handles = guihandles(fig);
	guidata(fig, handles);

	if nargout > 0
		varargout{1} = fig;
	end

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

	try
		if (nargout)
			[varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
		else
			feval(varargin{:}); % FEVAL switchyard
		end
	catch
		disp(lasterr);
	end

end
%| ABOUT CALLBACKS:
%| GUIDE automatically appends subfunction prototypes to this file, and 
%| sets objects' callback properties to call them through the FEVAL 
%| switchyard above. This comment describes that mechanism.
%|
%| Each callback subfunction declaration has the following form:
%| <SUBFUNCTION_NAME>(H, EVENTDATA, HANDLES, VARARGIN)
%|
%| The subfunction name is composed using the object's Tag and the 
%| callback type separated by '_', e.g. 'slider2_Callback',
%| 'figure1_CloseRequestFcn', 'axis1_ButtondownFcn'.
%|
%| H is the callback object's handle (obtained using GCBO).
%|
%| EVENTDATA is empty, but reserved for future use.
%|
%| HANDLES is a structure containing handles of components in GUI using
%| tags as fieldnames, e.g. handles.figure1, handles.slider2. This
%| structure is created at GUI startup using GUIHANDLES and stored in
%| the figure's application data using GUIDATA. A copy of the structure
%| is passed to each callback.  You can store additional information in
%| this structure at GUI startup, and you can change the structure
%| during callbacks.  Call guidata(h, handles) after changing your
%| copy to replace the stored original so that subsequent callbacks see
%| the updates. Type "help guihandles" and "help guidata" for more
%| information.
%|
%| VARARGIN contains any extra arguments you have passed to the
%| callback. Specify the extra arguments by editing the callback
%| property in the inspector. By default, GUIDE sets the property to:
%| <MFILENAME>('<SUBFUNCTION_NAME>', gcbo, [], guidata(gcbo))
%| Add any extra arguments after the last argument, before the final
%| closing parenthesis.

% --------------------------------------------------------------------
function varargout = leftendpt_Callback(h, eventdata, handles, varargin)
try
    disp(get(h,'Value'));
    set(h,'Value',eval(get(h,'String')));
catch
    set(h,'String','0.0');
    set(h,'Value',0.0);
end

% --------------------------------------------------------------------
function varargout = rightendpt_Callback(h, eventdata, handles, varargin)
try
    set(h,'Value',eval(get(h,'String')));
catch
    set(h,'String','1.0');
    set(h,'Value',1.0);
end

% --------------------------------------------------------------------
function varargout = dfnfx_Callback(h, eventdata, handles, varargin)
try
    testfcn = inline(get(h,'String'));
    testfcn(eval(get(handles.leftendpt,'String')));
catch
    set(h,'String','sin(pi*x)');
end

% --------------------------------------------------------------------
function varargout = nodes_Callback(h, eventdata, handles, varargin)

try
    set(h,'Value',max(4,round(eval(get(h,'String')))));
    set(h,'String',num2str(get(h,'Value')));
catch
    set(h,'String','4');
    set(h,'Value',4);
end

% --------------------------------------------------------------------
function varargout = noiselevel_Callback(h, eventdata, handles, varargin)
try
    set(h,'Value',max(0,min(100,eval(get(h,'String')))));
    set(h,'String',num2str(get(h,'Value')));
catch
    set(h,'String','0.0');
    set(h,'Value',0.0);
end

% ----------------------------------------------------------
function varargout = cleargraphbutton_Callback(h, eventdata, handles, varargin)
set(handles.noderr,'String','Max Node Error: 0');
set(handles.alpha,'String','');

cla;


% --------------------------------------------------------------------
function varargout = drawgraphbutton_Callback(h, eventdata, handles, varargin)
leftendpoint = get(handles.leftendpt,'Value');
rightendpoint = get(handles.rightendpt,'Value');
nodes = get(handles.nodes,'Value');
noiselevel = get(handles.noiselevel,'Value');
if (leftendpoint >= rightendpoint)
    disp('Endpoints are reversed. Fix this.');
    return;
end
xnodes = linspace(leftendpoint,rightendpoint,nodes);
ynodes = xnodes;
rfcn = inline('40*(x+1/4)*(x-1/3)*(x-2/3)');
ufcn = inline('exp(x*(x*(x*(x*10-10)-5/9)+20/9))');
fcn = inline(get(handles.dfnfx,'String'));
for ii = (1:nodes)
   ynodes(ii) = fcn(xnodes(ii));
end
rand('state',nodes); % initial seed depends only on node size
% now the problem of how to specify error...1st is % of val, 2nd % of max
%rndynodes = ynodes.*(1+(noiselevel/100)*(2*rand(size(ynodes))-1));
rndynodes = ynodes+norm(ynodes,inf)*noiselevel/100*(2*rand(size(ynodes))-1);
mnodes = 10*(nodes - 1)+1; % densify existing nodes by factor of 10
mxnodes = linspace(leftendpoint,rightendpoint,mnodes);
mynodes = mxnodes;
for ii = (1:mnodes)
    mynodes(ii) = fcn(mxnodes(ii)); % calculate fcn value
end
% calculate derivative of fcn value reasonable accurately
mynodesp = [-3*mynodes(1)+4*mynodes(2)-mynodes(3),...
            (mynodes(3:mnodes)-mynodes(1:mnodes-2)),...
            mynodes(mnodes-2)-4*mynodes(mnodes-1)+3*mynodes(mnodes)]...
        /(2*(mxnodes(2)-mxnodes(1)));
set(handles.alpha,'String','');
switch get(handles.graphlistbox,'Value')
case 1 % graph function
    plot(mxnodes,mynodes,'r-');
    set(handles.noderr,'String','Max Node Error: 0');
    set(handles.alpha,'String','');
case 2 % graph function with noise
    plot(xnodes,rndynodes,'b-');
    set(handles.noderr,'String',strcat('Max Node Error: ',...
        num2str(norm(rndynodes-mynodes(1:10:mnodes),inf))));
    set(handles.alpha,'String','');
case 3 % graph natural cubic spline through function with noise
    pp = NCSpp(xnodes,rndynodes);
    pynodes = PPeval(pp,mxnodes);
    plot(mxnodes,pynodes,'m-');
    set(handles.noderr,'String',strcat('Max Node Error: ',...
        num2str(norm(pynodes -mynodes,inf))));
    set(handles.alpha,'String','');
case 4 % graph function derivative
    plot(mxnodes,mynodesp,'r-');
        set(handles.noderr,'String','Max Node Error: 0');
        set(handles.alpha,'String','');
case 5 % graph forward difference approximate derivative with noise
    tmpnodes = [(rndynodes(2:nodes)-rndynodes(1:nodes-1)),...
        rndynodes(nodes)-rndynodes(nodes-1)]/(xnodes(2)-xnodes(1));
    plot(xnodes,tmpnodes,'b-');
    set(handles.noderr,'String',strcat('Max Node Error: ',...
        num2str(norm(tmpnodes -mynodesp(1:10:mnodes),inf))));
    set(handles.alpha,'String','');
case 6  % graph centered difference approximate derivative with noise
    tmpnodes = [-3*rndynodes(1)+4*rndynodes(2)-rndynodes(3),...
            (rndynodes(3:nodes)-rndynodes(1:nodes-2)),...
            rndynodes(nodes-2)-4*rndynodes(nodes-1)+3*rndynodes(nodes)]...
            /(2*(xnodes(2)-xnodes(1)));
    plot(xnodes,tmpnodes,'g-');
    set(handles.noderr,'String',strcat('Max Node Error: ',...
        num2str(norm(tmpnodes - mynodesp(1:10:mnodes),inf))));
    set(handles.alpha,'String','');
case 7  % graph cubic spline approximate derivative with noise
    pp = NCSpp(xnodes,rndynodes);
    ppp = PPder(pp);
    mynodes = PPeval(ppp,mxnodes);
    plot(mxnodes,mynodes,'m-');
    set(handles.noderr,'String',strcat('Max Node Error: ',...
        num2str(norm(mynodes -mynodesp,inf))));
    set(handles.alpha,'String','');
case 8 % regularized cubic spline approximate derivative with noise
    alp = inputdlg('Input regularization parameter: ',...
        'Regularized Natural Cubic Spline');
    alp = str2num(alp{1});
    pp = RNCSpp(xnodes,rndynodes,alp);
    ppp = PPder(pp);
    mynodes = PPeval(ppp,mxnodes);
    plot(mxnodes,mynodes,'m-');
    set(handles.noderr,'String',strcat('Max Node Error: ',...
        num2str(norm(mynodes -mynodesp,inf))));
    set(handles.alpha,'String',strcat('alpha= ',...
        num2str(alp)));
case 9 % least squares polynomial fit
    dgr = inputdlg('Input degree of polynomial fit: ',...
        'Polynomial Fitting of Data');
    dgr = max(1,str2num(dgr{1}));
    pp = polyfit(xnodes,rndynodes,dgr);
    ppp = polyder(pp);
    mynodes = polyval(ppp,mxnodes);
    plot(mxnodes,mynodes,'c-');
    set(handles.noderr,'String',strcat('Max Node Error: ',...
        num2str(norm(mynodes -mynodesp,inf))));
    set(handles.alpha,'String',strcat('degree= ',...
        num2str(dgr)));
case 10 % population u(t) with noise
    for ii = (1:nodes)
       ynodes(ii) = ufcn(xnodes(ii));
    end
    rndynodes = ynodes+norm(ynodes,inf)*noiselevel/100*(2*rand(size(ynodes))-1);
    plot(xnodes,rndynodes,'b-');
    set(handles.noderr,'String','Noisy u(t)');
    set(handles.alpha,'String','');
case 11 % estimate rate r(t)=u'/u, approx by reg spline
    alp = inputdlg('Input regularization parameter: ',...
        'Regularized Natural Cubic Spline');
    alp = str2num(alp{1});
    for ii = (1:nodes)
       ynodes(ii) = ufcn(xnodes(ii));
    end
    rndynodes = ynodes+norm(ynodes,inf)*noiselevel/100*(2*rand(size(ynodes))-1);
    pp = RNCSpp(xnodes,rndynodes,alp);
    ppp = PPder(pp);
    mynodes = PPeval(ppp,mxnodes)./PPeval(pp,mxnodes);
    plot(mxnodes,mynodes,'m-');
    set(handles.noderr,'String','(du/dt)/u by regularized NCS for u');
    set(handles.alpha,'String',strcat('alpha= ',...
        num2str(alp)));
case 12  % graph centered difference approximate derivative with noise
    for ii = (1:nodes)
       ynodes(ii) = ufcn(xnodes(ii));
    end
    rndynodes = ynodes+norm(ynodes,inf)*noiselevel/100*(2*rand(size(ynodes))-1);
    tmpnodes = [-3*rndynodes(1)+4*rndynodes(2)-rndynodes(3),...
            (rndynodes(3:nodes)-rndynodes(1:nodes-2)),...
            rndynodes(nodes-2)-4*rndynodes(nodes-1)+3*rndynodes(nodes)]...
            /(2*(xnodes(2)-xnodes(1)));
    tmpnodes = tmpnodes./rndynodes;
    plot(xnodes,tmpnodes,'g-');
    set(handles.noderr,'String','(du/dt)/u by centered diffs for u');
    set(handles.alpha,'String','');
case 13  % least squares poly fit
    dgr = inputdlg('Input degree of polynomial fit: ',...
        'Polynomial Fitting of Data');
    dgr = max(1,str2num(dgr{1}));
    for ii = (1:nodes)
       ynodes(ii) = ufcn(xnodes(ii));
    end
    rndynodes = ynodes+norm(ynodes,inf)*noiselevel/100*(2*rand(size(ynodes))-1);
    pp = polyfit(xnodes,rndynodes,dgr);
    ppp = polyder(pp);
    mynodes = polyval(ppp,mxnodes);
    plot(mxnodes,mynodes,'c-');
    set(handles.noderr,'String','(du/dt)/u by LS poly fit for u');
    set(handles.alpha,'String',strcat('degree= ',...
        num2str(dgr)));
    

otherwise
    ;
end

% --------------------------------------------------------------------
function varargout = graphlistbox_Callback(h, eventdata, handles, varargin)

% --------------------------------------------------------------------
function varargout = figurenumderiv_CloseRequestFcn(h, eventdata, handles, varargin)
closereq;

function pp = PP(order,knum,knots,coefs)
% usage: pp = PP(order,knum,knots,coefs)
% description: given input of at least positive integers order and
% knum, this function constructs a PP object with this order and
% knum knots. Last two arguments are optional, but if entered,
% user must define knots before coefs. If these are not supplied,
% the  the coef matrix is initialized to zeros and knots to 0:knum-1.

% local variables:
% j: loop variable

if (nargin < 2 | nargin >4)
  disp('Incorrect paramater list')
  return;
end;
pp.order = order;
pp.knum = knum;
pp.knots = 0:(knum-1);
pp.coefs = zeros(knum-1,order);
if (nargin >= 3)
  if (sum([1,knum] ~=size(knots)))
    disp('Incorrect knots vector')
    clear pp;
    return;
  end;
  for j = 2:knum
    if (knots(j) <= knots(j-1))
      disp('Knots vector is not strictly increasing')
      clear pp;
      return;
    end;
  end;
  pp.knots = knots;
end;
if (nargin == 4)
  if (sum([knum-1,order] ~=size(coefs)))
    disp('Incorrect coefs matrix')
    clear pp;
    return;
  end;
  pp.coefs = coefs;
end;

%

function y = PPeval(pp,x)
% usage: y = PPeval(pp,x)
% description:  given input pp a PP structure and x a vector of
% absicssas, this routine returns the vector of values of the PP
% structure pp at the coordinates of x.
 
% local variables: 
% j,k,m:  index variables
% n: size of node array
% klim: number of knots
% ord: order of pp
% xtmp: temporary variable

n = length(x);

y = zeros(size(x));
klim = pp.knum;
ord = pp.order;
for j = 1:n % loop over entries of x 
  xtmp = x(j); 
  % find index of interval containing entry
  for k = 1:klim-1
    if (xtmp <  pp.knots(k+1))
      break;
    end;
  end;
  dx = xtmp - pp.knots(k);
  % use Horner's method to evaluate the polynomial
  y(j) = pp.coefs(k,1);
  for m = 2:ord
    y(j) = dx*y(j)+pp.coefs(k,m);
  end;
end;

%

function dd = PPder(pp)
% usage: dd = PPder(pp)
% description: given input pp a PP structure, this function
% returns the PP structure representing the derivative of the p.p.
% function represented by pp.

% local variables: 
% j,k:  index variables
% klim: number of knots
% ord: order of pp

dd.knum = pp.knum;
dd.knots = pp.knots;
dd.order =pp.order -1;
dd.coefs = zeros(dd.knum -1, dd.order);
deg = pp.order - 1;
for j = 1:dd.knum-1
  for k = 1:dd.order
    dd.coefs(j,k) = (deg-k+1)*pp.coefs(j,k);
  end;
end;

%

function pp = HCpp(x,f,fp) 
% usage: pp = HCpp(x,f,fp) 
% description: this function returns a PP structure representing
% the Hermite cubic p.p. with node row vector x, corresponding
% function values in the row vector f and derivative values in the
% row vector fp.


% local variables: 
% n,m: size of node and value arrays
% j: index variable
% h: stepsize

n = length(x);
m= size(f);
m = m(1)+m(2)-1;
if (m ~= n) % error check for size and order
  disp('Size mismatch of abscissa and function ordinate');
  return;
end;
m = length(fp);
if (m ~= n) % error check for size and order
  disp('Size mismatch of abscissa and derivative ordinate');
  return;
end;
pp = PP(4, n, x);
% determine the poly coefs in each knot subinterval
for j = 1:(n-1)
  % set up and solve the system on this subinterval
  h = x(j+1) - x(j);
  pp.coefs(j,:) = ([0, 0, 0, 1; h^3, h^2, h, 1; 0 0 1 0;
3*h^2, 2*h, 1, 0] \ [f(j); f(j+1); fp(j); fp(j+1)] )';  
end;

%

function pp = NKCSpp(xx,ff)
% usage: pp = NKCSpp(x,f)
% description: this function returns the not-a-knot cubic spline
% with nodes at the entries of the row vector of abscissas x and
% corresponding ordinates in the row vector f, excluding the 
% second and penultimate entries of each.

% local variables:
% x,f: vectors with not-a-knot entries removed
% nn,m,n: size of node and value arrays
% j: index variable
% h: stepsize vector
% d: divided difference vector
% Coef: coefficient matrix for the system
% Rhs: right hand side vector
% h0, h1, hnm2, hnm1: extra step sizes
% M: vector of second derivative values
% coefs: coefficient matrix of PP structure to be returned

nn = length(xx);
m= length(ff);
if (m ~= nn | nn<4)  % error check for size and order
  disp('Size mismatch of abscissa and function ordinate or size problem');
  return;
end;
x = [xx(1), xx(3:nn-2), xx(nn)];
f = [ff(1), ff(3:nn-2), ff(nn)];
n = nn - 2; % adjust size to reflect size of x
% construct the stepsize and divided difference vectors
h = x(2:n) - x(1:n-1);
d = (f(2:n) - f(1:n-1))./h;
hh = xx(2:nn) - xx(1:nn-1);
% now set up the system as in Atkinson, p. 169, row by row
Coef = zeros(n,n); 
Rhs = zeros(n,1);
% first row (use equation of exercise 3.34)
Coef(1, 1:2) = [hh(1)+2*hh(2), 2*hh(1)+hh(2)];
Rhs(1) = 6*(ff(1)/hh(1) - h(1)*ff(2)/(hh(1)*hh(2)) + ff(3)/hh(2));
for j = 2:(n-1) % interior rows
  Coef(j, (j-1):(j+1)) = [h(j-1)/6, (h(j-1)+h(j))/3, h(j)/6];
  Rhs(j) = d(j) - d(j-1);
end;
% last row (use equation of exercise 3.34)
Coef(n, (n-1):n) = [hh(nn-2)+2*hh(nn-1), 2*hh(nn-2)+hh(nn-1)]; 
Rhs(n) = 6*(ff(nn-2)/hh(nn-2) - h(n-1)*ff(nn-1)/(hh(nn-2)*hh(nn-1)) + ff(nn)/hh(nn-1));
% solve for M
M = Coef\Rhs;
% finally construct the PP structure to be returned
coefs = zeros(n-1,4);
% given si(x)=ai*(x-xi)^3+bi*(x-xi)^2+ci*(x-xi)+di 
% we deduce from si''(x)=6*ai*(x-xi)+2*bi that
% bi = Mi/2 and ai = (M(i+1)-Mi)/(6*hi)
coefs(:,2) = M(1:n-1)/2;
coefs(:,1) = (M(2:n) - M(1:n-1))./(6*h');
% evaulation gives di = yi and ci = (y(i+1) - ai*hi^3-bi*hi^2-di)/hi
coefs(:,4) = f(1:n-1)';
coefs(:,3) = (f(2:n)' - coefs(:,1).*h'.^3 - coefs(:,2).*h'.^2 - coefs(:,4)) ./ h';
pp = PP(4,n,x,coefs);

%

function pp = NCSpp(x,f)
% usage: pp = NCSpp(x,f)
% description: this function returns the natural cubic spline
% with nodes at the entries of the row vector of abscissas x and
% corresponding ordinates in the row vector f. 

% local variables:
% n,m: size of node and value arrays
% j: index variable
% h: stepsize vector
% d: divided difference vector
% Coef: coefficient matrix for the system
% Rhs: right hand side vector
% M: vector of second derivative values
% coefs: coefficient matrix of PP structure to be returned

n = length(x);
m= length(f);
if (m ~= n) % error check for size and order
  disp('Size mismatch of abscissa and function ordinate');
  return;
end;
% construct the stepsize and divided difference vectors
h = x(2:n) - x(1:n-1);
d = (f(2:n) - f(1:n-1))./h;
% now set up the system as in Atkinson, p. 169, row by row
Coef = zeros(n,n); 
Rhs = zeros(n,1);
% first row
Coef(1, 1) = 1;
Rhs(1) = 0;
for j = 2:(n-1) % interior rows
  Coef(j, (j-1):(j+1)) = [h(j-1)/6, (h(j-1)+h(j))/3, h(j)/6];
  Rhs(j) = d(j) - d(j-1);
end;
% last row
Coef(n, n) = 1; 
Rhs(n) = 0;
% solve for M
M = Coef\Rhs;
% finally construct the PP structure to be returned
coefs = zeros(n-1,4);
% given si(x)=ai*(x-xi)^3+bi*(x-xi)^2+ci*(x-xi)+di 
% we deduce from si''(x)=6*ai*(x-xi)+2*bi that
% bi = Mi/2 and ai = (M(i+1)-Mi)/(6*hi)
coefs(:,2) = M(1:n-1)/2;
coefs(:,1) = (M(2:n) - M(1:n-1))./(6*h');
% evaulation gives di = yi and ci = (y(i+1) - ai*hi^3-bi*hi^2-di)/hi
coefs(:,4) = f(1:n-1)';
coefs(:,3) = (f(2:n)' - coefs(:,1).*h'.^3 - coefs(:,2).*h'.^2 - coefs(:,4)) ./ h';
pp = PP(4,n,x,coefs);

function pp = RNCSpp(x,a,alp)
% usage: pp = RNCSpp(x,a,alp)
% description: this function returns the regularized natural 
% cubic spline with nodes at the entries of the vector
% of abscissas x and corresponding ordinates in the 
% vector a. Here alp is  the regularization parameter.
% Given that the unregularized system is R*c=Q*a, where c is
% the column of moments (2nd deriv values), with 
% both R and Q tridiagonal and diagonally dominant, 
% the regularized system is 
%             (R + alp*Q*Q^{T}*R^{-1})*c = Q*a
% which is equivalent to the interpolating natural cubic 
% spline when alp=0. 

% local variables:
% n,m: size of node and value arrays
% j: index variable
% h: stepsize vector
% R: coefficient matrix for the system
% Q: coefficient matrix for right hand side vector
% c: vector of second derivative values
% coefs: coefficient matrix of PP structure to be returned

% columnize the input data
[m,n]=size(x);
if (n > 1)
    x = x';
end;
[m,n] = size(a);
if (n > 1)
    a = a';
end;
n = length(x);
m= length(a);
if (m ~= n) % error check for size and order
  disp('Size mismatch of abscissa and function ordinate');
  return;
end;
% construct the stepsize and divided difference vectors
h = x(2:n) - x(1:n-1);
% now set up the system as in my notes along lines of de Boor
R = zeros(n,n);
Q = R;
% first row
R(1,1) = 1;
Q(1,1) = 1;
for j = 2:(n-1) % interior rows
  R(j,(j-1):(j+1)) = [h(j-1), (h(j-1)+h(j))*2, h(j)];
  Q(j,(j-1):(j+1)) = [3/h(j-1), -3/h(j-1)-3/h(j), 3/h(j)];
end;
% last row
R(n,n) = 1;
Q(n,n) = 1;
% fix R to be symmetric
R(2,1) = 0;
R(n-1,n) = 0;
% solve for c
c = (R+alp*Q*Q'*inv(R))\(Q*a);
% find new nodes for spline
a = Q\(R*c);
% fix first and last entries, currently spline values at end nodes
c(1) = 0;
c(n) = 0;
% finally construct the PP structure to be returned
coefs = zeros(n-1,4);
% given si(x)=ai+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3 
% we deduce from si''(x)=6*di*(x-xi)+2*ci that
% ai = (ci(+1)-ci)/(3*hi)
coefs(:,2) = c(1:n-1);
coefs(:,1) = (c(2:n) - c(1:n-1))./(3*h);
% evaulation gives ai = yi and bi = (y(i+1) - di*hi^3-ci*hi^2-ai)/hi
coefs(:,4) = a(1:n-1);
coefs(:,3) = (a(2:n) - coefs(:,1).*h .^3 - coefs(:,2).*h .^2 - ...
    coefs(:,4)) ./ h;
pp = PP(4,n,x',coefs);