% Simple test
% d=2;
% p=[2;1];
% a=[1;1;1];
% m=[2.1;2.3;2.7];
% c=[3;3;3];
% lb=[1;1;1];
% ub=[2;2;2];
% t=sum(c)*1.4;
%
% Test from Mathis and Mathis
%
% Original data
oij=[16.35;74.12;39.26;48.20;12.63;...
  13.22;12.67;22.71;5.44;7.24;...
  12.66;6.84;3.72;6.66;15.89;...
  33.03;0.18;5.23;7.60;67.54];
mij=[48.54;10.11;2.34;5.69;7.20;...
  48.45;20.94;64.95;22.88;17.49;...
  13.60;6.97;22.39;55.18;63.44;...
  14.04;7.07;1.33;41.53;23.74];
cij=[64.89;84.23;41.60;53.89;19.83;...
  61.66;33.61;87.66;28.32;24.73;...
  26.26;13.81;26.11;61.84;79.33;...
  47.07;7.25;6.56;49.13;91.28];
Ci=[185.68;118.90;134.21;166.56];
q=0.10;
d=4;
p=[5;5;5;5];
nb=sum(p);
rij=zeros(nb,1);
uij=ones(nb,1)*0.15;
lij=ones(nb,1)*(-0.05);
% Mathis' solution
r=[.15;.15;.15;.15;.15;...
  .05;.15;.15;-.05;.15;...
  .15;.15;-.05;-.05;-.05;...
  .15;.15;.15;.15;.15];
x=ones(nb,1)+r;

f=oij'*x + ...
	Ci(1)*((mij(1:5)'*x(1:5))/(cij(1:5)'*x(1:5)))+...
	Ci(2)*((mij(6:10)'*x(6:10))/(cij(6:10)'*x(6:10)))+...
	Ci(3)*((mij(11:15)'*x(11:15))/(cij(11:15)'*x(11:15)))+...
	Ci(4)*((mij(16:20)'*x(16:20))/(cij(16:20)'*x(16:20)));

disp('Mathis optimal value');
f
disp('Mathis feasibility...next value should be zero...');
cij'*r-q*sum(cij)


%
% Transformed data
%
c=cij;
m=[mij(1:5)*Ci(1);mij(6:10)*Ci(2);mij(11:15)*Ci(3);mij(16:20)*Ci(4)];
a=oij;
lb=ones(nb,1)+lij;
ub=ones(nb,1)+uij;
t=(1+q)*sum(cij);

optx=plsimplex(d,p,c,m,a,lb,ub,t)
% solution in original form...
rx=optx-ones(nb,1);

disp('My optimal value');
x=optx;
r=optx-ones(nb,1);
f=oij'*x + ...
	Ci(1)*((mij(1:5)'*x(1:5))/(cij(1:5)'*x(1:5)))+...
	Ci(2)*((mij(6:10)'*x(6:10))/(cij(6:10)'*x(6:10)))+...
	Ci(3)*((mij(11:15)'*x(11:15))/(cij(11:15)'*x(11:15)))+...
	Ci(4)*((mij(16:20)'*x(16:20))/(cij(16:20)'*x(16:20)))
disp('My feasibility...next value should be zero...');
cij'*r-q*sum(cij)

