% Script C8
% Checks out NearestKP

clc
randn('seed',0)
rand('seed',0)

% Small Example...
p = 7; q = 6;
A = randn(p*q,p*q);
A = A.*(rand(p*q,p*q)<=.3);
A = sparse(A);
spy(A)
shg
pause(1)
opts = struct('disp',0);
[B1,C1,B2,C2] = NearestKP(A,p,q,opts);
clc
format short
B1 = B1
C1 = C1
B2 = B2
C2 = C2
Residual = norm(A - kron(B1,C1) - kron(B2,C2))

% Larger example...
A = delsq(numgrid('H',90));
[n,n] = size(A);
spy(A)
shg
pause(1)
opts = struct('disp',0);
[B1,C1,B2,C2] = NearestKP(A,50,64,opts);
threshold = .001;
close all
B1 = sparse(B1.*(abs(B1)>threshold));
figure
spy(B1)
title('B1')
C1 = sparse(C1.*(abs(C1)>threshold));
figure
spy(C1)
title('C1')
B2 = sparse(B2.*(abs(B2)>threshold));
figure
spy(B2)
title('B2')
C2 = sparse(C2.*(abs(C2)>threshold));
figure
spy(C2)
title('C2')