function [solmean,solcvar,alpha,ymean,bias,variance] =  IPtoyRegFunc(eps,C,figno);
%eps = 0.202;
A = [1 1+eps; 1 - eps 1];
%C = eye(2);
logalpha = [-5:0.3:-2];
%logalpha = [-2:0.3:1];
aa = 5*10.^logalpha;
nalpha = length(aa);

[U,W,V] = svd(A);
disp(['det A : ',num2str(det(A))]);
disp(['condition number : ',num2str(W(1,1)/W(2,2))]);
nsample = 2000;
sol = zeros(2,nsample);
dat = zeros(2,nsample);
dat0 = [10;20];
sol0 = inv(A)*dat0;
disp(['noise free solution ',num2str(sol0')]);
disp(['svd of forward operator ',num2str(W(1,1)),' ',num2str(W(2,2))]);
%cols = {'r+','g+','b+','m+','c+'};
cols = {'r+','g+','b+','m+','c+','y+','ro','go','bo','mo','co','yo'};

figure(figno); clf; 
subplot(2,2,1); hold on; plot(dat0(1),dat0(2),'k+');title('data + noise');
subplot(2,2,2); hold on; plot(sol0(1),sol0(2),'k+');title('solutions');
subplot(2,2,3); hold on; plot(dat0(1),dat0(2),'k+');title('projected solution');

solcvar = zeros(2,2,nalpha);
solmean = zeros(2,nalpha);

y = solmean;
for j = 1:nalpha
    alpha = aa(j);
    disp(['regularisation parameter alpha =',num2str(alpha)]);
    for i = 1:nsample
        dat(:,i) = dat0 + randn(2,1);
        sol(:,i) = inv(A'*A + alpha*C)*A' * dat(:,i);
        %   sol(:,i) = inv(A'*A + alpha*[1 10; 10 1])*A' * dat(:,i); % rotated basis
    end
    subplot(2,2,1); hold on; plot(dat(1,:),dat(2,:),cols{j}); 
    subplot(2,2,2); hold on; plot(sol(1,:),sol(2,:),cols{j}); %hold on; plot(sol0(1),sol0(2),'g+');
    %plot(sol0(1)+V(2,1)*[-200,200],sol0(2)+V(2,2)*[-200,200]);
    %plot(sol0(1)+V(1,1)*[-20,20],sol0(2)+V(1,2)*[-20,20]);
    y = A*sol;
    subplot(2,2,3); hold on;  plot(y(1,:),y(2,:),cols{j});
    solmean(:,j) = sum(sol')/nsample;
    ymean(:,j) = sum(y')/nsample;
    solcvar(:,:,j) = sol*sol' /nsample - solmean(:,j) * solmean(:,j)';
    [UP,SP] = eig(solcvar(:,:,j));
    sd = 1./diag(sqrt(SP));
%    disp(['mean solution ',num2str(solmean(:,j)')]);
%    disp(['sd of covariance matrices ',num2str(sd(1)),' ',num2str(sd(2))]);
    subplot(2,2,1); hold on; plot(dat0(1),dat0(2),'k+');
    subplot(2,2,2); hold on; plot(sol0(1),sol0(2),'k+');
    subplot(2,2,3); hold on; plot(dat0(1),dat0(2),'k+');
    bias(j) = norm(sol0-solmean(:,j))^2;
    variance(j) = trace(solcvar(:,:,j));
    subplot(2,2,4); hold on; plot((variance(j)),(bias(j)),cols{j});title('bias-variance');
   
    pause(5);
end
    solmean = [sol0, solmean];
    ymean = [dat0, ymean];
    subplot(2,2,2); hold on; plot([ solmean(1,:)],[ solmean(2,:)],'k');
    subplot(2,2,3); hold on; plot([ ymean(1,:)],[ ymean(2,:)],'k');
    subplot(2,2,4); hold on; plot(variance,bias,'k');
%     figure(figno+1); clf;
%     hold on; loglog(variance,bias,'k');ylabel('bias'); xlabel('variance');
%     for j = 1:nalpha
%         loglog(variance(j),bias(j),cols{j});
%     end