function X_complete = spoc(X,M,r)
% 
% SPOC 2.0  
% January, 2008
%
% ------------------------------------------------------------------------
% Copyright (2007): Pedro M.Q. Aguiar, Joćo  M.F. Xavier, Marko Stosic
%
% Permission to use, copy, modify, and distribute this software for any
% purpose without fee is hereby granted, provided that this entire notice
% is included in all copies of any software which is or includes a copy or
% modification of this software and in all copies of the supporting
% documentation for such software. This software is provided "as is",
% without any express or implied warranty. In particular, the authors do
% not make any representation or warranty of any kind concerning the
% merchantability of this software or its fitness for any particular
% purpose.
% ------------------------------------------------------------------------
% 
% ------------------------------------------------------------------------
% This MATLAB file contains an implementation of the SPOC (SPectrally
% Optimal Completion) algorithm described and illustrated in the papers:
% (available at http://users.isr.ist.utl.pt/~aguiar/pubs.html)
%
% "Spectrally optimal factorization of incomplete matrices",
% by Pedro M. Q. Aguiar, Marko Stosic, and Joćo Xavier,
% IEEE Conf. on Computer Vision and Pattern Recognition CVPR’08,
% Anchorage, Alaska, USA, June 2008
%
% "Globally optimal solution to exploit rigidity when recovering structure
% from motion under occlusion",
% Pedro M. Q. Aguiar, Joćo M. F. Xavier, and Marko Stosic,
% IEEE Int. Conf. on Image Processing ICIP’08,
% S. Diego CA, USA, October 2008
%
% "On singular values of partially prescribed matrices"
% Pedro M. Q. Aguiar, Marko Stosic, and Joćo Xavier
% ELSEVIER Linear Algebra and its Applications,
% vol. 429, 8-9, October 2008
% ------------------------------------------------------------------------
%
% ------------------------------------------------------------------------
% This code is (permanently?) in development stage; any comments or bug
% reports are very welcome.
% Contact: aguiar@isr.ist.utl.pt
% ------------------------------------------------------------------------
%
% ------------------------------------------------------------------------
% X_complete = spoc(X,M,r)
%
% Input:
% X - incomplete matrix (missing entries form a top left Young diagram).
%     Ex: X=[?   ?   ?   ?   ?   x16 x17 x18
%            ?   ?   ?   ?   x25 x26 x27 x28
%            ?   ?   x33 x34 x35 x36 x37 x38
%            ?   x42 x43 x44 x45 x46 x47 x48
%            x51 x52 x53 x54 x55 x56 x57 x58
%            x61 x62 x63 x64 x65 x66 x67 x68]
% M - Young diagram mask (binary matrix, where: 1-known, 0-unknow).
%     Ex: M=[0   0   0   0   0   1   1   1
%            0   0   0   0   1   1   1   1
%            0   0   1   1   1   1   1   1
%            0   1   1   1   1   1   1   1
%            1   1   1   1   1   1   1   1
%            1   1   1   1   1   1   1   1]
% r - order of singular value to minimize.
%
% Output:
% X_complete - completion of X with minimum sigma_r.

[m,n] = size(X);
X_complete=X;
for i = m:-1:1
    for j = n:-1:1        
        if M(i,j) == 0
     %       fprintf('(%d,%d): ',i,j);
            x = one_missing(X(i:end,j:end),r);
            X_complete(i,j) = x;
        end;
    end;
end;


function x = one_missing(M,r)

% input
% M = [ ? v“ ;  u A ] with missing entry in (1,1)
% r = order of the singular value

% output
% x is optimal for min_t sigma_r([ t v' ; u A]) 

[m,n] = size(M);

%fprintf('Size %d,%d ',m,n);
u = M(2:end,1);
v = M(1,2:end); v = v';
A = M(2:end,2:end);

B = [ u A ];
sB = svd(B);

if r <= length(sB)
    boundB = sB(r);
else boundB = 0;
end;

C = [ v' ; A ];
sC = svd(C);

if r <= length(sC)
    boundC = sC(r);
else boundC = 0;
end;

s = max([ boundB boundC ])+0.0001;



xs = -4:0.1:4;
z = [];
for aux = 1:length(xs)
    x_aux = xs(aux);
    z = [ z , det([ x_aux^2 + v'*v-s^2 x_aux*u'+v'*A' ; x_aux*u+A*v u*u'+A*A'-s^2*eye(size(A,1)) ]) ];
end;
%figure(2); clf; plot(xs,z); 
    
x_aux = 0;
p0 = det([ x_aux^2 + v'*v-s^2 x_aux*u'+v'*A' ; x_aux*u+A*v u*u'+A*A'-s^2*eye(size(A,1)) ]);

x_aux = -1;
pm1 = det([ x_aux^2 + v'*v-s^2 x_aux*u'+v'*A' ; x_aux*u+A*v u*u'+A*A'-s^2*eye(size(A,1)) ]);

x_aux = 1;
p1 = det([ x_aux^2 + v'*v-s^2 x_aux*u'+v'*A' ; x_aux*u+A*v u*u'+A*A'-s^2*eye(size(A,1)) ]);

c = p0;
a = 0.5*(p1+pm1)-p0;
b = 0.5*(p1-pm1);

x = -b/(2*a); 


% s_aux = svd( [ x v' ; u A ] );
% if r <= length(s_aux)
%     achieved = s_aux(r);
% else achieved = 0;
% end;
% 
% 
% %fprintf('[ bound = %f, achieved = %f, optimal = %f ]\n',s,achieved,x);
% 
% t = -4:0.1:4;
% y = [];
% for aux = 1:length(t) sz = svd([ t(aux) v' ; u A ]); if r <= length(sz) new = sz(r); else new = 0; end;  y = [ y , new ]; end; 
% %figure(1); clf; plot(t,y);
% %pause;
%