Skip to content

Commit

Permalink
Reduces "bayesShadedPlot" to use only the "mean" fit option (#277)
Browse files Browse the repository at this point in the history 
* Removes "fit" option from "bayesShadedPlot"

* Tidies up "bayesShadedPlot"

* Addresses review comments
  • Loading branch information
@DrPaulSharp
DrPaulSharp authored Nov 11, 2024
1 parent 2d2a212 commit c19b35d
Show file tree
Hide file tree
Showing 5 changed files with 76 additions and 146 deletions.
Binary file modified examples/miscellaneous/convertRasCAL1Project/convertRascal.mlx
View file
Binary file not shown.
14 changes: 8 additions & 6 deletions examples/normalReflectivity/customXY/customXYDSPCScript.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,7 +188,7 @@
% ..and plot this out....

figure(30); clf;
bayesShadedPlot(problem, results,'fit','average','KeepAxes',true,'interval',65,'q4',false)
bayesShadedPlot(problem, results,'KeepAxes',true,'interval',65,'q4',false)

h3 = figure(40); clf
plotHists(results,h3,'smooth',true)
Expand Down Expand Up @@ -271,12 +271,14 @@

% In RAT, there is a useful function called 'shade' that we can use
% here.....
fillType = [1 2;2 1];
fillAlpha = 0.3;
cols = get(gca,'ColorOrder');
shade(z,ciSi(1,:),z,ciSi(2,:),'FillColor',cols(1,:),'FillType',[1 2;2 1],'FillAlpha',0.3);
shade(z,ciOxide(1,:),z,ciOxide(2,:),'FillColor',cols(2,:),'FillType',[1 2;2 1],'FillAlpha',0.3);
shade(z,ciHeadL(1,:),z,ciHeadL(2,:),'FillColor',cols(3,:),'FillType',[1 2;2 1],'FillAlpha',0.3);
shade(z,ciTails(1,:),z,ciTails(2,:),'FillColor',cols(4,:),'FillType',[1 2;2 1],'FillAlpha',0.3);
shade(z,ciHeadR(1,:),z,ciHeadR(2,:),'FillColor',cols(5,:),'FillType',[1 2;2 1],'FillAlpha',0.3);
shade(z,ciSi(1,:),z,ciSi(2,:),'FillColor',cols(1,:),'FillType',fillType,'FillAlpha',fillAlpha);
shade(z,ciOxide(1,:),z,ciOxide(2,:),'FillColor',cols(2,:),'FillType',fillType,'FillAlpha',fillAlpha);
shade(z,ciHeadL(1,:),z,ciHeadL(2,:),'FillColor',cols(3,:),'FillType',fillType,'FillAlpha',fillAlpha);
shade(z,ciTails(1,:),z,ciTails(2,:),'FillColor',cols(4,:),'FillType',fillType,'FillAlpha',fillAlpha);
shade(z,ciHeadR(1,:),z,ciHeadR(2,:),'FillColor',cols(5,:),'FillType',fillType,'FillAlpha',fillAlpha);
title('Volume Fractions');
%%
% .. and we are done.
Binary file modified examples/normalReflectivity/customXY/customXYDSPCSheet.mlx
View file
Binary file not shown.
202 changes: 65 additions & 137 deletions utilities/plotting/bayesShadedPlot.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ function bayesShadedPlot(problem,result,varargin)

% Plot the shaded reflectivities from Bayes output
% from RAT
sf = result.contrastParams.scalefactors;

if isa(problem,'domainsClass')
isDomains = true;
else
Expand All @@ -11,60 +11,55 @@ function bayesShadedPlot(problem,result,varargin)

% Parse the input options
if ~isempty(varargin)
defaultq4 = 0;
defaultFit = 'mean';

defaultq4 = false;
defaultKeep = false;
defaultInterval = 95;


allIntervals = [65 95];

p = inputParser;
addOptional(p, 'q4', defaultq4, @islogical);
addOptional(p, 'fit', defaultFit, @(x) any(strcmpi(x,{'mean','average','all'})));
addOptional(p, 'KeepAxes', defaultKeep, @islogical);
addOptional(p, 'interval', defaultInterval, @(x) ismember(x,allIntervals));

parse(p,varargin{:});
inputBlock = p.Results;

q4 = inputBlock.q4;
fit = inputBlock.fit;
keepAx = inputBlock.KeepAxes;
interval = inputBlock.interval;

else

q4 = false;
fit = 'average';
keepAx = false;
interval = 65;
end

%mean %Average % All
showWhichCurves = [false false false];
interval = 95;

switch fit
case 'mean'
showWhichCurves(1) = 1;
case 'average'
showWhichCurves(2) = 1;
case 'all'
showWhichCurves(1:2) = 1;
end

f = gcf;

if ~ keepAx
if ~keepAx
clf; hold on; box on
end

pLims = result.predictionIntervals;
refPlims = pLims.reflectivity;
sldPlims = pLims.sld;
switch interval
case 95
vals = [1 5];
case 65
vals = [2 4];
end

fillColor = [0.7 0.7 0.7];
fillType = [1 2;2 1];
fillAlpha = 0.3;

% Get the reflectivities and SLDs
bestReflectivity = result.reflectivity;
bestSld = result.sldProfiles;

% Get the reflectivities for mean...
bestRefMean = result.reflectivity;
bestSldMean = result.sldProfiles;
reflectivityLimits = result.predictionIntervals.reflectivity;
sldLimits = result.predictionIntervals.sld;

shiftedData = result.shiftedData;
numberOfContrasts = length(shiftedData);
Expand All @@ -76,154 +71,87 @@ function bayesShadedPlot(problem,result,varargin)

for i = 1:numberOfContrasts

%thisRef = reflect{i};
thisData = shiftedData{i};
thisRefMean = bestRefMean{i};
thisSf = sf(i);
reflectivity = bestReflectivity{i};

mult = 2^(4*i);
switch q4
case true
thisQ4 = thisData(:,1).^4;
case false
if i == 1
mult = 1;
else
mult = 2^(4*i);
end
otherwise
mult = 2^(4*i);
thisQ4 = thisData(:,1).^4;
end

% Get the limits and fits
theseLims = refPlims{i};

switch interval
case 95
vals = [1 5];
case 65
vals = [2 4];
end

thisMin = theseLims(vals(1),:)./mult;
thisMax = theseLims(vals(2),:)./mult;

thisRefAvg = theseLims(3,:)./mult;
limits = reflectivityLimits{i};

thisRefMean(:,2) = thisRefMean(:,2)./mult;
%thisRefMax(:,2) = thisRefMax(:,2)./mult;
min = limits(vals(1),:)./mult;
max = limits(vals(2),:)./mult;

reflectivity(:,2) = reflectivity(:,2)./mult;

thisDataX = thisData(:,1);
thisDataY = thisData(:,2)./mult;
thisDataErr = thisData(:,3)./mult;
dataX = thisData(:,1);
dataY = thisData(:,2)./mult;
dataErr = thisData(:,3)./mult;

thisSimX = result.reflectivity{i}(:,1);
thisSimQ4 = thisSimX.^4;
refXValues = result.reflectivity{i}(:,1);
thisSimQ4 = refXValues.^4;

switch q4
case true
thisMin = thisMin .* thisSimQ4;
thisMax = thisMax .* thisSimQ4;
thisRefMean(:,2) = thisRefMean(:,2) .* thisQ4;
thisRefAvg = thisRefAvg .* thisSimQ4;
thisDataY = thisDataY(:) .* thisQ4;
thisDataErr = thisDataErr(:) .* thisQ4;
min = min .* thisSimQ4;
max = max .* thisSimQ4;
reflectivity(:,2) = reflectivity(:,2) .* thisQ4;
dataY = dataY(:) .* thisQ4;
dataErr = dataErr(:) .* thisQ4;
end

errorbar(thisDataX,thisDataY,thisDataErr,'.');
errorbar(dataX,dataY,dataErr,'.');
shade(refXValues,min,refXValues,max,'FillColor',fillColor,'FillType',fillType,'FillAlpha',fillAlpha);
plot(reflectivity(:,1),reflectivity(:,2),'b-');

%plot(thisData(:,1),thisMin,'-','color',[0.7 0.7 0.7]);
%plot(thisData(:,1),thisMax,'-','color',[0.7 0.7 0.7]);
%plot(thisData(:,1),thisRef,'LineWidth',2.0);
%fillyy(thisData(:,1),thisMin,thisMax,[0.8 0.8 0.8]);
shade(thisSimX,thisMin,thisSimX,thisMax,'FillColor',[0.7 0.7 0.7],'FillType',[1 2;2 1],'FillAlpha',0.3);

% Plot the requested fit lines;
if showWhichCurves(1)
% Plot the mean
plot(thisRefMean(:,1),thisRefMean(:,2),'b-');
end

if showWhichCurves(2)
% Plot the average
plot(thisDataX,thisRefAvg,'r-');
end

end

% Now plot the SLD's
% Now plot the SLDs
subplot(1,2,2); hold on; box on

if ~isDomains
for i = 1:numberOfContrasts

thisSldMean = bestSldMean{i};
%thisSldMax = bestSld_max{i};

theseLims = sldPlims{i};

thisSldX = result.sldProfiles{i}(:,1);

switch interval
case 95
vals = [1 5];
case 65
vals = [2 4];
end
for i = 1:numberOfContrasts

thisMin = theseLims(vals(1),:);
thisMax = theseLims(vals(2),:);
sld = bestSld{i};
limits = sldLimits{i};
sldXValues = result.sldProfiles{i}(:,1);

thisSldAvg = theseLims(3,:);
min = limits(vals(1),:);
max = limits(vals(2),:);

if showWhichCurves(1)
% Plot the mean
plot(thisSldMean(:,1),thisSldMean(:,2),'b-');
end
plot(sld(:,1),sld(:,2),'b-');
shade(sldXValues,min,sldXValues,max,'FillColor',fillColor,'FillType',fillType,'FillAlpha',fillAlpha);

if showWhichCurves(2)
% Plot the max
plot(thisSldX,thisSldAvg,'r-');
end

shade(thisSldX,thisMin,thisSldX,thisMax,'FillColor',[0.7 0.7 0.7],'FillType',[1 2;2 1],'FillAlpha',0.3);
end
else
for i = 1:numberOfContrasts

thisSldMean = bestSldMean(i,:);
%thisSldMax = bestSld_max{i};

theseLims = sldPlims(i,:);

thisSldX = result.sldProfiles(i,:);

switch interval
case 95
vals = [1 5];
case 65
vals = [2 4];
end
else

for m = 1:2
thisMin = theseLims{m}(vals(1),:);
thisMax = theseLims{m}(vals(2),:);
for i = 1:numberOfContrasts

thisDomainSldX = thisSldX{m}(:,1);
thisSldAvg = theseLims{m}(3,:);
sld = bestSld(i,:);
limits = sldLimits(i,:);
sldXValues = result.sldProfiles(i,:);

if showWhichCurves(1)
% Plot the mean
plot(thisSldMean{m}(:,1),thisSldMean{m}(:,2),'b-');
end
for j = 1:2
min = limits{j}(vals(1),:);
max = limits{j}(vals(2),:);

if showWhichCurves(2)
% Plot the max
plot(thisDomainSldX,thisSldAvg,'r-');
end
thisDomainSldXValues = sldXValues{j}(:,1);

shade(thisDomainSldX,thisMin,thisDomainSldX,thisMax,'FillColor',[0.7 0.7 0.7],'FillType',[1 2;2 1],'FillAlpha',0.3);
plot(sld{j}(:,1),sld{j}(:,2),'b-');
shade(thisDomainSldXValues,min,thisDomainSldXValues,max,'FillColor',fillColor,'FillType',fillType,'FillAlpha',fillAlpha);
end
end

end

end
end
6 changes: 3 additions & 3 deletions utilities/plotting/plotBayes.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,13 +3,13 @@ function plotBayes(problem, results)
figure(10); clf; plotRefSLD(problem,results)

figure(30); clf;
bayesShadedPlot(problem,results,'fit','mean','KeepAxes',true,'interval',95,'q4',false)
bayesShadedPlot(problem,results,'KeepAxes',true,'interval',95,'q4',false)

h3 = figure(40); clf
plotHists(results,h3,'smooth',true)

% h4 = figure(5); clf
% cornerPlot(results,h4,'smooth',false)
% h4 = figure(5); clf
% cornerPlot(results,h4,'smooth',false)

figure(60); clf
mcmcplot(results.chain,[],results.fitNames,'chainpanel');
Expand Down

0 comments on commit c19b35d

Please sign in to comment.