
% SCRIPT: generate and solve random matrix fractional problems

list_n=[];
list_iter=[];
list_std=[];

for p=1:1:20,

n=p*2;

iter_1=[];
for k=1:11,

% random matrix-fractional problem
if ~exist('new'), new=1; end;
if new,
	F=random('Normal',0,1,n,p);
	g=random('Normal',0,1,n,1);
	P=[];
	for i=0:p,
		r=n;	% use this to test P's w/ different ranks
		if i==0,
			V=random('Normal',0,1,n,r);
			% V=zeros(n,1);	% to switch from affine to linear
		else
			V=random('Normal',0,1,n,r);
		end;
		P=[P; V*V'];
	end;
end;

% cast as SOCP
f=[g; 1];
N=(n+1)*ones(p+1,1);
Psqrt=sqrtm(P(1:n,:));
A=[1/2*Psqrt zeros(n,1); zeros(1,n) -1/2];
C=[zeros(1,n) 1/2];
b=[zeros(n,1); 1/2];
d=1/2;
for i=1:p,
	Psqrt=sqrtm(P(i*n+1:i*n+n,:));
	A=[A; 1/2*Psqrt zeros(n,1); 1/2*F(:,i)' 0];
	C=[C; -1/2*F(:,i)' 0];
	b=[b; zeros(n,1); 1/2];
	d=[d; 1/2];
end;

% solve
maxgap=1e-6;
reltol=0;
maxiter=100;
[x,info,z,w,gap,time]=...
	socp(f,A,b,C,d,N,[],[],[],maxgap,reltol,0,maxiter,[]);
iter=sum(time(:,3));

% retrieve solution to original problem
xp=zeros(p,1);
for i=1:p,
	xp(i)=1/2*(z((i+1)*(n+1))-w(i+1));
end;
socp_obj=-f'*x;

iter_1=[iter_1 iter];

end;	% of for k=...

% performance info and plots
list_n=[list_n n];
list_iter=[list_iter mean(iter_1)];
list_std=[list_std std(iter_1)];
figure(1);
if 1,
	errorbar(list_n,list_iter,list_std,list_std,'o');
else
	plot(list_n,list_iter,'o');
end;
axis([0 max(list_n) 0 max(list_iter)+2]);
xlabel('n');
ylabel('i');
grid on;

%% save mf.mat list_n list_iter list_std

end;	% of for p=...



if 0,

if p==2,

% mesh plot of objective around optimum

t=(0.1:0.1:1.9)';
u=xp;
D=zeros(length(t),length(t));
Z=zeros(length(t),length(t));
for i=1:length(t),
	u(1)=t(i)*xp(1);
	for j=1:length(t),
		u(2)=t(j)*xp(2);
		Pu=P(1:n,:);
		for l=1:p,
			Pu=Pu+u(l)*P(l*n+1:l*n+n,:);
		end;
		D(i,j)=min(svd(Pu));
		Z(i,j)=(F*u+g)'*inv(Pu)*(F*u+g);
	end;
end;
m=(length(t)+1)/2;
comp_obj=Z(m,m);
disp('values at minimum =');
disp(Z(m-1:m+1,m-1:m+1));
figure(1);
mesh(t,t,Z);
grid on;

% plot directional derivative as fct. of angle

a=0:5:360;
h=0.1;
dx=zeros(n,length(a));
dx(1,:)=h*cos(a/360*2*pi);
dx(2,:)=h*sin(a/360*2*pi);
u=xp*ones(1,length(a))+dx;
R=zeros(length(a),1);
for i=1:length(a),
	Pu=P(1:n,:);
	for l=1:p,
		Pu=Pu+u(l,i)*P(l*n+1:l*n+n,:);
	end;
	R(i)=(F*u(:,i)+g)'*inv(Pu)*(F*u(:,i)+g);
end;
R=(R-Z(m,m))./sqrt(sum(dx.^2))';
figure(2);
plot(a,R);
grid on;
R=R(find(all(u>=0)));

else	% p>2

% compute directional derivatives in random directions and plot distribution

m=100;
h=1e-3*norm(xp);
D=[];
Z=[];
dx=[];
for i=1:m,
	if i==1,
		u=xp;
	else
		dx(:,i)=random('Normal',0,1,p,1);
		dx(:,i)=dx(:,i)/norm(dx(:,i));
		u=xp+h*dx(:,i);
	end;
	u=abs(u);	% make all components positive
	Pu=P(1:n,:);
	for l=1:p,
		Pu=Pu+u(l)*P(l*n+1:l*n+n,:);
	end;
	D=[D min(svd(Pu))];
	Z=[Z (F*u+g)'*inv(Pu)*(F*u+g)];
end;
D=D(2:m);
comp_obj=Z(1);
R=(Z(2:m)-Z(1))/h;
figure(1);
hist(R);

end;	% of if p==2,

disp('minimum directional derivative =');
disp(min(R));
disp('minimum singular value =');
disp(min(min(D)));

disp('socp objective / computed objective - 1 =');
disp(socp_obj/comp_obj-1);
disp('optimal x =');
disp(xp);
disp('outer iterations =');
disp(iter);

end;	% of if 0,
