function demo_standard_projector;
% standard projection (quasiinterpolant) of degree 2
% to a function, subdivision until error bound is reached

% Initialize data
global AMB_Data;
AMB_Data=[];
AMB_Data.npts = 801;  
AMB_Data.Fun = @(x) sin(1./(1.1+sin(4*pi*x/2)));  
AMB_Data.Funstring = '@(x) sin(1./(1.1+sin(4*pi*x/2)))';  
AMB_Data.error = inf;
AMB_Data.errorbound = 1e-5;
AMB_Data.knots = [-1;0;0;1;1;2];
f=AMB_Data.Fun([0;1/2;1]);
AMB_Data.c=([f(1);2*f(2)-(f(1)+f(3))/2;f(3)]);

% Initialize GUI
clf reset;

% Text for title
uicontrol('Units','normalized', ...
  	'BackgroundColor',get(gcf,'Color'), ...
	'HorizontalAlignment','center', ...
        'FontWeight','bold',...
        'FontSize',20,...
	'Position',[.05,.92,.9,.07],...
	'String','Standard Projector',...
	'Style','text');

% Text for function
uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'FontWeight','bold',...
        'HorizontalAlignment','left', ...
	'Position',[.05,.86,.3,.05],...
	'String','Function:',...
	'Style','text');


% Input for function
AMB_Data.uicontrols(1)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'FontWeight','bold',...
        'HorizontalAlignment','left', ...
	'Position',[.2,.86,.4,.05],...
	'String',sprintf('%s',AMB_Data.Funstring),...
	'Style','edit',...
	'Callback',@function_update);

% Text for error bound
uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.05,.81,.2,.05],...
	'String','Error bound:',...
	'Style','text');

% Input for error bound
AMB_Data.uicontrols(2)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.2,.81,.15,.05],...
	'String',sprintf('%8.3g',AMB_Data.errorbound),...
	'Style','edit',...
	'Callback',@errorbound_update);

% Text for error
AMB_Data.uicontrols(3)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.4,.81,.2,.05],...
	'String',sprintf('Error: %8.3g',AMB_Data.error),...
	'Style','text');

% Next button
AMB_Data.uicontrols(5)=uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'String','Next', ...
	'Position',[.85,.81,.1,.05], ...
	'Callback',@nextstep,...
	'Style','pushbutton');

% Restart button
uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'String','Restart', ...
	'Position',[.85,.86,.1,.05], ...
	'Callback',@restart,...
	'Style','pushbutton');

% Text for error bound reached
AMB_Data.uicontrols(5)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
        'FontWeight','bold', ...
        'ForegroundColor','red', ...
	'Position',[.7,.81,.3,.05],...
	'String','Error bound reached',...
        'visible','off',...
	'Style','text');

% Axis for function display
AMB_Data.axis(1)=subplot('position',[.05,.30,.9,.5]);

% Axis for error display
AMB_Data.axis(2)=subplot('position',[.05,.05,.9,.2]);

% Axis for xi
AMB_Data.axis(3)=subplot('position',[.05,.25,.9,.05]);

% Display first approximation
myplot;

%%% Callback Functions %%%

% new function input
function function_update(varargin)
global AMB_Data;
AMB_Data.Funstring = get(AMB_Data.uicontrols(1),'String');
AMB_Data.Fun= eval(AMB_Data.Funstring);
restart;

% new error bound input
function errorbound_update(varargin)
global AMB_Data;
AMB_Data.errorbound = str2double(get(AMB_Data.uicontrols(2),'String'));
myplot;

% button next pressed
function nextstep(varargin)
global AMB_Data;
if AMB_Data.error <  AMB_Data.errorbound return,end 

xin=zeros(0,1);
xi=AMB_Data.knots;
c=AMB_Data.c;
for k1=3:length(xi)-3
  x=linspace(xi(k1),xi(k1+1),101)';
  p=spline_evaluation(c,xi,x);
  err=max(abs(p-AMB_Data.Fun(x)));
  if err>AMB_Data.errorbound
    xin=[xin;(xi(k1)+xi(k1+1))/2];
  end    
end
xi=sort([xi;xin]);
f1=-AMB_Data.Fun(xi(2:end-1))/2;
f2=AMB_Data.Fun((xi(2:end-2)+xi(3:end-1))/2);
AMB_Data.c=(f1(1:end-1)+f1(2:end)+2*f2);
AMB_Data.knots=xi;
myplot;

% button restart pressed
function restart(varargin)
global AMB_Data;
AMB_Data.knots = [-1;0;0;1;1;2];
f=AMB_Data.Fun([0;1/2;1]);
AMB_Data.c=([f(1);2*f(2)-(f(1)+f(3))/2;f(3)]);
AMB_Data.error=inf;
myplot;

%%% Internal Functions %%%
function myplot;
global AMB_Data;
subplot(AMB_Data.axis(1));cla;hold on;box on;

x=linspace(0,1,AMB_Data.npts)';
xi=AMB_Data.knots;
c=AMB_Data.c;
p=spline_evaluation(c,xi,x);
y=AMB_Data.Fun(x);
y0=min([y;p]);y1=max([y;p]);
axis([0,1,y0-0.1*(y1-y0),y1+0.1*(y1-y0)]);

% plot function and visualize prescribed tolerance
plot(x,y,'-','linewidth',2,'color',[0,0,.6]);
% plot projection
plot(x,p,'k-');
set(gca,'Xtick',[]);

% plot error
subplot(AMB_Data.axis(2));cla;
err=abs(y-p);
semilogy(x,err,'-k');
hold on;
plot([0,1],[1,1]*AMB_Data.errorbound,':k');
set(gca,'Xtick',[0,1]);

% display error and check for bound
AMB_Data.error=max(err);
set(AMB_Data.uicontrols(3),'String',sprintf('Error: %8.3g',AMB_Data.error)),
if AMB_Data.error <  AMB_Data.errorbound 
  set(AMB_Data.uicontrols(5),'visible','on');
  set(AMB_Data.uicontrols(4),'visible','off');
else
  set(AMB_Data.uicontrols(5),'visible','off');
  set(AMB_Data.uicontrols(4),'visible','on');
end 

% plot ticks
subplot(AMB_Data.axis(3));cla;hold on;axis([0,1,-.1,.1]);axis off;
plot([0,1],[0,0],'-k');
for k=xi(3:end-2)'
  plot(k*[1,1],[-.08,.08],'-k');
end

%%% Functions from AMB_Programs (included for standalone use) %%% 

function p = spline_evaluation(c,xi,x) 
% SPLINE_EVALUATION: evaluation of a spline function
% p = spline_evaluation(c,xi,x) 
% values p(j) at x(j) of a spline function  
% with coefficients c and knots xi

% change to column vectors (if necessary)
c = c(:); xi = xi(:); x = x(:);

% degree
n = length(xi)-length(c)-1;

% add points and knots to avoid exceptional cases
c = [zeros(n+1,1); c; zeros(n+1,1)];
xi_max = (2*norm([xi; x],inf)+1)*ones(n+1,1);
xi = [-xi_max; xi; xi_max];

for j = 1:length(x)
   % interval: xi(kj) <= x(j) < xi(kj+1)
   kj = min(find(xi>x(j)))-1;
   % relevant control points
   pxj = c(kj-n:kj);
   % triangular scheme
   for m = 1:n 
      ind = [kj-n+m:kj];
      gamma = ((x(j) - xi(ind)) ./ (xi(ind+n+1-m) - xi(ind)));
      pxj = gamma.*pxj(2:end) + (1-gamma).*pxj(1:end-1);
   end 
   p(j) = pxj;
end
p = p(:); 






