% gbp_demo.m
%
% Script file that runs GBP, an Interior Point method, and the Simplex
% method on NPROBLEMS randomly generated problems in D dimensions using
% a K times overcomplete dictionary. Outputs the resulting errors,
% L1-norms, and running times.
%
% Assumes that you are running Matlab 7.
%
fprintf('\nStarting GBP demo...');

nproblems = 8; % the number of problems to run
d = 64; % dimension of the signal
k = 4; % overcompleteness of the dictionary
nAtoms = k*d; % number of atoms in the dictionary

% variables to store running times, errors, and L1-norms
% for each method tested
time_gbp = zeros(nproblems, 1);
time_ip = zeros(nproblems, 1);
time_simplex = zeros(nproblems, 1);

error_gbp = zeros(nproblems, 1);
error_ip = zeros(nproblems, 1);
error_simplex = zeros(nproblems, 1);

l1norm_gbp = inf(nproblems, 1);
l1norm_ip = inf(nproblems, 1);
l1norm_simplex = inf(nproblems, 1);

for iproblem = 1:nproblems,
    fprintf('\n----------------------------------------------------------------------\n');
    fprintf('Problem %i (of %i). Dimension: %i, Atoms: %i (%i-times overcomplete).\n', iproblem, nproblems, d, nAtoms, k);
    %---------------------------------------------
    % Initialize
    clear D;
    % generate a random problem
    D = zeros(2*nAtoms, d);
    [D_pos, x] = random_bp_problem(d, nAtoms);
    % "double" the dictionary to ensure that coefficients are positive
    D_neg = -D_pos;
    % perturb atoms for sake of LP solvers
    % (otherwise matrix inverse can fail)
    epsilon = 0.000001;
    for i = 1 : nAtoms,
        v = rand(1,d) - 0.5;
        v = v / norm(v);
        D_pos(i,:) = D_pos(i,:) + epsilon*v;
        D_pos(i,:) = D_pos(i,:) / norm(D_pos(i,:));
        v = rand(1,d) - 0.5;
        v = v / norm(v);
        D_neg(i,:) = D_neg(i,:) + epsilon*v;
        D_neg(i,:) = D_neg(i,:) / norm(D_neg(i,:));
    end
    D(1:nAtoms,:) = D_pos;
    D(nAtoms+1:2*nAtoms, :) = D_neg;
    %---------------------------------------------
    % Greedy Basis Pursuit
    fprintf('\nRunning Greedy Basis Pursuit...\n');
    % Run
    time = cputime;
    [i_list_gbp, alpha_list_gbp] = gbp_safe(D, x); % safe version
    % [i_list_gbp, alpha_list_gbp] = gbp(D, x); % raw version
    time = cputime - time;
    % Record results
    time_gbp(iproblem) = time;
    error_gbp(iproblem) = sqrt(sum((x - alpha_list_gbp*D(i_list_gbp,:)).^2));
    l1norm_gbp(iproblem) = sum(alpha_list_gbp);
    %---------------------------------------------
    % Basis Pursuit via Linear Programming
    % Identifications for Matlab:
    %  x = linprog(f,A,b,Aeq,beq,lb,ub,x0,options)
    %  x <-- alpha_list
    %  f <-- ones(n,1)
    %  Aeq <-- D'
    %  beq <-- x'
    %---------------------------------------------
    % Linear Programming -- Interior Point Method (Matlab 7)
    fprintf('\nRunning Interiot Point method (Matlab 7)...\n');
    % Set options of Interior Point method (LIPSOL)
    options = optimset('LargeScale','on');
    % Run
    time = cputime;
    alpha_list_ip = linprog(ones(2*nAtoms,1),[],[],D',x',zeros(2*nAtoms,1),inf*ones(2*nAtoms,1),[],options);
    time = cputime - time;
    % Record results
    time_ip(iproblem) = time;
    error_ip(iproblem) = sqrt(sum((x - alpha_list_ip'*D).^2));
    l1norm_ip(iproblem) = sum(alpha_list_ip);
    %---------------------------------------------
    % Linear Programming -- Simplex Method (Matlab 7)
    fprintf('\nRunning Simplex method (Matlab 7)...\n');
    % Set options for Simplex method
    options = optimset('LargeScale','off','Simplex','on');
    % Run
    time = cputime;
    alpha_list_simplex = linprog(ones(2*nAtoms,1),[],[],D',x',zeros(2*nAtoms,1),inf*ones(2*nAtoms,1),[],options);
    time = cputime - time;
    % Record results
    time_simplex(iproblem) = time;
    error_simplex(iproblem) = sqrt(sum((x - alpha_list_simplex'*D).^2));
    l1norm_simplex(iproblem) = sum(alpha_list_simplex);    
    %---------------------------------------------
    % print results to screen
    fprintf('\n*Results*\n\n');
    fprintf(' Simplex:\n');
    fprintf('     Error: %f\n', error_simplex(iproblem));
    fprintf('     L1-norm: %f\n', l1norm_simplex(iproblem));
    fprintf('     Time: %f\n', time_simplex(iproblem));
    fprintf(' Interior Point:\n');
    fprintf('     Error: %f\n', error_ip(iproblem));
    fprintf('     L1-norm: %f\n', l1norm_ip(iproblem));
    fprintf('     Time: %f\n', time_ip(iproblem));
    fprintf(' GBP:\n');
    fprintf('     Error: %f\n', error_gbp(iproblem));
    fprintf('     L1-norm: %f\n', l1norm_gbp(iproblem));
    fprintf('     Time: %f\n', time_gbp(iproblem));
    %---------------------------------------------
end
%---------------------------------------------
% print summary to screen
fprintf('\n----------------------------------------------------------------------\n');
fprintf('*Summary*\n\n');
fprintf('%i problems. Dimension: %i, Atoms: %i (%i-times overcomplete).\n', nproblems, d, nAtoms, k);
fprintf('\n __________Error_____________________\n')
fprintf('                   min |  max | mean \n');
fprintf('         Simplex: %3.2f | %3.2f | %3.2f \n', min(error_simplex), max(error_simplex), mean(error_simplex));
fprintf('  Interior Point: %3.2f | %3.2f | %3.2f \n', min(error_ip), max(error_ip), mean(error_ip));
fprintf('             GBP: %3.2f | %3.2f | %3.2f \n', min(error_gbp), max(error_gbp), mean(error_gbp));
fprintf('\n ________L1-norm_______________________________________\n');
fprintf('                         min |        max |       mean \n');
fprintf('         Simplex: %10.4f | %10.4f | %10.4f \n', min(l1norm_simplex), max(l1norm_simplex), mean(l1norm_simplex));
fprintf('  Interior Point: %10.4f | %10.4f | %10.4f \n', min(l1norm_ip), max(l1norm_ip), mean(l1norm_ip));
fprintf('             GBP: %10.4f | %10.4f | %10.4f \n', min(l1norm_gbp), max(l1norm_gbp), mean(l1norm_gbp));
fprintf('\n ___Running time_______________________________________\n');
fprintf('                         min |        max |       mean \n');
fprintf('         Simplex: %10.4f | %10.4f | %10.4f \n', min(time_simplex), max(time_simplex), mean(time_simplex));
fprintf('  Interior Point: %10.4f | %10.4f | %10.4f \n', min(time_ip), max(time_ip), mean(time_ip));
fprintf('             GBP: %10.4f | %10.4f | %10.4f \n', min(time_gbp), max(time_gbp), mean(time_gbp));
fprintf('\n----------------------------------------------------------------------\n');
fprintf('...GBP demo complete.\n');