% REPEATED_SEARCH_TEST  A test of the efficiency of the FBMP tree search 
% and the Larsson/Selen tree search in recovering all "highly probable" 
% mixture selection vectors.  This program will perform an exhaustive 
% search of all possible mixture selection vectors and determine their 
% posterior likelihood.  It will use summary statistics to capture 
% important performance metrics for both algorithms
%
% * Note that this function was coded to make use of parallel Matlab
% capabilities for faster execution, and was equipped with features that
% enable job recovery if the job terminated unexpectedly before completing,
% as well as a feature that suspends execution during work-day hours.
% Simple modifications can remove these features if not desired/available *
%

%
% Coded by: Justin Ziniel, Ohio State Univ.
% Email: ziniel.1@osu.edu
% Last modified: April 2, 2009
% FBMP version 1.2 
% Copyright (c) Philip Schniter, Lee C. Potter, and Justin Ziniel, 2008
%

clear; clc
randn('state', 0);      % Return random normal generator to default state
rand('twister', 5489);  % Return random generator to default state

%% Declare important constants

datafilename = 'REPEATED_SEARCH_RESULTS.mat';   % Filename for saving results
nice = 0;       % If nice == 1, pause execution during the hours 6am - 6pm

try
    load(datafilename);             % Attempt to load an existing data file
    out_start = floor(min(find(MMSE_error_fbmp == inf))/12);  % Start trial count where left off
    rand('twister', rand_state);    % Set random number generator state
    randn('state', randn_state);    % Set normal random state
    time = datevec(datenum(clock) - datenum([0000, 0, 0, 0, 0, est_time*...
                    60*60*12*out_start/(T - 12*out_start)]));
    disp(sprintf('Resuming search efficiency test at trial %d...', ((12*out_start) + 1)));
catch   % This is a new test, thus initialize variables

    % These parameters can be manipulated by the user to try various setups
    T = 204;          	% # of trials    
    N = 24;   			% # parameters
    M = ceil(N/3);      % # observations
    p1 = 0.04;			% prior probability of non-zero tap
    SNRdB = 15;			% AWGN variance
    
    % Setup the variables that will hold outcomes
    MSE_fbmp_mmse = inf*ones(T,1);     % NMSE of FBMP's MMSE estimator
    MSE_fbmp_map = inf*ones(T,1);      % NMSE of FBMP's MAP estimator
    MSE_ls_mmse = inf*ones(T,1);       % NMSE of Larsson/Selen's MMSE estimator
    MSE_ls_map = inf*ones(T,1);        % NMSE of Larsson/Selen's MAP estimator
    MMSE_error_fbmp = inf*ones(T,1);    % Error between true MMSE estimator and FBMP's
    MMSE_error_ls = inf*ones(T,1);      % Error between true MMSE estimator and Larsson/Selen's
    mass_ratio_fbmp = inf*ones(T,1);    % Ratio of prob. mass returned by FBMP to total
    mass_ratio_ls = inf*ones(T,1);      % Ratio of prob. mass returned by Larsson/Selen to total
    
    % Reset index start values
    out_start = 0;
    if rem(T, 12) ~= 0
        error('# of trials should be evenly divisible by 12')
    else
        out_length = T/12;
    end
    
    
    % Reset random number generators
    rand('twister',5489)
    randn('state',0)
    
    time = clock;
end


%% Execute test in parallel

path(path, '../functions/')

matlabpool open

for out = out_start : out_length - 1
    parfor t = ((12*out) + 1):((12*out) + 12)  
        % Run a trial
        [MSE_fbmp_mmse_tr, MSE_fbmp_map_tr, MSE_ls_mmse_tr, ...
         MSE_ls_map_tr, MMSE_error_fbmp_tr, MMSE_error_ls_tr, ...
         mass_ratio_fbmp_tr, mass_ratio_ls_tr] = search_efficiency_fxn(N, M, ...
                                                                       p1, SNRdB);

        % Store outcome
        MSE_fbmp_mmse(t,1) = MSE_fbmp_mmse_tr;
        MSE_fbmp_map(t,1) = MSE_fbmp_map_tr;
        MSE_ls_mmse(t,1) = MSE_ls_mmse_tr;
        MSE_ls_map(t,1) = MSE_ls_map_tr;
        MMSE_error_fbmp(t,1) = MMSE_error_fbmp_tr;
        MMSE_error_ls(t,1) = MMSE_error_ls_tr;
        mass_ratio_fbmp(t,1) = mass_ratio_fbmp_tr;
        mass_ratio_ls(t,1) = mass_ratio_ls_tr;
    end

    %****************************************************************%
    %                         TRIAL CONTROL                          %
    %****************************************************************%

    % Provide an update of the estimated time to completion
    est_time = (T - (12*(out + 1)))*...
                (etime(clock, time)/(12*(out + 1)))/60/60;
    disp(sprintf(['Estimated time remaining: %3.2f hours. '...
                  ' |  Percentage complete: %2.1f  |  Printed: ' ...
                  datestr(clock, 21)], est_time, ...
                  100*(12*(out + 1))/T));

    % Save progress
    rand_state = rand('twister');
    randn_state = randn('state');
    save(datafilename);
    disp(sprintf('Data autosaved to: %s', datafilename));

    % Check for possible job suspension
    if (nice == 1)  % Suspension-mode active
        walltime = clock;
        while (walltime(4) >= 6) && (walltime(4) < 18)
            % Restricted hours, therefore, suspend execution unless
            % the file 'override.txt' exists in the current
            % directory, in which case, do not suspend
            fid = fopen('override.txt');
            if (fid == -1)  % Override file absent
                pause(60*5);   % Wait 5 minutes before re-check
                time = datevec(datenum(time) + ...
                        datenum([0000, 0, 0, 0, 5, 0]));
                walltime = clock;
            else
                fclose(fid);
                break
            end
        end
    end

    %****************************************************************%
    %                       END TRIAL CONTROL                        %
    %****************************************************************%
end

disp(sprintf('Total runtime: %3.2f hours', etime(clock, time)/60/60));

matlabpool close

% Save results one last time
clear MSE_fbmp_mmse_tr MSE_fbmp_map_tr MSE_ls_mmse_tr MSE_ls_map_tr ...
      MMSE_error_fbmp_tr MMSE_error_ls_tr mass_ratio_fbmp_tr mass_ratio_ls_tr
save(datafilename);


%% Calculate summary statistics

AMEF = 10*log10(mean(MMSE_error_fbmp));     % [dB]
AMEL = 10*log10(mean(MMSE_error_ls));       % [dB]

AMFMAP = 10*log10(mean(MSE_fbmp_map));   	% [dB]
AMFMMSE = 10*log10(mean(MSE_fbmp_mmse)); 	% [dB]
AMLMAP = 10*log10(mean(MSE_ls_map));        % [dB]
AMLMMSE = 10*log10(mean(MSE_ls_mmse));    	% [dB]

MRF = mean(mass_ratio_fbmp);
MRL = mean(mass_ratio_ls);


%% Display summary statistics

clc
disp('Summary statistics for the search efficiency trial')
disp('--------------------------------------------------')
disp(' ')
disp(sprintf('Average MMSE error for FBMP: %2.2f [dB]', AMEF));
disp(sprintf('Average MMSE error for LS: %2.2f [dB]', AMEL));
disp(' ')
disp(sprintf('Average MSE for FBMP''s MAP estimator: %2.2f [dB]', AMFMAP));
disp(sprintf('Average MSE for LS''s MAP estimator: %2.2f [dB]', AMLMAP));
disp(sprintf('Average MSE for FBMP''s MMSE estimator: %2.2f [dB]', AMFMMSE));
disp(sprintf('Average MSE for LS''s MMSE estimator: %2.2f [dB]', AMLMMSE));
disp(' ')
disp(sprintf('Average mass-ratio captured by FBMP: %2.2f', MRF));
disp(sprintf('Average mass-ratio captured by LS: %2.2f', MRL));
