function [qest1roots, qest2roots, qest3roots, qest4roots] = lieven4(C1, C2)
%------------------------------------------------------------------------------
%	lieven4.m
%	Lieven De Lathauwer's four methods for tensor decomposition based on two
%   tensor C_31 and C_22.
%
% Usage:
%	 [qest1roots, qest2roots, qest3roots, qest4roots] = lieven4(C1, C2);
% Where
%	C1        : The tensor C_31.
%	C2        : The tensor C_22.
% 
%  Modified by Binning Chen on June 12, 2001.
%-------------------------------------------------------------------------------
      
% METHOD 2: ALS ITERATIONS INITIALIZED WITH [1 0]

%C1 = cum40(OBS); C2 = cum22(OBS); % defined before method 1

%C1=C_31;
%C2=TC;

ii=sqrt(-1);

GAMMA1 = ...
[C1(1,1,1,1) C1(1,2,1,1) C1(2,2,1,1) C1(1,1,2,1) C1(1,2,2,1) C1(2,2,2,1)
 C1(1,1,1,2) C1(1,2,1,2) C1(2,2,1,2) C1(1,1,2,2) C1(1,2,2,2) C1(2,2,2,2)];
 
 
GAMMA2 = ...
[C2(1,1,1,1) C2(1,2,1,1) C2(2,2,1,1) C2(1,1,2,1) C2(1,2,2,1) C2(2,2,2,1)
 C2(1,1,1,2) C2(1,2,1,2) C2(2,2,1,2) C2(1,1,2,2) C2(1,2,2,2) C2(2,2,2,2)]; 
 
G11 = GAMMA1(:,1:3); G12  = GAMMA1(:,4:6);
G21 = GAMMA2(:,1:3); G22  = GAMMA2(:,4:6);

x = [1; 0]; Qest1 = [1 0;0 1;0 0;0 0]; sQest1 = 1;
i1 = 1; i1max = 1000; sQestmin = 1e-8;

while i1 < i1max & sQest1 > sQestmin
 [U,S,V] = svd( [x(1)*G11+x(2)*G12; conj(x(1))*G21+conj(x(2))*G22] );
 y = V(:,3);
 [U,S,V] = svd( [G11*y G12*y; conj(G21*y) conj(G22*y)] );
 x = V(:,2);
 qest1 = conv(x,y);
 Qest1 = [Qest1(:,2) qest1];
 SQest1 = svd(Qest1);
 sQest1 = SQest1(2);
 i1 = i1+1;
end

q1est = kron(x,y);

qest1roots = sort(roots(qest1))

% METHOD 3: ALS ITERATIONS INITIALIZED WITH [0 1]

x = [0; 1]; Qest2 = [1 0;0 1;0 0;0 0]; sQest2 = 1;
i1 = 1; i1max = 1000; sQestmin = 1e-8;

while i1 < i1max & sQest2 > sQestmin
 [U,S,V] = svd( [x(1)*G11+x(2)*G12; conj(x(1))*G21+conj(x(2))*G22] );
 y = V(:,3);
 [U,S,V] = svd( [G11*y G12*y; conj(G21*y) conj(G22*y)] );
 x = V(:,2);
 qest2 = conv(x,y);
 Qest2 = [Qest2(:,2) qest2];
 SQest2 = svd(Qest2);
 sQest2 = SQest2(2);
 i1 = i1+1;
end

q2est = kron(x,y);
qest2roots = sort(roots(qest2))

% METHOD 4: FOR A GIVEN XINIT, COMPUTE Y AND X (ALS)
%	    MINIMIZE THE VALUE OF THE EUSIPCO'00 COST FUNCTION 
%	    CORRESPONDING TO (X,Y) WITH fminsearch. 
%	    USE THE OPTIMUM TO INITIALIZE EXTRA ALS ITERATIONS

OLDOPTS = [];
OPTIONS = optimset(OLDOPTS,'display','off');
x_param = fminsearch('obj_sylv_p2_2v2',[0 0],OPTIONS,GAMMA1,GAMMA2);
		% REMARK: xinit = [1 0]
x = [cos(x_param(1)); sin(x_param(1))*exp(ii*x_param(2))];

Qest3 = [1 0;0 1;0 0;0 0]; sQest3 = 1;
i1 = 1; i1max = 1000; sQestmin = 1e-8;

while i1 < i1max & sQest3 > sQestmin
 [U,S,V] = svd( [x(1)*G11+x(2)*G12; conj(x(1))*G21+conj(x(2))*G22] );
 y = V(:,3);
 [U,S,V] = svd( [G11*y G12*y; conj(G21*y) conj(G22*y)] );
 x = V(:,2);
 qest3 = conv(x,y);
 Qest3 = [Qest3(:,2) qest3];
 SQest3 = svd(Qest3);
 sQest3 = SQest3(2);
 i1 = i1+1;
end

q3est = kron(x,y);	% norm(q3est) = 1;
qest3roots = sort(roots(qest3))

% METHOD 5: FOR A GIVEN XINIT, COMPUTE Y (ALS)
%	    MINIMIZE THE VALUE OF THE EUSIPCO'00 COST FUNCTION 
%	    CORRESPONDING TO (XINIT,Y) WITH fminsearch. 
%	    USE THE OPTIMUM TO INITIALIZE EXTRA ALS ITERATIONS


x_param = fminsearch('obj_sylv_p2_3v2',[0 0],OPTIONS,GAMMA1,GAMMA2);
		% REMARK: xinit = [1 0]
x = [cos(x_param(1)); sin(x_param(1))*exp(ii*x_param(2))];

Qest4 = [1 0;0 1;0 0;0 0]; sQest4 = 1;
i1 = 1; i1max = 1000; sQestmin = 1e-8;

while i1 < i1max & sQest4 > sQestmin
 [U,S,V] = svd( [x(1)*G11+x(2)*G12; conj(x(1))*G21+conj(x(2))*G22] );
 y = V(:,3);
 [U,S,V] = svd( [G11*y G12*y; conj(G21*y) conj(G22*y)] );
 x = V(:,2);
 qest4 = conv(x,y);
 Qest4 = [Qest4(:,2) qest4];
 SQest4 = svd(Qest4);
 sQest4 = SQest4(2);
 i1 = i1+1;
end

q4est = kron(x,y);	% norm(q4est) = 1;
qest4roots = sort(roots(qest4))