function OUTPUT= OD_computing(INPUT,OPTIONS)
% -------------------------------------------------------------------------
% INPUT PARAMETERS:     see OD_problem class
% OPTIONS PARAMETER:    see OD_options class
% OUTPUT PARAMETER:     see OD_result class
%
% CALLS:        ONB.m  
% DEMO:         OD_DEMO.m
% -------------------------------------------------------------------------
% The computation of an optimal design for a given function space on a
% given point set is attacked as optimization problem for the regularized 
% design energy functional as follows.
%
% Minimize E_eta(w):=-1/N log det G_deg(w) + eta\|pi_K w\|^2 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, and pi_K is the projection onto the kernel of U^t,
% where U is the Vandermonde matrix of the space of functions
% phi_i(x)phi_j(x).
% Optimization is carried out by integrating the gradient flow equation
%    d/dt z = -grad F_eta(z),
% where F_eta(z) := E_eta(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
%--------------------------------------------------------------------------
%% initializing
% local variables from input
V=INPUT.V;
V2=INPUT.V2;
nodes=INPUT.nodes;
weights=INPUT.weights;
OUTPUT.W=weights;
eta=INPUT.eta;
% local variables to be setted
gradflow_iteration=0;
gradient_res=OPTIONS.gradient_res_toll+1;
ind_positive_weights=ones(size(weights));
I=1:length(weights);
total_iterations=0; % counter of solutions of linear systems
[Q,R]=qr(V2); r= rank(R);
Q0=Q(:,r+1:end);
Projection=Q0*Q0';
% dimensions
N=size(INPUT.V,2);
M=size(INPUT.nodes,1);
M0=M;
% DECLARE OPTIMIZER OUTPUT
OUTPUT=OD_result;
% output variables that change size
OUTPUT.orthcheck=[];
OUTPUT.spectrum=[];
OUTPUT.gradient_res_vec=[];
OUTPUT.spectrum=[];
OUTPUT.kkt_res=[];
% PRE-PROCESSING (ONLY IF OPTIONS.PREMULT>0)
% INITIAL GUESS BY MULTIPLICATIVE ALGORITHM
for k=1: OPTIONS.premult
    [~,B,~] = ONB(INPUT.V,weights.^(1/2));
    weights=weights.*B/N;
end
OUTPUT.W=[OUTPUT.W,weights];
% first orthogonalization
[K,B,~] = ONB(INPUT.V,weights.^(1/2));
if isempty(INPUT.tau)
    tau=max((N/max(B))/2,OPTIONS.tau_min);
else
    tau=INPUT.tau;
end
z=weights.^(1/2);
%% algorithm body
% BACKWARD EULER METHOD
while and(gradflow_iteration < OPTIONS.gradflow_max_iteration,...
        gradient_res>OPTIONS.gradient_res_toll)
    
    % INIZIALIZE OR UPDATE LOCAL VARIABLES OF THE CYCLE
    z_old=z;
    z=z_old;
    restart = 1; % flag for restart newton solver with smaller time step
    n_restart =0; % newton restart counter for local time step
    continue_newton=1;
    
    % CYCLE FOR ADAPTIVE TIME STEP
    while restart==1 % cycle for Newton solver restart
        Newton_iteration=0; % counter of Newton It in the BE t-step
        
        % NEWTON SOLVER FOR EACH TIME STEP
        while and(continue_newton,Newton_iteration<OPTIONS.Newton_max_iteration)
            A0=4* diag(z)*((K.^2)./N+2*eta*Projection)*diag(z)+2*diag(1-B/N+2*eta*Projection*(z.^2));
            A=A0+1/tau*eye(length(z));
            b=-(2*diag(z)*(1-B/N+2*eta*Projection*(z.^2))+(z-z_old)./tau);
            % EVENTUALLY COMPUTES SMALL EIGENVALUES OF THE HESSIAN
            if OPTIONS.checkspectrum==1
                OUTPUT.spectrum=[OUTPUT.spectrum,svds(A0,OPTIONS.n_eigs,'smallest')];
            end
            % SOLVE NEWTON LINEAR SYSTEM
            switch OPTIONS.solver
                case 'backslash'
                    step=A\b;
                case 'pcg'
                    step=pcg(A,b,OPTIONS.lin_toll,OPTIONS.lin_maxit);
            end
            Newton_res=A*step-b;
            z=z+step;
            [K,B,V_ONB] = ONB(V,z);
            Newton_iteration=Newton_iteration+1;
            total_iterations=total_iterations+1;
            % CHECK IF NEWTON STOPPING CRITERION IS MET 
            if min(abs(Newton_res)<=OPTIONS.Newton_res_toll*abs(z-z_old))
                continue_newton=0;
            end
        end
        % CHECK IF NEWTON EXITED WITHOUT MEETING STOPPING CRITERION
        if  not(continue_newton)
            restart=0;
            gradflow_iteration=gradflow_iteration+1;
            tau=OPTIONS.alpha*tau;
            ind_zeroweights=(z<=(OPTIONS.kkt_toll).^(1/2));
            ind_posweights=(z>(OPTIONS.kkt_toll).^(1/2));
            gradE=1-B/N+2*eta*Projection*(z.^2);
        else 
            % RESTART WITH SMALLER TIME STEP
            if n_restart< OPTIONS.n_restart
                tau=tau/OPTIONS.beta;
                n_restart= n_restart+1;
                fprintf('Newton solver failed\n')
                fprintf('restart number %d, time step reduction: tau =%g\n',[n_restart,tau])
            else
                fprintf('maximum number %d of restarts reached\n', n_restart)
                fprintf('Backward Euler wrong update\n');
                gradflow_iteration=gradflow_iteration+1;
                break
            end
        end
        fprintf('GF it %d\t tot it %d\t tau %g\t residual %g\n',...
        [gradflow_iteration,total_iterations,tau,gradient_res])
        
    end
    % CUTOFF OF VERY SMALL WEIGHTS (ONLY IF OPTIONS.CUTOFF = 1)
    if OPTIONS.cutoff==1
        ind_positive_weights=(z>(OPTIONS.cutoff_toll/length(z)).^(1/2));
        I=I(ind_positive_weights);
        OUTPUT.W=[OUTPUT.W,zeros(size(OUTPUT.W,1),1)];
        M=sum(ind_positive_weights);
        fprintf('weights cutoff to %d\n',M);
        nodes=nodes(ind_positive_weights,:);
        z=z(ind_positive_weights);
        z=z/norm(z);
        V=V(ind_positive_weights,:);
        [K,B,V_ONB] = ONB(V,z);
        OUTPUT.W(I,end)=z.^2;
        OUTPUT.kkt_res=[OUTPUT.kkt_res,...
                max(abs(min(min(gradE(ind_zeroweights)),0)),...
                max(abs(gradE(ind_posweights))))];
    else
        OUTPUT.W=[OUTPUT.W,z.^2];
        OUTPUT.kkt_res=[OUTPUT.kkt_res,...
                max(abs(min(min(gradE(ind_zeroweights)),0)),...
                max(abs(gradE(ind_posweights))))];
    end
    if or(total_iterations>=OPTIONS.total_max_iterations, tau< OPTIONS.tau_min)
        fprintf('global optimization failed: too many Newton iterations\n')
        break
    end
    gradient_res=norm(2*z.*(B/N-1-2*eta*Projection*(z.^2)));
    OUTPUT.gradient_res_vec=[OUTPUT.gradient_res_vec,gradient_res];
end
OUTPUT.V_ONB=V_ONB;
% CUTOFF OF VERY SMALL WEIGHTS AFTER TIME STEPPING (ONLY IF OPTIONS CUTOFF=1)
if OPTIONS.final_cutoff==1
    ind_positive_weights=(abs(z)>(OPTIONS.final_cutoff_toll/length(z)).^(1/2));
    I=I(ind_positive_weights);
    OUTPUT.W=[OUTPUT.W,zeros(size(OUTPUT.W,1),1)];
    M=sum(ind_positive_weights);
    fprintf('weights cutoff to %d\n',M);
    nodes=nodes(ind_positive_weights,:);
    z=z(ind_positive_weights);
    z=z/sum(z);
    OUTPUT.V_ONB=V_ONB(ind_positive_weights,:);
    OUTPUT.W(I,end)=z.^2;
end
OUTPUT.weights=z.^2;
OUTPUT.nodes=nodes;
OUTPUT.B=B;
OUTPUT.N=N;