%Demo for 'derivativeNE.m'
%After 5 clicks on 'run' you get the results for sine curves
%No plot for dy1c(x1c), since dy1c is nearly on top of dy1
clc
close all
%Create parbola or sine
if exist('curveTypeCounter','var')
  curveTypeCounter=curveTypeCounter+1;
  if curveTypeCounter>10
    curveTypeCounter=1;
  end
else
  clear
  curveTypeCounter=1;
end
curveType=(curveTypeCounter>5);
%  curveType=1;%0=force parabola, 1=force sine
%Create parbola or sine
n=round(rand(1)*150+10);%Number of samples, min=10, max=160
% n=20;%Force n samples
if curveType==0
  %Parabola with spike at x=0
  xr=rand(n,1)*2;
  x=sort(xr)-1;
  %create spike
  indxL=find(x<-0.1,1,'last');
  indxR=find(x>0.15,1,'first');
  x=[x(1:indxL);-0.1;0;0.15;x(indxR:end)];
  x(1)=-1;
  x(end)=1;
  y=x.^2;
  y(indxL+2)=0.05;%Spike amplitude
  legendText1=sprintf('y=x^{2} with spike at x=0');
else
  %Sine
  xr=rand(n,1)*2*pi;
  x=sort(xr);
  x(1)=0;
  x(end)=2*pi;
  y=sin(x);
  legendText1='y=sin(x)';
end

[dy1,dy2,x1c,dy1c]=derivativeNE(x,y);%Calculate derivatives

%Plot results
nl=char(10);%new line, compatible to older Matlab releases
hfig=figure(1);
clf
hfig.Position(3)=hfig.Position(4)*1.8;
yyaxis left
p1=plot(x,y,'.-g');% original curve
xlabel('x')
if curveType==0
  %Parabola
  ylabel('y = x^2 and spike at x=0')
  xticks([-1,-.5,-.1,0,.15,.5,1])
  text(1.27,0.5,['The exact derivates' nl 'at x=-0.1, 0, 0.15' nl 'are not defined as numbers.' nl ...
    'They are for the 1st derivative' nl 'Heaviside functions' nl ...
    'and for the 2nd derivative' nl 'Dirac functions.'],'BackgroundColor',[1,.85,.85],'FontSize',8)
  yyaxis right
  xexact=-1:0.01:1;
  dy1exact=2*xexact;
  p2=plot(xexact(1:91),dy1exact(1:91),'color',[.7,.7,1],'LineWidth',4);%left part exact 1st derivative
  hold on
  indxx0=find(x==0);
  dyL=(y(indxx0)-y(indxx0-1))/(x(indxx0)-x(indxx0-1));
  dyR=(y(indxx0+1)-y(indxx0))/(x(indxx0+1)-x(indxx0));
  plot(xexact([91,91,101,101]),[dy1exact(91),dyL,dyL,dyR],'-','color',[.7,.7,1],'LineWidth',4);%left part of spike exact 1st derivative
  plot(xexact([101,116,116]),[dyR,dyR,dy1exact(116)],'-','color',[.7,.7,1],'LineWidth',4);%right part of spike exact 1st derivative
  plot(xexact(116:end),dy1exact(116:end),'-','color',[.7,.7,1],'LineWidth',4);%right part exact 1st derivative
  p3=plot(xexact([1,90]),[2,2],'-','color',[1,.7,.8],'LineWidth',4);%left part exact 2nd derivative
  plot(xexact([92,115]),[0,0],'-','color',[1,.7,.8],'LineWidth',4);%left and right part of spike exact 2nd derivative
  plot(xexact([117,end]),[2,2],'-','color',[1,.7,.8],'LineWidth',4);%right part exact 2nd derivative
  %Plot Diracs
  quiver(xexact(91),0,0,8,'-','Color',[1,0.7,0.7],'LineWidth',2.5,'MaxHeadSize',.06)
  quiver(xexact(101),0,0,-6,'-','Color',[1,0.7,0.7],'LineWidth',2.5,'MaxHeadSize',.08)
  quiver(xexact(116),0,0,8,'-','Color',[1,0.7,0.7],'LineWidth',2.5,'MaxHeadSize',.06)
  
  p4=plot(x,dy1,'b-');%1st numerical derivative interpolated
  % p4=plot(x1s,dy1c,'.y');%1st derivative centered
  p5=plot(x,dy2,'+r-');%2nd numerical derivative
  legend([p1,p4,p2,p5,p3],{sprintf('y=x^{2} with spike at x=0'),'1st derivative numerical','1st derivative exact',...
    '2nd derivative numerical','2nd derivative exact'},'Location','northeastoutside')
  ylim([-8,12])
  
  %Calculate Mean Absolute Errors and max errors without range of spike
  %1st derivative at x
  x4MAEdy1=[x(1:indxx0-2);x(indxx0+2:end)];%     x for MAE
  dy14MAE= [dy1(1:indxx0-2);dy1(indxx0+2:end)];% dy1 for MAE
  dy1AD=abs(dy14MAE-2*x4MAEdy1);%                Absolute differences
  [dy1Max,dy1Maxindx]=max(dy1AD);%               Max absolute difference of dy1-dy1exact
  dy1MAE= sum(dy1AD)/numel(x4MAEdy1);%           MAE dy1
  
  %1st derivative at x1c including spike
  x4MAEdy1c=[x1c(1:indxx0-2);x1c(indxx0+1:end)];%    x for MAE
  dy1c4MAE= [dy1c(1:indxx0-2);dy1c(indxx0+1:end)];%  dy1c for MAE
  dy1cAD=abs(dy1c4MAE-2*x4MAEdy1c);%                 Absolute differences
  [dy1cMax,dy1cMaxindx]=max(dy1cAD);%                Max absolute difference of dy1c-dy1exact
  dy1cMAE= sum(dy1cAD)/numel(x4MAEdy1c);%            MAE dy1c
  
  %2nd derivative
  x4MAEdy2=[x(2:indxx0-3);x(indxx0+3:end-1)];%      x for dy2 MAE
  dy24MAE=[dy2(2:indxx0-3);dy2(indxx0+3:end-1)];%   dy2 for MAE
  dy2AD=abs(dy24MAE-2);%                            Absolute differences
  [dy2Max,dy2Maxindx]=max(dy2AD);%                  Max absolute difference of dy2-dy2exact
  dy2MAE=sum(dy2AD)/numel(x4MAEdy2);%               MAE dy2
  
  %Output text for MAE
  text(1.27,-6,['Mean Abs. Error     Max abs. error' nl ...
     'dy1:  ' num2str(dy1MAE,'%.1e') '         ' num2str(dy1Max,2) ' @x=' num2str(x4MAEdy1(dy1Maxindx),2) nl ...
     'dy1c: ' num2str(dy1cMAE,'%.1e') '         ' num2str(dy1cMax,2) ' @x=' num2str(x4MAEdy1c(dy1cMaxindx),2) nl ...
    'dy2:  ' num2str(dy2MAE,'%.1e') '         ' num2str(dy2Max,2) ' @x=' num2str(x4MAEdy2(dy2Maxindx),2) nl ...
    'Range around spike is excluded'], ...
    'BackgroundColor',[0.98,.98,.98],'FontSize',8)

else
  %Sine
  ylim([-1.05,1.05])
  ylabel('y = sin(x)')
  xticks(0:pi/2:2*pi+eps)
  set(gca,'TickLabelInterpreter','latex','XTickLabel',{'0','$\frac{\pi}{2}$','$\pi$','$\frac{3\pi}{2}$','$2\pi$'})
%   set(gca,'XMinorGrid','on')%Only for test
  yyaxis right
  xexact=linspace(0,2*pi,200);
  p2=plot(xexact,cos(xexact),'color',[.8,.8,1],'LineWidth',4);%         exact 1st derivative
  hold on
  p3=plot(xexact,-sin(xexact),'-','color',[1,.8,0.8],'LineWidth',4);%   exact 2nd derivative
  p4=plot(x,dy1,'b-');%                                                 1st numerical derivative interpolated
  % p4=plot(x1s,dy1c,'y-');%1st derivative centered
  p5=plot(x,dy2,'+r-');%                                                2nd numerical derivative
  legend([p1,p4,p2,p5,p3],{'y = sin(x)','1st derivative numerical','1st derivative exact',...
    '2nd derivative numerical','2nd derivative exact'},'Location','northeastoutside')
  ylim([-1.05,1.05])
  xlim([0,2*pi])
  
  %Calculate Mean Absolute Errors and max errors 
  %1st derivative at x
  dy1AD=abs(dy1-cos(x));%            Absolute differences
  [dy1Max,dy1Maxindx]=max(dy1AD);%   Max absolute difference of dy1-dy1exact
  dy1MAE= sum(dy1AD)/numel(x);%      MAE dy1
  
  %1st derivative at x1c 
  dy1cAD=abs(dy1c-cos(x1c));%             Absolute differences
  [dy1cMax,dy1cMaxindx]=max(dy1cAD);%     Max absolute difference of dy1c-dy1exact
  dy1cMAE= sum(dy1cAD)/numel(x1c);%       MAE dy1c
  
  %2n derivative
  x4MAEdy2=x(2:end-1);%                x for dy2 MAE
  dy24MAE=dy2(2:end-1);%               dy2 for MAE
  dy2AD=abs(dy24MAE+sin(x4MAEdy2));%   Absolute differences
  [dy2Max,dy2Maxindx]=max(dy2AD);%     Max absolute difference of dy2-dy2exact
  dy2MAE=sum(dy2AD)/numel(x4MAEdy2);%  MAE dy2
  
  
  %Output text for MAE
    text(7.2,0,['Mean Abs. Error     Max abs. error' nl ...
     'dy1:  ' num2str(dy1MAE,'%.1e') '         ' num2str(dy1Max,'%.1e') ' @x=' num2str(x(dy1Maxindx),2) nl ...
     'dy1c: ' num2str(dy1cMAE,'%.1e') '         ' num2str(dy1cMax,'%.1e') ' @x=' num2str(x1c(dy1cMaxindx),2) nl ...
    'dy2:  ' num2str(dy2MAE,'%.1e') '         ' num2str(dy2Max,'%.1e') ' @x=' num2str(x4MAEdy2(dy2Maxindx),2)], ...
    'BackgroundColor',[0.98,.98,.98],'FontSize',8)

end
ylabel('1st and 2nd derivatives')
hold off
grid on

function [dy1,dy2,x1c,dy1c]=derivativeNE(x,y)
%Calculate the 1st and 2nd derivative for Not Equal spaced samples.
%
%INPUT
%  x as column or row vector, more than 2 values
%  y as column or row vector, same amount of values as for x
%
%OUTPUT
%  dy1: 1st derivative, y'(x)
%       
%  dy2: 2nd derivative, y"(x)
%       The first dy2 at x(1) and the last dy2 at x(end) don't exist and are set to NaN
%  x1c and dy1c return a more accurate (about 3 times better) 1st derivative dy1c(x1c), but with new x-values
%
%  If x or y is a row vector, all outputs are returned as row vectors.
%  On error NaN is returned
%
%Syntax examples:
%  dy1=derivativeNE([1;3;4;5],[1;9;16;25]);           % 1st derivative y'(x)
%  [~,dy2]=derivativeNE([1;3;4;5],[1;9;16;25]);       % 2nd derivative y"(x)
%  [~,~,x1c,y1c]=derivativeNE([1;3;4;5],[1;9;16;25]); % 1st derivative more accurate (about 3 times), y'(x1c)
%
%Peter Seibold, August 2023

%Check inputs
if size(x,2)>1 || size(y,2)>1
  %Either x or y is a row vector
  isRowVector=true;
  x=x(:);%convert to column vector
  y=y(:);
else
  isRowVector=false;
end
if numel(x)~=numel(y) || numel(x)<3
  dy1=NaN; dy2=NaN; x1c=NaN; dy1c=NaN;%Return all output as NaN
  disp('x and y must have the same size and at least 3 elements!')
  return
end

%1st derivative
dx=diff(x);
dy=diff(y);
dy1c=dy./dx;%          1st derivative for centered x
x1c=x(1:end-1)+dx/2;%  center each x (shift each x by half distance to next x)

%interpolate/extrapolate dy1c(x1c) to match x, result: dy1(x)
y2=[dy1c(2);dy1c(2:end);dy1c(end)];   %repeat border values for extrapolation
y1=[dy1c(1);dy1c(1:end-1);dy1c(end-1)];
x1=[x1c(1);x1c(1:end-1);x1c(end-1)];
x2=[x1c(2);x1c(2:end);x1c(end)];
dy1=(y2-y1)./(x2-x1).*(x-x2)+y2;      %interpolate/extrapolate

%2nd derivative
y32x32=dy1c(2:end);                   %1st derivatives at even positions
y21x21=dy1c(1:end-1);                 %1st derivatives at odd positions
dx31=y32x32;                          %create vector with dummy values, with two elements less than number of x-elements
dx31(1:2:end)=x(3:2:end)-x(1:2:end-2);%x-difference odd pairs
dx31(2:2:end)=x(4:2:end)-x(2:2:end-2);%x-difference even pairs
dy2c=2*(y32x32-y21x21)./dx31;         %2nd derivatives dy2 at center of the two surrounding samples
x2c=dy2c;                             %create vector with dummy values, with two elements less than number of x-elements
x2c(1:2:end)=(x(1:2:end-2)+x(3:2:end))/2;%x for dyc at center of the two surrounding samples, (overwrite dummy values)
x2c(2:2:end)=(x(2:2:end-2)+x(4:2:end))/2;%dito for even x centered, now we have all center x2c for 2nd derivatives

%interpolate dy2c(x2c) to match x, result: dy2(x)
if numel(x)>3
  dy2=[NaN;interp1(x2c,dy2c,x(2:end-1),'linear','extrap');NaN];%interpolate and pad unavailable border values
else
  %Interpolation needs at least two dy2c samples
  dy2=[NaN;dy2c;NaN];%    output only the single dy2(x)
end

%Convert back to user format
if isRowVector
  %Either x or y input was a row vector
  %Convert to row vector
  dy1=dy1';
  dy2=dy2';
  x1c=x1c';
  dy1c=dy1c';
end