N = 6;
[x,y] = meshgrid(-10:10, -13:13);
x = x/max(x(:));
y = y/max(y(:));
a = randn(1,(N+1)*(N+2)/2).*(1:((N+1)*(N+2)/2))/N^2;
fx =  a(1)*x + a(2)*y + ...
      a(3)*x.^2 +  a(4)*x.*y + a(5)*y.^2 + ...
      a(6)*x.^3 +  a(7)*x.^2.*y + a(8)*x.*y.^2 + a(9)*y.^3 + ...
      a(10)*x.^4 + a(11)*x.^3.*y + a(12)*x.^2.*y.^2 + a(13)*x.*y.^3 + a(14)*y.^4 ...
      a(15)*x.^5 + a(16)*x.^4.*y + a(17)*x.^3.*y.^2  + a(18)*x.^2.*y.^3 + a(19)*x.*y.^4 + a(20)*y.^5 ;
%     a(20)*x.^6 + a(21)*x.^5.*y + a(22)*x.^4.*y.^2 + a(23)*x.^3.*y.^3 + a(24)*x.^2.*y.^4 + a(25)*x.*y.^5 + a(26)*y.^6;

b = randn(1,(N+1)*(N+2)/2).*(1:((N+1)*(N+2)/2))/N^2;
fy = b(1)*x + b(2)*y + ...
     b(3)*x.^2 + b(4)*x.*y + b(5)*y.^2 + ...
     b(6)*x.^3 + b(7)*x.^2.*y + b(8)*x.*y.^2 + b(9)*y.^3 + ...
     b(10)*x.^4 + b(11)*x.^3.*y + b(12)*x.^2.*y.^2 + b(13)*x.*y.^3 + b(14)*y.^4 + ...
     b(15)*x.^5 + b(16)*x.^4.*y + b(17)*x.^3.*y.^2 + b(18)*x.^2.*y.^3 + b(19)*x.*y.^4 + b(20)*y.^5;
%     b(21)*x.^6 + b(22)*x.^5.*y + b(23)*x.^4.*y.^2 + b(24)*x.^3.*y.^3 + b(25)*x.^2.*y.^4 + b(26)*x.*y.^5 + a(27)*y.^6;

fx = fx + 0*randn(size(fx))/400;
fy = fy + 0*randn(size(fy))/400;
[ux,sx,vx] = svd(fx);
[uy,sy,vy] = svd(fy);

%% resample the columns up to the highest valid singular value
Nx = select_valid_singval(diag(sx));
Ny = select_valid_singval(diag(sy));
Nx = max(Nx,Ny);
Ny = Nx;

Ux = zeros(size(ux,1)*32, Nx);
Vx = zeros(size(vx,1)*32, Nx);

Uy = zeros(size(uy,1)*32, Ny);
Vy = zeros(size(vy,1)*32, Ny);

# Interpolated vx into Vx and Ux into Ux
yi = linspace(x(1,1), x(1,end), size(Ux,1));
xi = linspace(y(1,1), y(end,1), size(Vx,1));

# interpolate all columns of vx and ux
for n=1:Nx
  [p,r] = polyfit(y(:,1)', ux(:,n), Nx);
  Ux(:,n) = polyval(p, yi);
endfor
# interpolate all columns of vx and ux
for n=1:Nx
   [p,r] = polyfit(x(1,:)', vx(:,n), Nx);
   Vx(:,n) = polyval(p, xi);
endfor
fxi = Ux(:,1:Nx)*diag(diag(sx)(1:Nx))*(Vx(:,1:Nx))';
for n=1:Ny
  [p,r] = polyfit(y(:,1)', uy(:,n), Ny);
  Uy(:,n) = polyval(p, yi);
endfor
# interpolate all columns of vx and ux
for n=1:Ny
   [p,r] = polyfit(x(1,:)', vy(:,n), Ny);
   Vy(:,n) = polyval(p, xi);
endfor
fyi = Uy(:,1:Ny)*diag(diag(sy)(1:Ny))*(Vy(:,1:Ny))';
[x,y] = meshgrid(linspace(-1,1,size(fxi,2)), linspace(-1,1,size(fxi,1)));
x = x/max(x(:));
y = y/max(y(:));
## Recreate the resolved input without noise to get a benchmark
fxo = a(1)*x + a(2)*y + ...
      a(3)*x.^2 +  a(4)*x.*y + a(5)*y.^2 + ...
      a(6)*x.^3 +  a(7)*x.^2.*y + a(8)*x.*y.^2 + a(9)*y.^3 + ...
     a(10)*x.^4 + a(11)*x.^3.*y + a(12)*x.^2.*y.^2 + a(13)*x.*y.^3 + ...
     a(14)*x.^5 + a(15)*x.^4.*y + a(16)*x.^3.*y.^2  + a(17)*x.^2.*y.^3 + a(18)*x.*y.^4 + a(19)*y.^5 ;
fyo = b(1)*x + b(2)*y + ...
     b(3)*x.^2 + b(4)*x.*y + b(5)*y.^2 + ...
     b(6)*x.^3 + b(7)*x.^2.*y + b(8)*x.*y.^2 + b(9)*y.^3 + ...
     b(10)*x.^4 + b(11)*x.^3.*y + b(12)*x.^2.*y.^2 + b(13)*x.*y.^3 + ...
     b(14)*x.^5 + b(15)*x.^4.*y + b(16)*x.^3.*y.^2  + b(17)*x.^2.*y.^3 + b(18)*x.*y.^4 + b(19)*y.^5;

##subplot(1,3,1),plot(log(diag(sx))); hold on,plot(log(diag(sy))),hold off,grid on
##subplot(1,3,2), imagesc((fxi - fxo)/max(fxo(:))),axis equal,colorbar
##subplot(1,3,3), imagesc((fyi - fyo)/max(fyo(:))),axis equal,colorbar
e = [e std(fxi(:) - fxo(:))/std(fxo(:))];
e = [e std(fyi(:) - fyo(:))/std(fyo(:))];
