% test of KWW 2012
clear all
%define the working deriction
workhouse='E:\Research\GVCs\data\KWW12\BEC-classicification\only_2007\constrainted_optimisation\'
cd(workhouse)

%inmput constraint parameters
% 1.data from IO table
x =xlsread('x.xlsx');       % gross output of product
z =xlsread('z.xlsx');       % goods i used as intermediate inputs in j
v =xlsread('v.xlsx');       % value add in secter j
m =xlsread('m.xlsx');       % value add in secter j
y =xlsread('y.xlsx');       % total final demand except good i

%data from trade statistics
mp =xlsread('m_p.xlsx');
md =xlsread('m_d.xlsx');
ep =xlsread('e_p.xlsx');
en =xlsread('e_d.xlsx');

% data of final use
yd =xlsread('y_d.xlsx');   %final demand of demestic product
ym =xlsread('y_m.xlsx');   %final demand of import product

% optimisation process
A=[];
b=[];                     % without liner constraint%
lb=zeros(32560,1);                    % without underground
ub=[];                    % without ceiling
nonlcon=[] ;              % without non-liner constraint

% Constraint of Optimization
% basic opreator of martix adding by "j"
for i=1:90;
a1(:,:,i)=eye(90);
end
addop1=[]
for i=1:90;
addop1=cat(2,addop1,a1(:,:,i));
end
delop1=zeros(90,8100);
delopv1=zeros(90,180);


% basic otreator of vectors in add by "j"
%---------(5.21)---------%
O_21=[delop1 delop1 addop1 addop1 delopv1];
%---------(5.22)---------%
O_22=[addop1 addop1 delop1 delop1 delopv1];
%---------(5.25)---------%
O_25=[delop1 addop1 delop1 delop1 delopv1];
%---------(5.26)---------%
O_26=[addop1 delop1 delop1 delop1 delopv1];
%---------(5.27)---------%
O_27=[delop1 delop1 addop1 addop1 delopv1];
clear delop1 addop1 delopv1 addopv1;

% basic opreator of martix adding by "i"
%attention : defination of A(:) is very important
I2=eye(90);
E2=ones(1,90);
addop2=kron(I2,E2);        % add   opreator
addopv2=eye(90);
delop2=zeros(90,8100);     % delet opreator
delopv2=zeros(90);


%---------(5.23)---------%
%--------vd then  vp------%
O_23=[delop2 addop2 addop2 delop2 addopv2 delopv2];
%---------(5.24)---------%
O_24=[addop2 delop2 delop2 addop2 delopv2 addopv2];
clear delop2 addop2 delopv2 addopv2 I2 E2;


%---------(5.28)---------%
a2=eye(8100);
b2=zeros(8100,180);
O_28=[a2 a2 a2 a2 b2];
clear a2 b2;


%---------(5.29)--------%
a3=[eye(90) eye(90)];
b3=zeros(90,32400);
O_29=[b3 a3];
clear a3 b3;

Aeq=[O_21;O_22;O_23;O_24;O_25;O_26;O_27;O_28;O_29] ;
clear O_21 O_22 O_23 O_24 O_25 O_26 O_27 O_28 O_29;
beq=[x-ep-en-yd;m-ym;x-ep;ep;md;mp;sum(z,2)-(md-mp);z(:);v];
A=[];
b=[];


options = optimoptions( 'fmincon', 'StepTolerance',1.0000e-100000,...
                        'MaxFunctionEvaluations',100000000,...
                        'display','iter-detailed',...
                        'ConstraintTolerance',1.0000e-08,...
                        'OptimalityTolerance',1.0000e-10,... 
                        'Algorithm','Interior-Point',...
                        'FunctionTolerance',1.0000e-03,...
                        'PlotFcn',{@optimplotfval,@optimplotfirstorderopt},...
                        'HessianApproximation','lbfgs',...
                        'MaxIterations',200)

                    
Z_g=xlsread('guess_value.xlsx');          %the parameters of Z_g
dZ_g=xlsread('guess_dvalue.xlsx');        %the parameters of dZ_g
final=@(Z_t)kww2012(Z_t,Z_g,dZ_g);        %final will be annonymous function of kww12
%Z_t0=zeros(80,322);                      % initial value of Z_t
Z_t0=Z_g; 
%Z_t0=xlsread('fguess_value.xlsx');

[Z_t,fval,exitflag,output]=fmincon(final,Z_t0,A,b,Aeq,beq,lb,ub,nonlcon,options);  


%Z_t = patternsearch(final,Z_t0,A,b,Aeq,beq,lb,ub,nonlcon,options);

% local optimazations

% f = zeros(size(Z_t0)); 
% 
% options = optimoptions('linprog','MaxIterations',50000)
% 
% [Z_t,fval,exitflag]= linprog(f,A,b,Aeq,beq,lb,ub,Z_t0,options);

%[Z_t,fval,exitflag] = ga(final,25760,A,b,Aeq,beq);

%[x,fval,exitflag] = ga(final,25760,A,b,Aeq,beq,lb,ub,nonlcon,options);   

% global optimization with GA

xlswrite('ture_value_without_90.xlsx',Z_t);

% read data from ture value
%Z_t=xlsread('ture_value_handled.xlsx');

Z_t=xlsread('ture_value_without_90.xlsx');

zmp=Z_t(:,1:90);
zmd=Z_t(:,91:180);
zdd=Z_t(:,181:270);
zdp=Z_t(:,271:360);
vd=Z_t(:,361);
vp=Z_t(:,362);

% IO coefficients of new IO model
Adv=vd./(x-ep);
Apv=vp./ep;

alter=x-ep;
for i=1:90
   Add(i,:)=zdd(i,:)./alter';
end

for i=1:90
   Amd(i,:)=zmd(i,:)./alter';
end

for i=1:90
   Amp(i,:)=zmp(i,:)./ep';
end

for i=1:90
   Adp(i,:)=zdp(i,:)./ep';
end

% VSSd and VSSp
u=ones(1,90);

VSSd=(u*Amd*inv(eye(90)-Add))';

VSSp=(u*Amd*inv(eye(90)-Add)*Adp+u*Amp)';

t2=Amd*inv(eye(90)-Add)*Adp;

t1=Amd*inv(eye(90)-Add);

t3=Amp;

tpindre=(u*Amd*inv(eye(90)-Add)*Adp)';

tpdre=(u*Amp)';

tn=(u*Amd*inv(eye(90)-Add))';

tdindre=(u*Amd*inv(eye(90)-Add))'-(u*Amd)';

xlswrite('foreign_content_indrect_p.xlsx',tpindre);
xlswrite('foreign_content_indrect_d.xlsx',tdindre);




