s2let_demo_ridgelet_evaluate

 s2let_demo_ridgelet_evaluate
 Evaluate timing and error of ridgelet transform.

 S2LET package to perform Wavelets on the Sphere.
 Copyright (C) 2012-2015  Boris Leistedt & Jason McEwen
 See LICENSE.txt for license details

0001 
0002 % s2let_demo_ridgelet_evaluate
0003 % Evaluate timing and error of ridgelet transform.
0004 %
0005 % S2LET package to perform Wavelets on the Sphere.
0006 % Copyright (C) 2012-2015  Boris Leistedt & Jason McEwen
0007 % See LICENSE.txt for license details
0008 
0009 clear;
0010 
0011 N_test = 3
0012 
0013 B = 2;
0014 J_min = 0;
0015 spin = 2;
0016 reality = false;
0017 sampling_method = 'MWSS';
0018 save_plots = false;
0019 
0020 Ls = [32 64] % 128 256 512]
0021 
0022 err = zeros(N_test, length(Ls));
0023 time_analysis = zeros(N_test, length(Ls));
0024 time_synthesis = zeros(N_test, length(Ls));
0025 
0026 el_ind = 0;
0027 for L = Ls
0028    el_ind = el_ind + 1
0029    
0030    for n = 1:N_test
0031    
0032       n
0033       
0034       % Generate band-limited function in harmonic space
0035       flm = zeros(L^2,1);
0036       flm = rand(size(flm)) + sqrt(-1)*rand(size(flm));
0037       flm = 2.*(flm - (1+sqrt(-1))./2);
0038       
0039       % For spin signals, set coefficients for el < spin to zero.
0040       if spin > 0
0041          for el = 0:abs(spin)-1
0042             for m = -el:el
0043                ind = ssht_elm2ind(el, m);
0044                flm(ind) = 0;
0045             end
0046          end
0047       end
0048       
0049       % Impose condition to ensure invertible
0050       % (imposes antipodal symmetry for scalar signals).
0051       for el = max([0 abs(spin)]):L-1         
0052          if mod(el + spin, 2) == 1
0053             for m = -el:el
0054                ind = ssht_elm2ind(el, m);
0055                flm(ind) = 0;
0056             end
0057          end
0058       end
0059             
0060       % Impose antipodal symmetry for spin signals.
0061 %       for el = 0:L-1
0062 %          ind = el*el + el + 1;
0063 %          if mod(el + spin, 2) == 1
0064 %             flm(ind,1) = 0;
0065 %          end
0066 %          for m = 1:el
0067 %             ind_pm = el*el + el + m + 1;
0068 %             ind_nm = el*el + el - m + 1;
0069 %             flm(ind_nm,1) = (-1).^el .* (-1)^m .* (-1)^spin * conj(flm(ind_pm,1));
0070 %          end
0071 %       end
0072       
0073       % Construct the corresponding signal on the sphere
0074       f = ssht_inverse(flm, L, 'Method', sampling_method, 'Spin', spin);                  
0075       
0076       tstart = tic;
0077       [f_ridgelet_wav, f_ridgelet_scal] = s2let_ridgelet_analysis(f, ...
0078          'B', B, 'J_min', J_min, ...
0079          'Upsample', true, 'Spin', spin, 'Reality', reality, ...
0080          'Sampling', sampling_method);
0081       time_analysis(n, el_ind) = toc(tstart);
0082       
0083       tstart = tic;
0084       f_recov = s2let_ridgelet_synthesis(f_ridgelet_wav, f_ridgelet_scal, ...
0085          'L', L, 'B', B, 'J_min', J_min, ...
0086          'Upsample', true, 'Spin', spin, 'Reality', reality, ...
0087          'Sampling', sampling_method);
0088       time_synthesis(n, el_ind) = toc(tstart);
0089       
0090       err(n, el_ind) = max(abs(f(:) - f_recov(:)));
0091       
0092    end
0093 
0094 end
0095 
0096 err_mean = mean(err);
0097 err_std = std(err);
0098 err_std_log = std(log10(err));
0099 
0100 time_analysis_mean = mean(time_analysis);
0101 time_analysis_std = std(time_analysis);
0102 
0103 
0104 time_synthesis_mean = mean(time_synthesis);
0105 time_synthesis_std = std(time_synthesis);
0106 
0107 
0108 time_total = time_analysis + time_synthesis;
0109 time_total_mean = mean(time_total);
0110 time_total_std = std(time_total);
0111 time_total_std_log = std(log10(time_total));
0112 
0113 
0114 
0115 %% Define plotting parameters.
0116 
0117 istart = 1;
0118 iend = length(Ls);
0119 
0120 line_width = 1.8;
0121 line_width_thick = 2.5;
0122 marker_size = 7;
0123 marker_type = 'o';
0124 green_light = [0.2 0.6 0.4];% x339966
0125 green_dark = [0 0.4 0.2];   % x006633
0126 blue_light = [0.2 0.4 0.8]; % x3366CC
0127 blue_dark = [0 0 1];        % x0000FF
0128 red_light = [1 0.4 0.2];    % xFF6633
0129 red_dark = [0.8 0.2 0];     % xCC3300
0130 
0131 
0132 % Plot error
0133 figure;
0134    
0135 plot(log2(Ls(istart:iend)), ...
0136    log10(err_mean(istart:iend)), ...
0137    'Color', green_dark, ...
0138    'Marker', marker_type, ...
0139    'MarkerSize', marker_size, ...
0140    'MarkerFaceColor', green_light, ...
0141    'MarkerEdgeColor', green_dark, ...
0142    'LineStyle', '--', ...
0143    'LineWidth', line_width);
0144 
0145 hold on;
0146 
0147 plot(log2(Ls(istart:iend)), ...
0148    log10(Ls(istart:iend).^2)-15, ...
0149    'r', ...
0150    'LineWidth', line_width_thick);
0151 a = axis;
0152 set(gca,'XTick',a(1):1:a(2));
0153 set(gca,'XTickLabel', 2.^[a(1):1:a(2)]);
0154 set(gca,'YTick',floor(a(3)):1:ceil(a(4)));
0155 set(gca,'YTickLabel',{10.^[floor(a(3)):ceil(a(4))]});
0156 xlabel('L');
0157 ylabel('E');
0158 
0159 set(gca, 'LineWidth', 3)
0160 set(gca, 'FontSize', 20)
0161 
0162 if save_plots, print('-depsc2', 'plots/ridgelet_error.eps'); end
0163 
0164 % Plot error with error bars
0165 figure
0166   
0167 errorbar(log2(Ls(istart:iend)), ...
0168    log10(err_mean(istart:iend)), ...
0169    err_std_log(istart:iend), ...
0170    'Color', green_dark, ...
0171    'Marker', marker_type, ...
0172    'MarkerSize', marker_size, ...
0173    'MarkerFaceColor', green_light, ...
0174    'MarkerEdgeColor', green_dark, ...
0175    'LineStyle', '--', ...
0176    'LineWidth', line_width);
0177 
0178 hold on;
0179 
0180 plot(log2(Ls(istart:iend)), ...
0181    log10(Ls(istart:iend).^2)-15, ...
0182    'r', ...
0183    'LineWidth', line_width_thick);
0184 a = axis;
0185 set(gca,'XTick',a(1):1:a(2));
0186 set(gca,'XTickLabel', 2.^[a(1):1:a(2)]);
0187 set(gca,'YTick',floor(a(3)):1:ceil(a(4)));
0188 set(gca,'YTickLabel',{10.^[floor(a(3)):ceil(a(4))]});
0189 xlabel('L');
0190 ylabel('E');
0191 
0192 set(gca, 'LineWidth', 3)
0193 set(gca, 'FontSize', 20)
0194 
0195 if save_plots, print('-depsc2', 'plots/ridgelet_error_errorbars.eps'); end
0196 
0197 
0198 % Plot timing
0199 figure;
0200    
0201 plot(log2(Ls(istart:iend)), ...
0202    log10(time_total_mean(istart:iend)), ...
0203    'Color', green_dark, ...
0204    'Marker', marker_type, ...
0205    'MarkerSize', marker_size, ...
0206    'MarkerFaceColor', green_light, ...
0207    'MarkerEdgeColor', green_dark, ...
0208    'LineStyle', '--', ...
0209    'LineWidth', line_width);
0210 
0211 hold on;
0212 
0213 plot(log2(Ls(istart:iend)), ...
0214    log10(Ls(istart:iend).^3)-5, ...
0215    'r', ...
0216    'LineWidth', line_width_thick);
0217 a = axis;
0218 set(gca,'XTick',a(1):1:a(2));
0219 set(gca,'XTickLabel', 2.^[a(1):1:a(2)]);
0220 set(gca,'YTick',floor(a(3)):1:ceil(a(4)));
0221 set(gca,'YTickLabel',{10.^[floor(a(3)):ceil(a(4))]});
0222 xlabel('L');
0223 ylabel('T');
0224 
0225 set(gca, 'LineWidth', 3)
0226 set(gca, 'FontSize', 20)
0227 
0228 if save_plots, print('-depsc2', 'plots/ridgelet_timing.eps'); end
0229 
0230 
0231 % Plot timing with error bars
0232 figure;
0233    
0234 errorbar(log2(Ls(istart:iend)), ...
0235    log10(time_total_mean(istart:iend)), ...
0236    time_total_std_log(istart:iend), ...
0237    'Color', green_dark, ...
0238    'Marker', marker_type, ...
0239    'MarkerSize', marker_size, ...
0240    'MarkerFaceColor', green_light, ...
0241    'MarkerEdgeColor', green_dark, ...
0242    'LineStyle', '--', ...
0243    'LineWidth', line_width);
0244 
0245 hold on;
0246 
0247 plot(log2(Ls(istart:iend)), ...
0248    log10(Ls(istart:iend).^3)-5, ...
0249    'r', ...
0250    'LineWidth', line_width_thick);
0251 a = axis;
0252 set(gca,'XTick',a(1):1:a(2));
0253 set(gca,'XTickLabel', 2.^[a(1):1:a(2)]);
0254 set(gca,'YTick',floor(a(3)):1:ceil(a(4)));
0255 set(gca,'YTickLabel',{10.^[floor(a(3)):ceil(a(4))]});
0256 xlabel('L');
0257 ylabel('T');
0258 
0259 set(gca, 'LineWidth', 3)
0260 set(gca, 'FontSize', 20)
0261 
0262 if save_plots, print('-depsc2', 'plots/ridgelet_timing_errorbars.eps'); end
0263 
0264