%%% This is an adaption of the Matlab code from 
%%% https://web.stanford.edu/~jhain/synthpage.html
%%% 
%%% From readme_synth_matlab.txt:
%%% Synth MATLAB Code (11/07/2006) written for MATLAB 7.0
%%% by Alberto Abadie, Alexis Diamond, and Jens Hainmueller (all Harvard University)
%%% Contact: jhainm@harvard.edu

clear all;
diary main;

%% Get Data
load main_X0.txt;
load main_X1.txt;
load main_Z0.txt;
load main_Z1.txt;
X0 = main_X0;
X1 = main_X1;
Z0 = main_Z0;
Z1 = main_Z1;

% Normalization (probably could be done more elegantly)
bigdata = [X0,X1];
divisor = std(bigdata');
scamatrix = (bigdata' * diag(( 1./(divisor) * eye(size(bigdata,1))) ))';
X0sca = scamatrix([1:size(X0,1)],[1:size(X0,2)]);
X1sca = scamatrix(1:size(X1,1),[size(scamatrix,2)]);
X0 = X0sca;
X1 = X1sca;
clear divisor X0sca X1sca scamatrix bigdata;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now we implement Optimization

% Check and maybe adjust optimization settings if necessary
options = optimset('fmincon')

% Get Starting Values 
s = std([X1 X0]')';
s2 = s; s2(1)=[];
s1 = s(1);
v20 =((s1./s2).^2);

[v2,fminv,exitflag] = fmincon('loss_function',v20,[],[],[],[],...
   zeros(size(X1)),[],[],options,X1,X0,Z1,Z0);
display(sprintf('%15.4f',fminv));
v = [1;v2];
% V-weights
v

% Now recover W-weights
D = diag(v);
H = X0'*D*X0;
f = - X1'*D*X0;
options = optimset('quadprog')
[w,fval,e]=quadprog(H,f,[],[],ones(1,length(X0)),1,zeros(length(X0),1),ones(length(X0),1),[],options);
w = abs(w); 

% W-weights
w
save -ascii main_w.txt w
save -ascii main_v.txt v
loss = fminv/size(Z0,1);
save -ascii main_loss.txt loss

diary off;
