function fig8_5

%  poly fit: normal vs QR, regularization, and overfitting plots
%  also include figures 8.6, 8.7, and 8.8
%  also inclues table 8.2

% use the given training and testing data:  icase=1
% use given training data, but create and save testing data set:  icase=2
% create, and save, training and testing data:  icase=3

icase=1

lambda=1e-4

if icase==1

    X=importdata('training_data.txt');
    xd=X(:,1)';  yd=X(:,2)';
    nd=length(xd);

    X2=importdata('testing_data.txt');
    xd2=X2(:,1)';  yd2=X2(:,2)';
    nd2=length(xd2);

elseif icase==2

    X=importdata('training_data.txt');
    xd=X(:,1)';  yd=X(:,2)';
    nd=length(xd);

    %  generate testing data
    n=12;
    xs=linspace(0,1,n);
    ys=[2 8 8 8 3 2 1 1 1 5 5 5];
    xd2=sort(rand(1,nd));
    yd2=nspline(xs,ys,xd2);
    dis=1+0.5*(rand(1,nd)-0.5);
    yd2=dis.*yd2;
    X2=[xd2' yd2'];
    save('/Users/mark/Desktop/testing_data.txt','X2','-ascii','-double')

elseif icase==3

    %  generate spline
    n=12;
    xs=linspace(0,1,n);
    ys=[2 8 8 8 3 2 1 1 1 5 5 5];

    %  generate data for regression
    nd=20;
    xd=sort(rand(1,nd));
    yd=nspline(xs,ys,xd);
    dis=1+0.2*(rand(1,nd)-0.5);
    yd=dis.*yd;
    X=[xd' yd'];
    save('/Users/mark/Desktop/training_data.txt','X','-ascii','-double')

    %  generate testing data
    xd2=sort(rand(1,nd));
    yd2=nspline(xs,ys,xd2);
    dis=1+0.5*(rand(1,nd)-0.5);
    yd2=dis.*yd2;
    X2=[xd2' yd2'];
    save('/Users/mark/Desktop/testing_data.txt','X2','-ascii','-double')

end

% pseudo-inverse and Rm
m=4;
A=mat(xd,nd,m);
c4=fliplr((inv(A'*A)*A'*yd')');
[Q,R]=qr(A);
Rm=R(1:m,1:m);
Qy=Q'*yd';
ym=Qy(1:m,1);
r4=fliplr((inv(Rm)*ym)');
condc4=cond(A'*A,inf);
condr4=cond(Rm,inf);

m=8;
A=mat(xd,nd,m);
c8=fliplr((inv(A'*A)*A'*yd')');
[Q,R]=qr(A);
Rm=R(1:m,1:m);
Qy=Q'*yd';
ym=Qy(1:m,1);
r8=fliplr((inv(Rm)*ym)');
condc8=cond(A'*A,inf);
condr8=cond(Rm,inf);

m=12;
A=mat(xd,nd,m);
c12=fliplr((inv(A'*A)*A'*yd')');
[Q,R]=qr(A);
Rm=R(1:m,1:m);
Qy=Q'*yd';
ym=Qy(1:m,1);
r12=fliplr((inv(Rm)*ym)');
condc12=cond(A'*A,inf);
condr12=cond(Rm,inf);

cR12=fliplr((inv(A'*A+lambda*eye(12))*A'*yd')');
condcR12=cond(A'*A+lambda*eye(12),inf);

m=16;
A=mat(xd,nd,m);
cc16=(inv(A'*A)*A'*yd')';
c16=fliplr(cc16);
[Q,R]=qr(A);
Rm=R(1:m,1:m);
Qy=Q'*yd';
ym=Qy(1:m,1);
cond(Rm,inf);
rr16=(inv(Rm)*ym)';
r16=fliplr(rr16);
condc16=cond(A'*A,inf);
condr16=cond(Rm,inf);

cc16(1:4)
rr16(1:4)

cR16=fliplr((inv(A'*A+lambda*eye(16))*A'*yd')');
condcR16=cond(A'*A+lambda*eye(16),inf);

%  nn = number of points plotted
nn=200;
X = linspace(xd(1),xd(nd),nn);
Y4 = polyval(c4,X);
R4 = polyval(r4,X);
Y8 = polyval(c8,X);
R8 = polyval(r8,X);
Y12 = polyval(c12,X);
R12 = polyval(r12,X);
Y16 = polyval(c16,X);
R16 = polyval(r16,X);

YR12 = polyval(cR12,X);
YR16 = polyval(cR16,X);

% condition numbers
fprintf('%d &  %5.2e  & %5.2e \n',4,condc4,condr4)
fprintf('%d &  %5.2e  & %5.2e \n',8,condc8,condr8)
fprintf('%d &  %5.2e  & %5.2e \n',12,condc12,condr12)
fprintf('%d &  %5.2e  & %5.2e \n',16,condc16,condr16)

figure(1)
clf
% get(gcf)
set(gcf,'Position', [27 978 666 367])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
marker=18;
m_size=16;

clf
xc=0.12; yc=0;
subaxis(2,2,1,1,'MT',0.005,'MB',0.1,'MR',-0.01,'ML',0.04,'P',0.02)
% plot interpolation fn's and exact result
hold on
plot(X,R4,'r','LineWidth',1.6)
plot(X,Y4,'color',[0 0.5 0],'LineWidth',1.6)
text(xc, yc,['m = 4'],'FontSize',m_size,'FontWeight','bold')
plot(xd,yd,'.b','MarkerSize',marker)
axis([0 1 -3 12])
%xlabel('x-axis')
ylabel('y-axis')
set(gca,'XTick',[0 0.2 0.4 0.6 0.8 1.0])
set(gca,'YTick',[0 6 12])
grid on
set(gca,'FontSize',14,'FontWeight','bold')

subaxis(2,2,2,1)
hold on
plot(X,R8,'r','LineWidth',1.6)
plot(X,Y8,'color',[0 0.5 0],'LineWidth',1.6)
plot(xd,yd,'.b','MarkerSize',marker)
text(xc, yc,['m = 8'],'FontSize',m_size,'FontWeight','bold')
%xlabel('x-axis')
axis([0 1 -3 12])
set(gca,'XTick',[0 0.2 0.4 0.6 0.8 1.0])
set(gca,'YTick',[0 6 12])
grid on
set(gca,'FontSize',14,'FontWeight','bold')

subaxis(2,2,1,2)
hold on
plot(X,R12,'r','LineWidth',1.6)
plot(X,Y12,'color',[0 0.5 0],'LineWidth',1.6)
plot(xd,yd,'.b','MarkerSize',marker)
text(xc, yc,['m = 12'],'FontSize',m_size,'FontWeight','bold')
xlabel('x-axis')
axis([0 1 -3 12])
ylabel('y-axis')
set(gca,'XTick',[0 0.2 0.4 0.6 0.8 1.0])
set(gca,'YTick',[0 6 12])
grid on
set(gca,'FontSize',14,'FontWeight','bold')

subaxis(2,2,2,2)
hold on
plot(X,R16,'r','LineWidth',1.6)
plot(X,Y16,'color',[0 0.5 0],'LineWidth',1.6)
plot(xd,yd,'.b','MarkerSize',marker)
text(xc, yc,['m = 16'],'FontSize',m_size,'FontWeight','bold')
xlabel('x-axis')
axis([0 1 -3 12])
set(gca,'XTick',[0 0.2 0.4 0.6 0.8 1.0])
set(gca,'YTick',[0 6 12])
grid on
set(gca,'FontSize',14,'FontWeight','bold')

%exportgraphics(gcf,'/Users/mark/Desktop/p1.eps')


figure(2)

for im=2:16
    m=im;
    A=mat(xd,nd,m);
    cR=fliplr((inv(A'*A+lambda*eye(im))*A'*yd')');
    ELS_T(im-1)=norm(yd-polyval(cR,xd),2)^2;
    ELS_R(im-1)=norm(yd2-polyval(cR,xd2),2)^2;

    % check on differences in x^{12} coefficient
    if im==12 || im==16
        cR
        cR(m-11)
    end

end

% error for training set
for im=1:15
    E(im)=norm(yd-polyval(polyfit(xd,yd,im),xd),2)^2;
end

clf
% get(gcf)
set(gcf,'Position', [27 674 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.16,'MR',-0.02,'ML',0.04,'P',0.04)

hold on
plot(2:16,ELS_T,'--ob','LineWidth',1.6,'MarkerSize',8)
plot(2:16,ELS_R,'--or','LineWidth',1.6,'MarkerSize',8)
legend(' Training',' Testing','Location','North','AutoUpdate','off','FontSize',16,'FontWeight','bold')
xlabel('m-axis')
ylabel('Error')
grid on
set(gca,'FontSize',16,'FontWeight','bold')

%exportgraphics(gcf,'/Users/mark/Desktop/p2.eps')


figure(3)

% error for training set
for im=1:15
    E(im)=norm(yd-polyval(polyfit(xd,yd,im),xd),2)^2;
end

% error for testing set
for im=1:15
    EE(im)=norm(yd2-polyval(polyfit(xd,yd,im),xd2),2)^2;
end
marker=20;
clf
set(gcf,'Position', [27 232 666 367])
%subaxis(1,1,1,1,'MT',0.01,'MB',0.16,'MR',-0.02,'ML',0.04,'P',0.04)
subaxis(2,1,1,1,'MT',-0.01,'MB',0.07,'MR',-0.03,'ML',0.03,'P',0.045)
hold on
plot(xd,yd,'.b','MarkerSize',20)
plot(xd2,yd2,'+r','MarkerSize',12,'LineWidth',2)
legend(' Training Data',' Testing Data','Location','NorthEast','AutoUpdate','off','FontSize',16,'FontWeight','bold')
xlabel('x-axis')
ylabel('y-axis')
grid on
set(gca,'FontSize',14,'FontWeight','bold')

subaxis(2,1,1,2)
hold on
plot(2:16,E,'-ob','LineWidth',1.6,'MarkerSize',6)
plot(2:16,EE,'-+r','LineWidth',1.6,'MarkerSize',10)
legend(' Training',' Testing','Location','NorthEast','AutoUpdate','off','FontSize',16,'FontWeight','bold')
axis([2 16 0 200])
xlabel('m-axis')
ylabel('Error')
%set(gca,'XTick',[0 0.2 0.4 0.6 0.8 1.0])
grid on
set(gca,'FontSize',14,'FontWeight','bold')

%exportgraphics(gcf,'/Users/mark/Desktop/p3.eps')


figure(4)
clf
% get(gcf)
set(gcf,'Position', [695 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
marker=18;
m_size=16;

clf
xc=0.12; yc=0;
subaxis(1,2,1,1,'MT',0.005,'MB',0.18,'MR',-0.01,'ML',0.04,'P',0.02)
% plot interpolation fn's and exact result
hold on
plot(X,YR12,'k','LineWidth',1.6)
plot(X,R12,'r','LineWidth',1.6)
text(xc, yc,['m = 12'],'FontSize',m_size,'FontWeight','bold')
plot(xd,yd,'.b','MarkerSize',marker)
axis([0 1 -3 12])
xlabel('x-axis')
ylabel('y-axis')
set(gca,'XTick',[0 0.2 0.4 0.6 0.8 1.0])
set(gca,'YTick',[0 6 12])
grid on
set(gca,'FontSize',14,'FontWeight','bold')

subaxis(1,2,2,1)
hold on
plot(X,YR16,'k','LineWidth',1.6)
plot(X,R16,'r','LineWidth',1.6)
text(xc, yc,['m = 16'],'FontSize',m_size,'FontWeight','bold')
plot(xd,yd,'.b','MarkerSize',marker)
axis([0 1 -3 12])
xlabel('x-axis')
%ylabel('y-axis')
set(gca,'XTick',[0 0.2 0.4 0.6 0.8 1.0])
set(gca,'YTick',[0 6 12])
grid on
set(gca,'FontSize',14,'FontWeight','bold')

%exportgraphics(gcf,'/Users/mark/Desktop/p4.eps')

function g=mat(x,n,m)
g=zeros(n,m);
for i=1:n
    for j=1:m
        g(i,j)=x(i)^(j-1);
    end
end



















