% -------------------------------------------------------------------------
% -----------COMPUTATION OF OPTIMAL EXPERIMANTAL DESIGN BY-----------------
% -------NEGATIVE LOG-DETERMINANT OPTIMIZATION BY GRADIENT FLOW------------
% -------------------------------------------------------------------------
% The computation of an optimal design for a given function space on a
% given point set is attacked as optimization problem for the design
% energy functional as follows.
%
% Minimize E(w):=-1/N log det G_deg(w) on w>=0, sum(w)=1, 
% where G_deg is the Gram matrix of degree deg matrix w.r.t. fixed points 
% and with weights w
% Optimization is carried out by integrating the gradient flow equation
%    d/dt z = -grad F(z),
% where F(z) := E(z.^2).
%
% The time integration is approximated by Backward Euler scheme and
% Newton's Method with a stopping criterion suitable for the problem.
% Global convergence of the algorithm has been proved [1] for small values 
% of Newton_res_toll and alpha (see below).
% 
% [1] F. Piazzon "Computing Optimal Experimental Designs on Finite Sets by
% Log-determinant Gradient Flow" 2021 ArXiv preprint
%--------------------------------------------------------------------------
%%
clc
close all
clear all
warning off

%% VARIABLES INITIALIZATION------------------------------------------------

deg=10;
wam_deg_factor=4; 
wam_deg=floor(deg*wam_deg_factor);  
[weights,wam,V]=BuildMesh(deg,wam_deg,'square');

%% SET OPTIMIZER INPUT----------------------------------------------------- 

INPUT=OD_problem;
INPUT.tau=[];
INPUT.nodes=wam;
INPUT.weights=weights;
INPUT.V=V;

%% SET OPTIMIZER OPTIONS--------------------------------------------------- 

OPTIONS=OD_options;

% the class OD_options has default values for options parameters,
% uncomment and edit lines below for changing default options

 OPTIONS.alpha=1.1;
 OPTIONS.beta=1.5;
OPTIONS.gradient_res_toll=10^-13 ;          
OPTIONS.gradflow_max_iteration=200;         
OPTIONS.Newton_max_iteration=5;  
OPTIONS.Newton_res_toll=10^-3;
% OPTIONS.solver='pcg';           
% OPTIONS.lin_toll=10^-8;               
% OPTIONS.lin_maxit=200;
OPTIONS.n_restart=3;
 OPTIONS.total_max_iterations=400;           
% OPTIONS.tau_min=10^-6;    
  OPTIONS.premult=10;                         
% OPTIONS.cutoff=1;                          
% OPTIONS.cutoff_toll=10^-4;                  
% OPTIONS.final_cutoff=1;
% OPTIONS.final_cutoff_toll=10^-10;
% OPTIONS.checkspectrum=1;
% OPTIONS.n_eigs=2;

%% CALL THE OPTIMIZER------------------------------------------------------
tic;
OUTPUT= OD_computing(INPUT,OPTIONS);
time=toc;
N=size(OUTPUT.V_ONB,2);


%% VIDEO OUTPUT
indposweights=OUTPUT.W(:,end)>1e-12;
indzeroweights=OUTPUT.W(:,end)<=1e-12;
grad=1-OUTPUT.B/OUTPUT.N;
gradE_wpos=grad(indposweights);
gradE_wzero=grad(indzeroweights);
fprintf('Function space dimension %d\n', N)
fprintf('Design space dimension %d\n',size(INPUT.nodes,1))
fprintf('Optimal design support cardinality %d\n',size(OUTPUT.nodes,1))
fprintf('KKT conditions satisfied up to %g\n',...
    max(max(gradE_wpos),abs(min(0,min(gradE_wzero)))));
fprintf('Non-linear residual on objective function %g\n',OUTPUT.gradient_res_vec(end))
fprintf('Computation time %f\n',time)

%% PLOTS-------------------------------------------------------------------
figure(1);
semilogy(OUTPUT.W');
title('weights vs iteration')
%%
s=(sum(diff(OUTPUT.W').^2,2)).^(1/2);
figure(2);
subplot(1,2,1);semilogy(s,'k')
title('steps')
subplot(1,2,2);semilogy(OUTPUT.gradient_res_vec,'k--')
hold on
semilogy(OUTPUT.kkt_res,'k-.')
title('residuals')
legend('|grad F|','KKT residual')
%%
figure(3);
plot(INPUT.nodes(:,1),INPUT.nodes(:,2),'.')
hold on
plot(OUTPUT.nodes(:,1),OUTPUT.nodes(:,2),'o')
title('design support')


if OPTIONS.checkspectrum==1
    S=OUTPUT.spectrum';
    figure(6);
    plot(S)
end


