%--------------------------------------------------------------------------
% LP_Solver_1d.m
% David Shirokoff, NJIT Dept. Mathematical Sciences
% Sept 2016
%
%
% Uses MATLAB's LP solver to compute the solution to the 
% pimal-dual problems (R) and (D). 
%
% Variables to change:
% W(x) --- input function. This should be symmetric (W(-x) = W(x)).
% n --- number of gridpoints.
%
%--------------------------------------------------------------------------
function LP_Solver_1d

    n = 300;                                    %number gridpoints <----- Change this number to increase numerical resolution.
    L = 1;                                      %domain length
    x = (0:n-1)*L/n;                            %grid mesh
    dx = L/n;                                   %grid spacing
    
                                                %Function to be decomposed:
    Wx = exp(cos(10*pi*x)) - exp(cos(2*pi*x));  %<-----change this for different function decompositions. 
       
                                                %Function from the paper [Bandegi, Shirokoff, SIADS 2018]
    plot_SIADS_function = 0;                    %<---- set = 1, to plot the function from the SIADS paper.
    if (plot_SIADS_function == 1)
        alpha = 1/10;
        B1 = (1/(1 - exp(-1/(L*alpha)) ) )*(exp(-x/(L*alpha)) + exp(-(1-x)/(L*alpha) ) );
        B2 = (1/(1 - exp(-1/(alpha))))*(exp(-x/(alpha)) + exp(-(1-x)/(alpha)));
        Wx = -G*L*B1 + B2;
    end
    
                                                %Change this to improve error tolerances in MATLABs linprog solvers.
    options = optimoptions(@linprog, 'TolFun', 10^-6, 'MaxIter',500, 'Algorithm', 'interior-point')
    
                                                % Build and solve the linear programmin problem R.
    [fr,fval, Kx, Wplus, energy_R, Wx] =  SolveLP_Problem_1D(n, Wx, x, dx, options);
    
    
    %----------------------------------------------------------------------
    % Plot results:
    %
    LW = 'linewidth'; lw = 3; MS = 'markersize'; ms = 8;
    
    figure;
    plot(x, fr, 'Color', 'k', LW, lw, MS, ms);
    title(strcat('F(x), grid size n=', num2str(n)));
    set(gca,'FontSize',30)
    grid on, box on, axis square
    
    figure
    plot(x, Wx, 'Color', 'k', LW, lw, MS, ms);
    title(strcat('W(x), grid size n=', num2str(n)));
    set(gca,'FontSize',30)
    axis square
    grid on, box on, axis square
    
    figure
    plot(x, Wplus, 'Color', 'k', LW, lw, MS, ms);
    title(strcat('W^+(x), grid size n=', num2str(n)));
    set(gca,'FontSize',30)
    axis square
    grid on, box on, axis square
    
    
    figure
    plot(x, Kx, 'Color', 'k', LW, lw, MS, ms);
    title(strcat('K(x), grid size n=', num2str(n)));
    set(gca,'FontSize',30)
    axis square
    grid on, box on, axis square
    %----------------------------------------------------------------------
    
    
end

%--------------------------------------------------------------------------
% Function: SolveLP_Problem
% Description: Builds and solves LP problem in 1D
% Inputs: 
% n -- grid size
% Wx -- interaction potential 
% x -- grid
% dx -- grid space
% option -- LP tolerance object.
%
% Outputs: 
% fr -- solution to (R), optimal autocorrelation.
% fval -- optimal objective function value
% Kx -- output function, non-negative Fourier modes.
% Wplus -- output function, non-negative function. 
% opt_E -- output constant
% Together: Wx = Wplus + Kx + opt_E
% (opt_E agrees with fval)
%
function [fr,fval, Kx, Wplus, opt_E, Wx] = SolveLP_Problem_1D(n, Wx, x, dx, options)
    
    %-----------------------------------------------------------
    % Build constraint and equality matrices
    % Constraints are of the form:
    % A_ineq f <= 0.
    % where A_ineq = [-A_Cos];
    % A_eq = [Mass constraint]
    % b_eq = [ 1 ]
    % lb = [zeros], Use Matlabs built in upper and lower bound
    % on primal variables.
    % We build the new objective function by using the fact that W(x) is 
    % even.
    % The objective function:
    % <F, W> = Sum_{j=0}^N-1 F_j W_j = F_0 W_0 + 2 F_1 W_1 + 2 F_2 W_2 + ... 
    % + 2 F_nhalf-1 W_nhalf-1 + F_nhalf W_nhalf.
    % We therefore define W_pot_half = [W_0 2W_1 2W_2 ... 2W_nhalf-1 W_nhalf]
    % as the new linear vector in the objective function.
    %
    
    nhalf = floor(n/2) + 1;    %We construct even symmetric solutions that have length nhalf        
    W_pot_half = Wx(1: nhalf);
    W_pot_half(2:end-1) = W_pot_half(2:end-1)*2;
    
    %----------------------------------------------------------------------
    % Inequality matrix:
    % i.e., build cosine matrix using fft.
    %
    A_ineq = zeros(nhalf-1, n);
    for kk = 2:nhalf
        v = zeros(n,1);
        v(kk) = 1;
        A_ineq(kk-1, :) = -real(fft(v));        
    end
    A_ineq = A_ineq(1:end, 1:nhalf);
    A_ineq(:, 2:end-1) = A_ineq(:, 2:end-1)*2;
    b_ineq = zeros(1, nhalf-1)';
    
    %----------------------------------------------------------------------
    % Equality matrix:
    %   
    A_eq = [ones(1, nhalf)];
    A_eq(2:end-1) = A_eq(2:end-1)*2;
    b_eq = [n];
    
    %----------------------------------------------------------------------
    % Non-negativity constraint. It is computationally faster to use 
    % MATLABs 'lower bound' constraint.  However, MATLAB then also requires
    % an upper bound constraint.  The upper bound constraint is redundant 
    % and so we set it equalt to 10*nhalf.  Note that for 'fr' to be a 
    % probability distribution, the max 'fr' can be is a discrete Delta with
    % magnitude ~1/n.
    %
    ub = (10*nhalf+1)*ones(1, nhalf)';      %Build in the upper and lower bound constraints explicitly.
    lb = zeros(1, nhalf)';
    init_vect = ones(1,nhalf)'/nhalf;       %Active set method requires a feasible point as initialization.
    
    [f_best,fval,exitflag,output,lambda] = linprog(W_pot_half, A_ineq, b_ineq, A_eq, b_eq, lb, ub, init_vect, options);

    
    %----------------------------------------------------------------------
    % Obtain the dual variables from lambda for plotting:
    % 
    w_plus_temp = lambda.lower/2;
    w_plus_temp(end) = 2*w_plus_temp(end);
    
    Wplus = zeros(n, 1);
    Wplus(1:nhalf) = w_plus_temp;
    Wplus(nhalf+1:end) = fliplr(w_plus_temp(2:end-1)')';
    
    Kk = zeros(n, 1);
    Kk(2:nhalf) = lambda.ineqlin(:);
    Kk(nhalf+1:end) = fliplr(lambda.ineqlin(1:end-1)');
    Kx = ifft(Kk)*n/2;
    
    %----------------------------------------------------------------------
    % Reconstruct the solution fr on the full domain by extending the 
    % solution f_best to be even. 
    fr = zeros(n, 1);
    fr(1:nhalf) = f_best;
    fr(nhalf+1:end) = fliplr(f_best(2:end-1)')';    
    opt_E = fval*dx;    
end




