
%% reproduce a superset of Fig 5. in "On the equivalence between a minimal 
% codomain cardinality Riesz basis construction, a system of Hadamard-Sylvester 
% operators, and a class of sparse, binary optimisation problems", 
% Nelson, IEEE TSP

%% preamble
clear all
close all
clc
home 

%% check existence of format_ticks.m
fprintf('Checking existence of format_ticks.m... ')
% --- true if format_ticks cannot be found:
if exist('format_ticks', 'file')~=2 
  fprintf('\n')
  emsg1 = 'Cannot find format_ticks.m. Either download it at:';
  emsg2 = ' www.mathworks.com/matlabcentral/';
  emsg3 = 'fileexchange/downloads/80726/download as advised by the read_me';
  emsg4 = ' or hack all the plotting functions yourself.';
  error([emsg1, emsg2, emsg3, emsg4])
else 
    fprintf('done!\n')
end
    
%% set up file output directories
fprintf('This is %s\n', mfilename)
dname = 'figures//ord_0';
CheckCreateDirs({dname}, 1);

fn.enorm0  = strcat(dname, '//enorm_0_%.1d_%.1d');
fn.ord0    = strcat(dname, '//ord_0_%.1d_%.1d');
fn.approx0 = strcat(dname, '//approx_0_%.1d_%.1d');

%% constant parameters
q.ord = 0;
% --- number of time domain samples 
q.N   = 420;
% --- frequency sampling rate 
q.pad = 1;
% --- "upsample" time plots by this factor
pad_t = 8; 

%% parameters used in main loops
% --- frequency baseband to compute (Hz) [try 3,4,5,6,7, etc] 
w_bb = [4:2:22];
% --- weighting function = freq^(-decay) 
weight_decay = [0,1,2,4];

%% gimmicky ascii waitbar 
ii = length(w_bb)*length(weight_decay);
fprintf('Constructing/plotting: \n[')
fprintf(repmat('-', 1, ii));
fprintf(']');
i = 0;

%% loop over params and construct optimal functions
for m = weight_decay
    for n = w_bb        
%--- set up output file names
        fn.enorm  = sprintf(fn.enorm0, m, n);
        fn.approx = sprintf(fn.approx0, m, n);
        fn.ord    = sprintf(fn.ord0, m, n);
        
% === construct helper params [2nd arg = q.w_bb is freq. baseband] ===
        q = getPar(q, n);
        
% === error norms ===
% --- exh. search; L1/2 norms ("correct" global minimum error at t = 1/6) 
        enorm_0 = getEnorm0(q, m); %[2nd arg = weight_decay]
        
% --- plot error norms 
        plotEnorm0L1L2(q, enorm_0)
        print('-depsc', fn.enorm)

% --- optimise (find minimum) 
        [~, argmin_L1] = min(enorm_0.L1(:));
        [~, argmin_L2] = min(enorm_0.L2(:));            
% === reconstruct, and plot, best estimate ===
% --- for L1 and L2
        ord_0_L1 = reconstructOrd0(q, argmin_L1, pad_t);
        ord_0_L2 = reconstructOrd0(q, argmin_L2, pad_t);

        plotOrd0L1L2(q, ord_0_L1, ord_0_L2, pad_t)
        print('-depsc', fn.ord)
        
% === plot "approx" ===
% --- (in freq. domain) [i.e. plot of amplitude spec of derivative] 
        plotApprox0L1L2(q, argmin_L1, argmin_L2, 21)
        print('-depsc', fn.approx)

% --- update gimmicy ascii waitbar        
        i = updateBar(i,ii);
    end
end
fprintf('\ndone!\n')

%%




%%















