%make_lpsc_2014_width_figures

load('parseTransect1ChannelWidthMeasurementsSecondPassForLPSC2014_v2_output_terse','s','cqual','cpstyle','cx','cy','cz','minx','miny')

scatter(cx(s(10).cwrange>5),cy(s(10).cwrange>5),1500,'*')

nboots = 1000; %Number of bootstraps

okq = 4; %threshhold for 'good quality' data (will be displayed seperately).

%PLANAR CORRECTED STRAT ELEVATION
%Find planar elevation (from f1_f2_contact_second_pass_b20g02_only_plane_fit_one_point_per_contact_pick_24feb14_rms_output ...)    
%of minx, miny

cetl(1) = 6026.681 + minx*(0.00648907) + miny*(-0.0051384);%contact elevation at top left (of transect)

%QUADRATIC CORRECTED STRAT ELEVATION
% Quadratic: coef #          coef
% ------ ----------------
cq(1) = 370424.00040486;
cq(2) = 0.428152656880215;
cq(3) = 0.20357483081356;
cq(4) = 1.33106752829426e-007;
cq(5) =  2.40680350734772e-008;
cq(6) = 2.39768092116689e-007;
% ------ ----------------
% RMS Error  = 17.014793661717
% Chi-Square  = 5500.56086366536
% (from f1_f2_contact_second_pass_b20g02_only_quadratic_fit_one_point_per_contact_pick_24feb14_rms_output.TXT)
% 
% per http://resources.arcgis.com/en/help/main/10.1/index.html#//009z0000007n000000 ('How Trend Works') this is 
% z = c0 + x c1 + y c2 + x^2 c3 + x*y (c4) + y^2 (c5)

%Kludge fix for mixed preservation styles
cpstyle(cpstyle==9) = 8;
cpstyle(cpstyle==10) = 8;

zstrat_method_list = {'planar';'quadratic';'IDW'}; %Methods of interpolating structure contours between picked contacts.

for zstrat_method_index = 2
    zstrat_method = zstrat_method_list(zstrat_method_index);
    
    disp(zstrat_method)

    clear pksnominal pksnominal_good
%********* POSTPROCESSING *******
if strcmp(zstrat_method,'planar') == 1
        czstrat = cz - cetl(1) - cx.*0.00648907 + cy.*0.0051384; %Adjust for planar elevation
        czstrat_planar = czstrat;
elseif strcmp(zstrat_method,'quadratic') == 1
        cx = cx + minx; cy = cy + miny;
        czstrat = cz - (cq(1) + cx.*cq(2) + cy.*cq(3) + (cx.^2).*cq(4) + (cx.*cy).*cq(5) + (cy.^2).*cq(6));
        cx = cx - minx; cy = cy - miny;
        czstrat_quadratic = czstrat;
elseif strcmp(zstrat_method,'IDW') == 1
        WhichPass = 3
        load_channel_widths
        %zstrat obtained here is point-by-point. Need to convert to channel
        %median:
        clear czstrat
        for nn = 1:(max(bid)/2)
                czstrat(nn)    = nanmedian(zstrat((bid==(nn*2) | bid==(nn*2)-1) ));
        end
        czstrat_IDW = czstrat;
        clear zstrat
else
end

    if strcmp(zstrat_method,'planar') == 1
            sigmax = 29.084761000706; %RMS error from f1_f2_contact_second_pass_b20g02_only_plane_fit_one_point_per_contact_pick_24feb14_rms_output
    elseif strcmp(zstrat_method,'quadratic') == 1
            sigmax = 17.01479366; %RMS error
    elseif strcmp(zstrat_method,'IDW') == 1
            sigmax = 27.7; %cross-validation
    else
        error
    end

%Improvements desired: - bootstrap in zstrat as well as height <-- Done!
%                      - consider non-linear fits to data
%                      - consider different structure contour
%                        interpolations (not just planar). <-- Done!

%Preprocessing part 1. Bootstrap linear regression of channel width on stratigraphic elevation
%(redraw with replacement)

%Bootstrap all data 
b_ols = regress(s(11).cw',[ones(size(czstrat')) czstrat']) %bootstrap ordinary least squares
line([min(czstrat),max(czstrat)], ...
     [b_ols(1)+b_ols(2).*min(czstrat),b_ols(1)+b_ols(2).*max(czstrat)]) 
hold on
b = robustfit(czstrat',s(11).cw') %boostrap robust regression (w/ MATLAB defaults)
line([min(czstrat),max(czstrat)], ...
     [b(1)+b(2).*min(czstrat),b(1)+b(2).*max(czstrat)],'Color','r')

for l = 1:nboots
    bootindices = round(0.5+rand(size(czstrat))*length(czstrat)); %resample with replacement
    
    bbs_ols(l,:) = regress(s(11).cw(bootindices)',...
                   [ones(size(czstrat(bootindices)')) czstrat(bootindices)']);
    random_zs_offset_this_bootstrap(l,:) =      sigmax.*randn(length(czstrat),1);
    random_cw_offset_this_bootstrap(l,:) = s(11).cwdev'.*randn(length(czstrat),1);
    bbs_no_zs_errors(l,:) = robustfit(czstrat(bootindices)',s(11).cw(bootindices)');
    
    bootstrapped_proxy_data_values = s(11).cw(bootindices)'+random_cw_offset_this_bootstrap(l,bootindices)';
    bootstrapped_proxy_data_values(bootstrapped_proxy_data_values<0.1) = 0.1; %No negative channel widths!
    
    bbs(l,:) = robustfit(czstrat(bootindices)'+random_zs_offset_this_bootstrap(l,bootindices)', ...
                         bootstrapped_proxy_data_values);
   % hold on
   % line([min(czstrat),max(czstrat)], ...
   %  [bbs(l,1)+bbs(l,2).*min(czstrat),bbs(l,1)+bbs(l,2).*max(czstrat)],'Color','r')
   % line([min(czstrat),max(czstrat)], ...
   %  [bbs_ols(l,1)+bbs_ols(l,2).*min(czstrat), ...
   %   bbs_ols(l,1)+bbs_ols(l,2).*max(czstrat)],'Color','k')
    wbsl = s(11).cw(bootindices); %width bootstrap, this l
    
    wunitbs(l,1,1) =    nanmean(wbsl(czstrat(bootindices) <0)); %Mean channel width above and below paleohydrologic break
    wunitbs(l,2,1) =    nanmean(wbsl(czstrat(bootindices)>=0)); %Mean channel width above and below paleohydrologic break
    wunitbs(l,1,2) =  nanmedian(wbsl(czstrat(bootindices) <0)); %Median channel width above and below paleohydrologic break
    wunitbs(l,2,2) =  nanmedian(wbsl(czstrat(bootindices)>=0)); %Median channel width above and below paleohydrologic break
    wunitbs(l,1,3) =     nanstd(wbsl(czstrat(bootindices) <0)); %s.d. of channel width above and below paleohydrologic break
    wunitbs(l,2,3) =     nanstd(wbsl(czstrat(bootindices)>=0)); %s.d. of channel width above and below paleohydrologic break


end

%Find 1 sigma and 2 sigma errors on robust linear trend
    zszp = 300;%zstratzeropoint
    bss = NaN.*zeros(length([-zszp:round(max(czstrat)-1)]),nboots);
for zstrat = ceil(min(czstrat)):floor(max(czstrat))
    sorted_predictors = sort(bbs(:,1)+bbs(:,2).*zstrat);
    bss(zstrat+zszp,:) = sorted_predictors; %bootstrap output sorted - 5th and 95th percentiles can then be reported as errors.
    sorted_predictors_no_zs_errors = sort(bbs_no_zs_errors(:,1)+bbs_no_zs_errors(:,2).*zstrat);
    bss_no_zs_errors(zstrat+zszp,:) = sorted_predictors_no_zs_errors;
end

%Preprocessing part 1b. Bootstrap only the data with pqual = 2('good') or pqual = 1('gold standard')

wgood = s(11).cw(cqual<=okq);
zgood = czstrat(cqual<=okq);

for l = 1:nboots
    bootindices = round(0.5+rand(size(zgood))*length(zgood)); %resample with replacement
    
    bbs_ols_good(l,:) = regress(wgood(bootindices)',...
                   [ones(size(zgood(bootindices)')) zgood(bootindices)']);

    bbs_good_no_zs_errors(l,:) = robustfit(zgood(bootindices)',wgood(bootindices)');
    %Reuse random_zs_offset_this_bootsrap
    bbs_good(l,:) = robustfit(zgood(bootindices)'+random_zs_offset_this_bootstrap(l,bootindices)',wgood(bootindices)');
   % hold on
   % line([min(czstrat),max(czstrat)], ...
   %  [bbs(l,1)+bbs(l,2).*min(czstrat),bbs(l,1)+bbs(l,2).*max(czstrat)],'Color','r')
   % line([min(czstrat),max(czstrat)], ...
   %  [bbs_ols(l,1)+bbs_ols(l,2).*min(czstrat), ...
   %   bbs_ols(l,1)+bbs_ols(l,2).*max(czstrat)],'Color','k')
end

%Find 1 sigma and 2 sigma errors on robust linear trend
    bss_good = NaN.*zeros(length([-zszp:round(max(czstrat)-1)]),nboots);
for zstrat = ceil(min(czstrat)):floor(max(czstrat))
    sorted_predictors = sort(bbs_good(:,1)+bbs_good(:,2).*zstrat);
    bss_good(zstrat+zszp,:) = sorted_predictors; %bootstrap output sorted
    sorted_predictors_no_zs_errors = sort(bbs_good_no_zs_errors(:,1)+bbs_good_no_zs_errors(:,2).*zstrat);
    bss_good_no_zs_errors(zstrat+zszp,:) = sorted_predictors_no_zs_errors;
end


%Preprocessing part 2. Use K-S to find the breakpoint (is it at the strat boundary?).
%Then find the mean and s.d. of channels above and below that breakpoint.
for zstrat = ceil(min(czstrat)):floor(max(czstrat));
    [null,pksnominal(zstrat+zszp),ks2stat(zstrat+zszp)] = ...
     kstest2(s(11).cw(czstrat>zstrat),s(11).cw(czstrat<=zstrat));
end
    pksnominal(pksnominal==0) = NaN;
    [Cnominal,Inominal] = min(pksnominal); %Position of paleohydrologic break 
    %Find unit averages for hydrologically-defined units
    wunit(1,1) =    nanmean(s(11).cw(czstrat <(Inominal-zszp))); %Mean channel width above and below paleohydrologic break
    wunit(2,1) =    nanmean(s(11).cw(czstrat>=(Inominal-zszp))); %Mean channel width above and below paleohydrologic break
    wunit(1,2) =  nanmedian(s(11).cw(czstrat <(Inominal-zszp))); %Median channel width above and below paleohydrologic break
    wunit(2,2) =  nanmedian(s(11).cw(czstrat>=(Inominal-zszp))); %Median channel width above and below paleohydrologic break
    wunit(1,3) =     nanstd(s(11).cw(czstrat <(Inominal-zszp))); %s.d. of channel width above and below paleohydrologic break
    wunit(2,3) =     nanstd(s(11).cw(czstrat>=(Inominal-zszp))); %s.d. of channel width above and below paleohydrologic break


    wunit0(1,1) =    nanmean(s(11).cw(czstrat <(0))); %Mean channel width above and below paleohydrologic break
    wunit0(2,1) =    nanmean(s(11).cw(czstrat>=(0))); %Mean channel width above and below paleohydrologic break
    wunit0(1,2) =  nanmedian(s(11).cw(czstrat <(0))); %Median channel width above and below paleohydrologic break
    wunit0(2,2) =  nanmedian(s(11).cw(czstrat>=(0))); %Median channel width above and below paleohydrologic break
    wunit0(1,3) =     nanstd(s(11).cw(czstrat <(0))); %s.d. of channel width above and below paleohydrologic break
    wunit0(2,3) =     nanstd(s(11).cw(czstrat>=(0))); %s.d. of channel width above and below paleohydrologic break
    
    
randompool = 0.5+rand(nboots,length(czstrat))*length(czstrat);
    
for l = 1:nboots
    l
    bootindices = round(randompool(l,:)); %resample with replacement
    bszs_no_zs_offset =  czstrat(bootindices)'; %bootstrap zstrat values
    random_zs_offset_this_bootstrap(l,:) = sigmax.*randn(length(czstrat),1);
    random_cw_offset_this_bootstrap(l,:) = s(11).cwdev'.*randn(length(czstrat),1);

    bszs =  czstrat(bootindices)'+random_zs_offset_this_bootstrap(l,bootindices)'; %bootstrap zstrat values
 
    bootstrapped_proxy_data_values = s(11).cw(bootindices)'+random_cw_offset_this_bootstrap(l,bootindices)';
    bootstrapped_proxy_data_values(bootstrapped_proxy_data_values<0.1) = 0.1; %No negative channel widths!
    bsw = bootstrapped_proxy_data_values;
    
    
    for zstrat = ceil(min(czstrat)):floor(max(czstrat));
        if sum(bszs>zstrat) > 0 & sum(bszs<=zstrat) > 0
        [null,pks(l,zstrat+zszp),ks2stat(zstrat+zszp)] = ...
         kstest2(bsw(bszs>zstrat),bsw(bszs<=zstrat));
        end
        %Control with no offset:
        if sum(bszs_no_zs_offset>zstrat) > 0 & sum(bszs_no_zs_offset<=zstrat) > 0
        [null,pks_no_zs_offset(l,zstrat+zszp),ks2stat_no_zs_offset(zstrat+zszp)] = ...
         kstest2(bsw(bszs_no_zs_offset>zstrat),bsw(bszs_no_zs_offset<=zstrat));
        end
    end
    pks(pks==0) = NaN;
    [C(l),I(l)] = min(pks(l,:)); %Position of paleohydraulic break *with bootstrapping*
    pks_no_zs_offset(pks_no_zs_offset==0) = NaN;
    [C_no_zs_offset(l),I_no_zs_offset(l)] = min(pks_no_zs_offset(l,:));
end

%Preprocessing part 2b. Repeat part 2, but for only the data with pqual = 2('good') or pqual = 1('gold standard')
for zstrat = ceil(min(czstrat)):floor(max(czstrat));
     if sum(zgood>zstrat) > 0 & sum(zgood<=zstrat) > 0
    [null,pksnominal_good(zstrat+zszp),ks2stat(zstrat+zszp)] = ...
     kstest2(wgood(zgood>zstrat),wgood(zgood<=zstrat));
     end
end
    pksnominal_good(pksnominal_good==0) = NaN;
    [Cnominal_good,Inominal_good] = min(pksnominal_good); %Position of paleohydrologic break 
    %Find unit averages for hydrologically-defined units 
    wunit_good(1,1) =    nanmean(wgood(zgood <(Inominal_good-zszp))); %Mean channel width above and below paleohydrologic break
    wunit_good(2,1) =    nanmean(wgood(zgood>=(Inominal_good-zszp))); %Mean channel width above and below paleohydrologic break
    wunit_good(1,2) =  nanmedian(wgood(zgood <(Inominal_good-zszp))); %Median channel width above and below paleohydrologic break
    wunit_good(2,2) =  nanmedian(wgood(zgood>=(Inominal_good-zszp))); %Median channel width above and below paleohydrologic break
    wunit_good(1,3) =     nanstd(wgood(zgood <(Inominal_good-zszp))); %s.d. of channel width above and below paleohydrologic break
    wunit_good(2,3) =     nanstd(wgood(zgood>=(Inominal_good-zszp))); %s.d. of channel width above and below paleohydrologic break


%********* END OF POSTPROCESSING *******


%********* FIGURES *******
assert(min(s(11).cw(s(11).cw>0) - s(11).cwdev(s(11).cw>0))>3.5) %Check xlim doesn't clip out narrow channels

%show bootstrap positions:
    figure
    plot(cumsum(hist(2.*sort(I-zszp),[-200:1:200]))./max(cumsum(hist(2.*sort(I-zszp),[0:1:800]))), [-200:1:200], ...
    'LineWidth',2); 
    set(gca,'FontSize',18); 
    title('Channel width breakpoints')
    xlabel({'probability of being';'above breakpoint'});%,'Interpreter','latex');
    %ylabel('elevation relative to F1/F2 contact (m)');%,'Interpreter','latex'); 
    ylim([min(czstrat)-5 max(czstrat)+5]);%    scatter(zs,2.*[rhm(:).lambda],30,log10([rhm(:).quality])./log(2),'Filled')
    box on;grid off; set(gca,'FontSize',18);
 %   xlabel('Corrected elevation');ylabel(' meander wavelength (m)')
 %    title({cell2mat(strcat(zstrat_method,' correction,')); ...
 %         'meander wavelength versus corrected elevation: color is log2(quality)'})
        saveas(gca,cell2mat(strcat('LPSC2014ChannelWidth_breakpoints_',zstrat_method,'.eps')),'epsc');


figure('units','normalized','outerposition',[0 0 1 1])
    herrorbar(s(11).cw,cz,s(11).cwdev,'xk')
    hold on;grid off; set(gca,'FontSize',18)
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    scatter(s(11).cw,cz,30,cx,'LineWidth',2); xlim([3.5 80]);
    title({cell2mat(strcat(zstrat_method,' raw elevation'));'Color corresponds to E position'})
    set(gca,'xscale','log'); set(gca,'xtick',[5,10,15,20,30,40,60,80]);
    xlabel('Channel width (m)'); ylabel('Raw elevation (m)')
    set(gca,'yminortick','on')
                    saveas(gca,cell2mat(strcat('LPSC2014Width_1_',zstrat_method,'.eps')),'epsc');
    
figure('units','normalized','outerposition',[0 0 1 1])
    herrorbar(s(11).cw, cz, s(11).cwdev,'xk')
    hold on;grid off;    set(gca,'FontSize',18)
    title({cell2mat(strcat(zstrat_method,' raw elevation'));'color is log2(pquality)'})
    ylabel('Raw elevation (m)'); xlabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    set(gca,'xscale','log'); set(gca,'xtick',[5,10,15,20,30,40,60,80]);
    set(gca,'yminortick','on')
    scatter(s(11).cw,cz,30,log(cqual)./log(2),'LineWidth',2); xlim([3.5 80]); 
                    saveas(gca,cell2mat(strcat('LPSC2014Width_2_',zstrat_method,'.eps')),'epsc');


figure('units','normalized','outerposition',[0 0 1 1])
    herrorbar(s(11).cw, czstrat, s(11).cwdev,'xk')
    hold on;%colorbar;
    grid off
    set(gca,'FontSize',18)
    title({cell2mat(strcat(zstrat_method,' corrected elevation'));'color is x position'})
    ylabel('Elevation relative to F1/F2 contact (m)'); xlabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    set(gca,'xscale','log'); set(gca,'xtick',[5,10,15,20,30,40,60,80]); set(gca,'yminortick','on')
    scatter(s(11).cw,czstrat ,30,cx,'LineWidth',2); xlim([3.5 80]);ylim([min(czstrat)-5 max(czstrat)+5]);  
                    saveas(gca,cell2mat(strcat('LPSC2014Width_3_',zstrat_method,'.eps')),'epsc');

    
%4: Channel width, show pqual.    
figure('units','normalized','outerposition',[0 0 1 1])
    %bootstrap linear regressions:    
    herrorbar(s(11).cw,czstrat,s(11).cwdev,'xk')
    hold on
    set(gca,'FontSize',18); ylim([min(czstrat)-5 max(czstrat)+5]); xlim([3.5 80]);
    title({cell2mat(strcat(zstrat_method,' correction, one point per pick on F1-F2 contact,'));'color is log2(quality)'})
    ylabel('Stratigraphic elevation relative to F1-F2 contact (m)'); xlabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    set(gca,'xscale','log'); set(gca,'xtick',[5,10,15,20,30,40,60,80]); set(gca,'yminortick','on')
    scatter(s(11).cw,czstrat,30,log(cqual)./log(2),'LineWidth',2)
                        saveas(gca,cell2mat(strcat('LPSC2014Width_4_no_trendlines_',zstrat_method,'.eps')),'epsc');
    plot(bss(:,round(size(bss,2)*0.50)),[-zszp:round(max(czstrat)-1)],'r','LineWidth',1)
    hold on;grid off;colorbar
    plot(bss(:,round(size(bss,2)*0.05)),[-zszp:round(max(czstrat)-1)],'k','LineWidth',1)
    plot(bss(:,round(size(bss,2)*0.95)),[-zszp:round(max(czstrat)-1)],'k','LineWidth',1)
%    line([0.01 100],[0 0])
%    line([nanmean(s(11).cw(czstrat< 0)) nanmean(s(11).cw(czstrat< 0))],[min(czstrat) 0])
%    line([nanmean(s(11).cw(czstrat>=0)) nanmean(s(11).cw(czstrat>=0))],[max(czstrat) 0])
     line([0.01 100],[Inominal-zszp Inominal-zszp])
     line([wunit(1,2) wunit(1,2)],[min(czstrat) Inominal-zszp])
     line([wunit(2,2) wunit(2,2)],[max(czstrat) Inominal-zszp])
    scatter(s(11).cw,czstrat,30,log(cqual)./log(2),'LineWidth',2)
                    saveas(gca,cell2mat(strcat('LPSC2014Width_4_',zstrat_method,'.eps')),'epsc');

%5: Channel width, show pstyle.    
figure('units','normalized','outerposition',[0 0 1 1])    
    herrorbar(s(11).cw,czstrat,s(11).cwdev,'xk')
    hold on
    set(gca,'FontSize',18); ylim([min(czstrat)-5 max(czstrat)+5]); xlim([3.5 80]);
    title({cell2mat(strcat(zstrat_method,' correction, one point per pick on F1-F2 contact,'));'color is log2(pstyle)'})
    ylabel('Stratigraphic elevation relative to F1-F2 contact (m)'); xlabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
    set(gca,'xscale','log'); set(gca,'xtick',[5,10,15,20,30,40,60,80]); set(gca,'yminortick','on')
    scatter(s(11).cw,czstrat,30,log(cpstyle)./log(2),'LineWidth',2)
                        saveas(gca,cell2mat(strcat('LPSC2014Width_5_no_trendlines_',zstrat_method,'.eps')),'epsc');
        %bootstrap linear regressions:
    plot(bss(:,round(size(bss,2)*0.50)),[-zszp:round(max(czstrat)-1)],'r','LineWidth',1)
    hold on;grid off;colorbar
    plot(bss(:,round(size(bss,2)*0.05)),[-zszp:round(max(czstrat)-1)],'k','LineWidth',1)
    plot(bss(:,round(size(bss,2)*0.95)),[-zszp:round(max(czstrat)-1)],'k','LineWidth',1)
    %line([0.01 100],[0 0])
    %line([nanmean(s(11).cw(czstrat< 0)) nanmean(s(11).cw(czstrat< 0))],[min(czstrat) 0])
    %line([nanmean(s(11).cw(czstrat>=0)) nanmean(s(11).cw(czstrat>=0))],[max(czstrat) 0])
     line([0.01 100],[Inominal-zszp Inominal-zszp])
     line([wunit(1,2) wunit(1,2)],[min(czstrat) Inominal-zszp])
     line([wunit(2,2) wunit(2,2)],[max(czstrat) Inominal-zszp])
                    saveas(gca,cell2mat(strcat('LPSC2014Width_5_',zstrat_method,'.eps')),'epsc');


%6: Channel width, pqual=1 or 2 only, compare all data to pqual={1,2} only, show pqual.
figure('units','normalized','outerposition',[0 0 1 1])
    %show fits for 'all data' grayed out:
    plot(bss(:,round(size(bss,2)*0.50)),[-zszp:round(max(czstrat)-1)],'Color',[0.6 0.6 0.6],'LineWidth',1)
    hold on;grid off;colorbar; set(gca,'yminortick','on')
    plot(bss(:,round(size(bss,2)*0.05)),[-zszp:round(max(czstrat)-1)],'Color',[0.6 0.6 0.6],'LineWidth',1)
    plot(bss(:,round(size(bss,2)*0.95)),[-zszp:round(max(czstrat)-1)],'Color',[0.6 0.6 0.6],'LineWidth',1)
%    line([0.01 100],[0 0],'Color',[0.6 0.6 0.6])
%    line([nanmean(s(11).cw(czstrat< 0)) nanmean(s(11).cw(czstrat< 0))],[min(czstrat) 0],'Color',[0.6 0.6 0.6])
%    line([nanmean(s(11).cw(czstrat>=0)) nanmean(s(11).cw(czstrat>=0))],[max(czstrat) 0],'Color',[0.6 0.6 0.6])
     line([0.01 100],[Inominal-zszp Inominal-zszp],'Color',[0.6 0.6 0.6])
     line([wunit(1,2) wunit(1,2)],[min(czstrat) Inominal-zszp],'Color',[0.6 0.6 0.6])
     line([wunit(2,2) wunit(2,2)],[max(czstrat) Inominal-zszp],'Color',[0.6 0.6 0.6])
    
    %show fits for pqual={1,2} only:
    plot(bss_good(:,round(size(bss_good,2)*0.50)),[-zszp:round(max(czstrat)-1)],'r','LineWidth',1)
    plot(bss_good(:,round(size(bss_good,2)*0.05)),[-zszp:round(max(czstrat)-1)],'k','LineWidth',1)
    plot(bss_good(:,round(size(bss_good,2)*0.95)),[-zszp:round(max(czstrat)-1)],'k','LineWidth',1)
%    line([0.01 100],[0 0])
%    line([nanmean(s(11).cw(czstrat< 0 & cqual<=okq)) nanmean(s(11).cw(czstrat< 0 & cqual<=okq))],[min(czstrat) 0])
%    line([nanmean(s(11).cw(czstrat>=0 & cqual<=okq)) nanmean(s(11).cw(czstrat>=0 & cqual<=okq))],[max(czstrat) 0])
     line([0.01 100],[Inominal_good-zszp Inominal_good-zszp])
     line([wunit_good(1,2) wunit_good(1,2)],[min(czstrat) Inominal_good-zszp])
     line([wunit_good(2,2) wunit_good(2,2)],[max(czstrat) Inominal_good-zszp])

    
    herrorbar(s(11).cw(cqual<=okq),czstrat(cqual<=okq), ...
              s(11).cwdev(cqual<=okq),'xk')
    hold on;grid off;colorbar; set(gca,'yminortick','on')
    set(gca,'FontSize',18); ylim([min(czstrat)-5 max(czstrat)+5]); xlim([3.5 80])
    title({cell2mat(strcat(zstrat_method,' correction, one point per pick on F1-F2 contact,')); ...
                    strcat('color is log2(quality). Pqual<=',num2str(okq),' only')})
ylabel('Stratigraphic elevation relative to F1-F2 contact (m)'); xlabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
        set(gca,'xscale','log'); set(gca,'xtick',[5,10,15,20,30,40,60,80]);
    scatter(s(11).cw(cqual<=okq),czstrat(cqual<=okq), ...
            50,log(cqual(cqual<=okq))./log(2),'LineWidth',2)
                    saveas(gca,cell2mat(strcat('LPSC2014Width_6_',zstrat_method,'.eps')),'epsc');

                    
%7: Channel width, pqual=1 or 2 only, compare all data to pqual={1,2} only, show pstyle.
figure('units','normalized','outerposition',[0 0 1 1])
    %show fits for 'all data' grayed out:
    plot(bss(:,round(size(bss,2)*0.50)),[-zszp:round(max(czstrat)-1)],'Color',[0.6 0.6 0.6],'LineWidth',1)
    hold on;grid off;colorbar; set(gca,'yminortick','on')
    plot(bss(:,round(size(bss,2)*0.05)),[-zszp:round(max(czstrat)-1)],'Color',[0.6 0.6 0.6],'LineWidth',1)
    plot(bss(:,round(size(bss,2)*0.95)),[-zszp:round(max(czstrat)-1)],'Color',[0.6 0.6 0.6],'LineWidth',1)
%    line([0.01 100],[0 0],'Color',[0.6 0.6 0.6])
%    line([nanmean(s(11).cw(czstrat< 0)) nanmean(s(11).cw(czstrat< 0))],[min(czstrat) 0],'Color',[0.6 0.6 0.6])
%    line([nanmean(s(11).cw(czstrat>=0)) nanmean(s(11).cw(czstrat>=0))],[max(czstrat) 0],'Color',[0.6 0.6 0.6])
    line([0.01 100],[Inominal-zszp Inominal-zszp],'Color',[0.6 0.6 0.6])
    line([wunit(1,2) wunit(1,2)],[min(czstrat) Inominal-zszp],'Color',[0.6 0.6 0.6])
    line([wunit(2,2) wunit(2,2)],[max(czstrat) Inominal-zszp],'Color',[0.6 0.6 0.6])
    %show fits for pqual={1,2} only:
    plot(bss_good(:,round(size(bss_good,2)*0.50)),[-zszp:round(max(czstrat)-1)],'r','LineWidth',1)
    plot(bss_good(:,round(size(bss_good,2)*0.05)),[-zszp:round(max(czstrat)-1)],'k','LineWidth',1)
    plot(bss_good(:,round(size(bss_good,2)*0.95)),[-zszp:round(max(czstrat)-1)],'k','LineWidth',1)
%    line([0.01 100],[0 0])
%    line([nanmean(s(11).cw(czstrat< 0 & cqual<=okq)) nanmean(s(11).cw(czstrat< 0 & cqual<=okq))],[min(czstrat) 0])
%    line([nanmean(s(11).cw(czstrat>=0 & cqual<=okq)) nanmean(s(11).cw(czstrat>=0 & cqual<=okq))],[max(czstrat) 0])
     line([0.01 100],[Inominal_good-zszp Inominal_good-zszp])
     line([wunit_good(1,2) wunit_good(1,2)],[min(czstrat) Inominal_good-zszp])
     line([wunit_good(2,2) wunit_good(2,2)],[max(czstrat) Inominal_good-zszp])

    
    herrorbar(s(11).cw(cqual<=okq),czstrat(cqual<=okq), ...
              s(11).cwdev(cqual<=okq),'xk')
    hold on;grid off;colorbar
    set(gca,'FontSize',18); ylim([min(czstrat)-5 max(czstrat)+5]); xlim([3.5 80])
    title({cell2mat(strcat(zstrat_method,' correction, one point per pick on F1-F2 contact,')); ...
           strcat('color is log2(pstyle). Pqual<=',num2str(okq),' only')})
    ylabel('Stratigraphic elevation relative to F1-F2 contact (m)'); xlabel('Channel width (m)')
    %herrorbar(cz,cw,czdev,'LineSpec',' ') %scatter in NOMINAL ELEVATIONS
    %(not krided stratigraphic position!) is usually negligible
            set(gca,'xscale','log'); set(gca,'xtick',[5,10,15,20,30,40,60,80]);
    scatter(s(11).cw(cqual<=okq),czstrat(cqual<=okq), ...
            50,log(cpstyle(cqual<=okq))./log(2),'LineWidth',2)
                    saveas(gca,cell2mat(strcat('LPSC2014Width_7_',zstrat_method,'.eps')),'epsc');



%error('Not doing probability summation while debugging')

%Quick summation of probability distribution (assuming uncorrelated, Guassian error
%bars) From Wikipedia, "Multivariate normal distribution"
p = zeros(601,401);
warning('Assuming fixed sigmax for debug')
for nn = 1:length(cz) %add the probability kernel of each observation
oldp = p;
if rem(nn,10)==0; disp(nn);end
    if strcmp(zstrat_method,'planar') == 1
            sigmax = 29.084761000706; %RMS error from f1_f2_contact_second_pass_b20g02_only_plane_fit_one_point_per_contact_pick_24feb14_rms_output
    elseif strcmp(zstrat_method,'quadratic') == 1
            sigmax = 17.01479366; %RMS error
    elseif strcmp(zstrat_method,'IDW') == 1
            sigmax = 27.7; %cross-validation
    else
        error
    end
            
            sigmay = s(11).cwdev(nn);
            normfactor = 1/(2*pi*sigmax*sigmay);
            mux = czstrat(nn);
            muy = s(11).cw(nn);
            
             if sum(isnan([sigmax,sigmay,mux,muy]))==0
                 for x = -300:1:300 %height above F1/F2 contact (m)
                     for y = 1:1:300 %river channel width (m)
                             p(x+301,y+101) = p(x+301,y+101) + ...
                             log10(y).*normfactor*exp(-1/2*(((x-mux)^2/(sigmax^2))+((y-muy)^2/(sigmay^2))));
                     end
                 end
             end

            
            
%             if sum(isnan([sigmax,sigmay,mux,muy]))==0
%                 for x = -300:1:300 %height above F1/F2 contact (m)
%                     for y = -100:1:300 %river channel width (m)
%                             p(x+301,y+101) = p(x+301,y+101) + ...
%                             normfactor*exp(-1/2*(((x-mux)^2/(sigmax^2))+((y-muy)^2/(sigmay^2))));
%                     end
%                 end
%             end
                %For debug only
                %check that the additional probability sums to 1 (note that
                %the discretisation and the cutoffs will introduce errors!)
                %if abs((sum(p(:)) - sum(oldp(:)))-1) > 0.01
                %    error('MVN Dist Checksum Error!')
                %end
end

[X,Y] = meshgrid(-300:1:300,-100:1:300);

pcs = cumsum(p,1);
ps = sum(p');
for l = 1:601 %cumsum normalized
pcsn(l,:) = p(l,:)./ps(l);
end 
 saturate = 0.05;
 pcsn(pcsn>saturate)= saturate;

 
%Normalized probability density
figure
contourf(Y,X,pcsn',[0:.0005:saturate]);grid off %show the sum of
shading flat;cmap=flipud(colormap(gray));
colormap(cmap);
ylim([min(czstrat)-5 max(czstrat)+5]); xlim([3.5 80]);
hold on
    set(gca,'FontSize',18)
    title(strcat('Probability density normalizing to sum(strat. elev.): ',zstrat_method))
        ylabel('Stratigraphic elevation relative to F1-F2 contact (m)'); xlabel('Channel width (m)')
    hold on;grid off;
    errorbar(60,max(czstrat)-sigmax*1.2,sigmax,'r')
    herrorbar(s(11).cw,czstrat,s(11).cwdev,'xr')
    scatter(s(11).cw,czstrat,30,'r*','LineWidth',1)
    scatter(s(11).cw,czstrat,30,'r*','LineWidth',1)
    set(gca,'xscale','log'); set(gca,'xtick',[5,10,15,20,30,40,60,80]);
    saveas(gca,cell2mat(strcat('LPSC2014Width_8_notrendlines',zstrat_method,'.eps')),'epsc');
    plot(bss(:,round(size(bss,2)*0.50)),[-zszp:round(max(czstrat)-1)],'c','LineWidth',1)
    plot(bss(:,round(size(bss,2)*0.05)),[-zszp:round(max(czstrat)-1)],'m','LineWidth',1)
    plot(bss(:,round(size(bss,2)*0.95)),[-zszp:round(max(czstrat)-1)],'m','LineWidth',1)
%    line([0.01 100],[0 0],'Color','r','LineWidth',1)
%    line([nanmean(s(11).cw(czstrat< 0)) nanmean(s(11).cw(czstrat< 0))],[min(czstrat) 0],'Color','r','LineWidth',1)
%    line([nanmean(s(11).cw(czstrat>=0)) nanmean(s(11).cw(czstrat>=0))],[max(czstrat) 0],'Color','r','LineWidth',1)
     line([0.01 100],[Inominal-zszp Inominal-zszp])
     line([wunit(1,2) wunit(1,2)],[min(czstrat) Inominal-zszp])
     line([wunit(2,2) wunit(2,2)],[max(czstrat) Inominal-zszp])

    saveas(gca,cell2mat(strcat('LPSC2014Width_8_',zstrat_method,'.eps')),'epsc');

ps = sum(p);
for l = 1:401 %cumsum normalized
pcsnw(:,l) = p(:,l)./ps(l);
end 
 saturate = 0.015;
 pcsnw(pcsnw>saturate)= saturate;     

 %Normalized probability density
figure
    scatter(s(11).cw,czstrat,30,'r*','LineWidth',1)
contourf(Y,X,pcsnw',[0:0.0005:saturate]);grid off %show the sum of
shading flat;cmap=flipud(colormap(gray));
colormap(cmap);
ylim([min(czstrat)-5 max(czstrat)+5]); xlim([3.5 80]);
hold on
    herrorbar(s(11).cw,czstrat,s(11).cwdev,'xr')
    errorbar(60,max(czstrat)-sigmax*1.2,sigmax,'r')
    scatter(s(11).cw,czstrat,30,'r*','LineWidth',1)
    set(gca,'xscale','log'); set(gca,'xtick',[5,10,15,20,30,40,60,80]);
    set(gca,'FontSize',18)
    title(strcat('Probability density normalizing to sum(width): ',zstrat_method))
        ylabel('Stratigraphic elevation relative to F1-F2 contact (m)'); xlabel('Channel width (m)')
    hold on;grid off;; set(gca,'yminortick','on')
    scatter(s(11).cw,czstrat,30,'r*','LineWidth',1)
    saveas(gca,cell2mat(strcat('LPSC2014Width_9_notrendlines',zstrat_method,'.eps')),'epsc');
    plot(bss(:,round(size(bss,2)*0.50)),[-zszp:round(max(czstrat)-1)],'c','LineWidth',1)
    plot(bss(:,round(size(bss,2)*0.05)),[-zszp:round(max(czstrat)-1)],'m','LineWidth',1)
    plot(bss(:,round(size(bss,2)*0.95)),[-zszp:round(max(czstrat)-1)],'m','LineWidth',1)
%    line([0.01 100],[0 0],'Color','r','LineWidth',1)
%    line([nanmean(s(11).cw(czstrat< 0)) nanmean(s(11).cw(czstrat< 0))],[min(czstrat) 0],'Color','r','LineWidth',1)
%    line([nanmean(s(11).cw(czstrat>=0)) nanmean(s(11).cw(czstrat>=0))],[max(czstrat) 0],'Color','r','LineWidth',1)
     line([0.01 100],[Inominal-zszp Inominal-zszp])
     line([wunit(1,2) wunit(1,2)],[min(czstrat) Inominal-zszp])
     line([wunit(2,2) wunit(2,2)],[max(czstrat) Inominal-zszp])

    saveas(gca,cell2mat(strcat('LPSC2014Width_9_',zstrat_method,'.eps')),'epsc');

 

%Raw probability density
figure
contourf(Y,X,p',[0:0.001:max(p(:))]);grid off %show the sum of
shading flat;cmap=flipud(colormap(gray));
colormap(cmap);
ylim([min(czstrat)-5 max(czstrat)+5]); xlim([3.5 80]);
hold on
    set(gca,'FontSize',18)
    title(strcat('Raw probability density: ',zstrat_method))
        ylabel('Stratigraphic elevation relative to F1-F2 contact (m)'); xlabel('Channel width (m)')
    hold on;grid off; set(gca,'yminortick','on')
    herrorbar(s(11).cw,czstrat,s(11).cwdev,'xr')
    scatter(s(11).cw,czstrat,30,'r*','LineWidth',1)
    errorbar(60,max(czstrat)-sigmax*1.2,sigmax,'r')

    saveas(gca,cell2mat(strcat('LPSC2014Width_10_notrendlines',zstrat_method,'.eps')),'epsc');
    plot(bss(:,round(size(bss,2)*0.50)),[-zszp:round(max(czstrat)-1)],'c','LineWidth',1)
    plot(bss(:,round(size(bss,2)*0.05)),[-zszp:round(max(czstrat)-1)],'m','LineWidth',1)
    plot(bss(:,round(size(bss,2)*0.95)),[-zszp:round(max(czstrat)-1)],'m','LineWidth',1)
%    line([0.01 100],[0 0],'Color','r','LineWidth',1)
%    line([nanmean(s(11).cw(czstrat< 0)) nanmean(s(11).cw(czstrat< 0))],[min(czstrat) 0],'Color','r','LineWidth',1)
%    line([nanmean(s(11).cw(czstrat>=0)) nanmean(s(11).cw(czstrat>=0))],[max(czstrat) 0],'Color','r','LineWidth',1)
     line([0.01 100],[Inominal-zszp Inominal-zszp])
     line([wunit(1,2) wunit(1,2)],[min(czstrat) Inominal-zszp])
     line([wunit(2,2) wunit(2,2)],[max(czstrat) Inominal-zszp])
    set(gca,'xscale','log'); set(gca,'xtick',[5,10,15,20,30,40,60,80]);
    saveas(gca,cell2mat(strcat('LPSC2014Width_10_',zstrat_method,'.eps')),'epsc');

    
figure;plot([0:1:60],hist(squeeze(wunitbs(:,1,1)),[0:1:60]),'LineWidth',2);hold on;
plot([0:1:60],hist(squeeze(wunitbs(:,2,1)),[0:1:60]),'c','LineWidth',2);hold on; 
set(gca,'xminortick','on'); set(gca,'yminortick','on')
line([nanmean(s(11).cw(czstrat>=0)) nanmean(s(11).cw(czstrat>=0))],[0 300],'Color','c','LineWidth',2)
line([nanmean(s(11).cw(czstrat< 0)) nanmean(s(11).cw(czstrat< 0))],[0 300],'LineWidth',2)

set(gca,'FontSize',18)
xlabel('Mean channel width (m)')
ylabel('# of bootstraps')
     saveas(gca,cell2mat(strcat('LPSC2014Width_11_',zstrat_method,'.eps')),'epsc');
   
end
