clear; clc; t1=tic;

%------------- Set parameters --------------------- 
x_num=10;  % The number of atoms in approximating distribution
y_num=100; %200; %500 % The number of atoms in target distribution to be used
iprint=0;       % iprint = 1/0 print or not output txt files
FileName='Solving_process_Log.txt'; % Output Log file

%------------- Load coordinates of atoms and generate probabilities -------------

coord=load('Stan_coordinates_ab_samples_Aquino_50000x2.txt');
y_full=size(coord,1);        % The number of atoms in data array coord
dim=size(coord,2);           % dimention of atom's vector

%if y_num>y_full, error('y_num>yf_full'); end

y_num=min(y_full,y_num);
i=floor(y_full/y_num); j=i*y_num;
% take y_num atoms evenly form coord
y=coord(1:i:j,:)';           % a set of atoms in fixed distribution
% set equal probability
q=ones(1,y_num)/y_num;       % probabilities of atoms in fixed distribution

%-------------------------------------------------------
%                        Run algorithm 
%-------------------------------------------------------
[d_optim,x,x_initial,p] = transport_problem_All_l2_VK(y,q,x_num,iprint,FileName);
d_minimal=d_optim(size(d_optim,2));
p_minimal=p(:,size(p,2));
x_minimal=x;
%-------------------------------------------------------

clear matrix_*

% print total time 
if iprint>0, fid=fopen(FileName,'a');
  i=fprintf(fid, '                   Total time               %9.2f\n',toc(t1));
  fclose(fid);
end

%-------- Print init data and rezults in files ---------
if iprint>0, 
  fid=fopen('Coordinats_of_Atoms_in_Fixed_Distribution.txt','w');
  i=fprintf(fid, '%7i x %2i\n',y_num,dim);
  for j=1:y_num, i=fprintf(fid, '%9.6f  %9.6f\n',y(:,j)); end
  fclose(fid);
  fid=fopen('Initial_solution_X.txt','w');
  i=fprintf(fid, 'Number %5i     Objective %17.10g\nProbabilities  Coordinates\n',x_num,d_optim(1));
  for j=1:x_num, i=fprintf(fid, '%9.6f      %9.6f  %9.6f\n',p(j,1),x_initial(:,j)); end
  fclose(fid);
  fid=fopen('Final_solution_X.txt','w');
  i=fprintf(fid, 'Number %5i     Objective %17.10g\nProbabilities  Coordinates\n',x_num,d_optim(size(d_optim,2)));
  for j=1:x_num, i=fprintf(fid, '%9.6f      %9.6f  %9.6f\n',p(j,size(p,2)),x(:,j)); end
  fclose(fid);
end

%------- Draw picture if dimention = 2 ----------

if size(y,1)==2

  dist_xy = [];
  for i=1:size(x,2)
    dist_xy_loc=[];
    for j=1:size(y,2)   
        dist_xy_loc = [dist_xy_loc,sqrt(sum((x(:,i)-y(:,j)).^2))];   
    end
    dist_xy=[dist_xy;dist_xy_loc];
  end

  pos1=[];
  for j=1:size(y,2)
    pos1=[pos1,find(dist_xy(:,j)==min(dist_xy(:,j)))];    
  end

  h2 = figure;
  title('Plot 1')

  miy1=min(y(1,:)); miy2=min(y(2,:));
  d1=max(y(1,:))-miy1; d2=max(y(2,:))-miy2;
  c1=miy1+d1/2; c2=miy2+d2/2;
  da=0.6*max(d1,d2); 

  axis([c1-da, c1+da, c2-da, c2+da]);
  %axis([c1-d1/2, c1+d1/2, c2-d2/2, c2+d2/2]);
  
  hold on
  MSpec='o';   %'vs^od+*.x'; %'ph'   % points type
  clrs='bmrgyc';    % 'kw'           % points color
  [yj,ii]=sort(x(2,:));
  for j=1:size(x,2); i=ii(j);
    pos2 = find(pos1==i);
    mrk = MSpec(mod(j-1,size(MSpec,2))+1);
    clr = clrs(mod(j-1,size(clrs,2))+1); %rand(1,3);
    loc=[mrk clr];
    scatter(y(1,pos2),y(2,pos2),50,loc,'filled');   % fixed distribution
    scatter(x(1,i),x(2,i),100,'xk','LineWidth',2)   % approximating (variable) distribution
  end
  legend('colored are fixed atoms','X are approximating atoms');

end %if size(y,1)==2

clear x p x_initial d_optim datax_loc i j point_* t1 varx_loc c1 c2 clr clrs d1 d2 da dist_xy dist_xy_loc h2 ii loc miy1 miy2 mrk MSpec pos1 pos2 yj