Index: /issm/trunk/src/m/classes/public/plot/plot_qmuhistnorm.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_qmuhistnorm.m (revision 345) +++ /issm/trunk/src/m/classes/public/plot/plot_qmuhistnorm.m (revision 345) @@ -0,0 +1,209 @@ +% +% plot a relative histogram along with the normal distribution. +% +% []=plot_hist_norm(rfunc1 ,rfunc2 ,hmin,hmax,hnint) +% []=plot_hist_norm(data,desc,mu,sigma,hmin,hmax,hnint) +% +function []=plot_qmuhistnorm(varargin) + +if ~nargin + help plot_hist_norm + return +end + +%% process input data and assemble into matrices as needed + +ivar=1; +if isstruct(varargin{ivar}) + rfunc1=varargin{ivar}; + ivar=ivar+1; + + desc=cell (1,length(rfunc1)); + ldata=zeros(1,length(rfunc1)); + for i=1:length(rfunc1) + ldata(i)=length(rfunc1(i).sample); + end + data=zeros(max(ldata),length(rfunc1)); + data(:,:)=NaN; + + for i=1:length(rfunc1) + desc(i)=cellstr(rfunc1(i).descriptor); + data(1:ldata(i),i)=rfunc1(i).sample; + end +else + data=varargin{ivar}; + ivar=ivar+1; + + ldata(1:size(data,2))=size(data,1); + + if ivar <= nargin && iscell(varargin{ivar}) + desc=varargin{ivar}; + ivar=ivar+1; + else + desc={}; + end +end + +if ivar <= nargin && isstruct(varargin{ivar}) + rfunc2=varargin{ivar}; + ivar=ivar+1; + + mu =zeros(1,length(rfunc2)); + sigma=zeros(1,length(rfunc2)); + + for i=1:length(rfunc2) + mu (i)=rfunc2(i).mean; + sigma(i)=rfunc2(i).stddev; + end +elseif ivar+1 <= nargin + if isnumeric(varargin{ivar}) + if ~isempty(varargin{ivar}) + mu =varargin{ivar}; + else + mu =mean(data); + display('Using calculated means.') + end + end + ivar=ivar+1; + if isnumeric(varargin{ivar}) + if ~isempty(varargin{ivar}) + sigma=varargin{ivar}; + else + sigma=std(data); + display('Using calculated standard deviations.') + end + end + ivar=ivar+1; +end + +%% generate the intervals + +if ivar <= nargin && ~isempty(varargin{ivar}) + hmin=varargin{ivar}; +else + hmin=min(min(data)); +end +ivar=ivar+1; +if ivar <= nargin && ~isempty(varargin{ivar}) + hmax=varargin{ivar}; +else + hmax=max(max(data)); +end +ivar=ivar+1; +if ivar <= nargin && ~isempty(varargin{ivar}) + hnint=varargin{ivar}; +else + hnint=50; +end +ivar=ivar+1; +edges=hmin:(hmax-hmin)/hnint:hmax; + +%% generate the histogram counts and make them relative + +% note that for the histc function: +% n(k) counts the value x(i) if edges(k) <= x(i) < edges(k+1). +% The last bin counts any values of x that match edges(end). +% Values outside the values in edges are not counted. +% Use -inf and inf in edges to include all non-NaN values. + +dhistc=histc(data,edges); +for i=1:size(data,2) + dbelow(i) =length(find(data(:,i)edges(end)))/ldata(i); +end + +if exist('mu','var') && exist('sigma','var') + ncol=size(data,2); + for i=1:ncol + dbelow(ncol+i)=normcdf(edges( 1),mu(i),sigma(i)); + dhistc(1:size(dhistc,1)-1,ncol+i)=... + normcdf(edges(2:end ),mu(i),sigma(i))-... + normcdf(edges(1:end-1),mu(i),sigma(i)); + dabove(ncol+i)=norminv(edges(end),mu(i),sigma(i)); + if exist('desc','var') + desc(ncol+i)={[desc{i} ' norm']}; + end + end +end + +% draw the bar plot + +hl1=bar(edges(1:end-1),dhistc(1:end-1,:)); +ax1=gca; + +% add the annotation + +title('Relative Frequency Histogram') +xlabel('Interval Edge Value') +ylabel('Relative Frequency') + +if exist('desc','var') + hleg1=legend(ax1,desc,'Location','NorthWest',... + 'Interpreter','none'); +else + hleg1=legend(ax1); +end + +%% generate the cumulative distribution functions + +% cdf=zeros(size(dhistc)); +% cdf(1,:)=dhistc(1,:); +% for i=2:size(dhistc,1) +% cdf(i,:)=cdf(i-1,:)+dhistc(i,:); +% end +cdf=cumsum(dhistc); +for i=1:size(dhistc,2) + cdf(:,i)=dbelow(i)+cdf(:,i); +end +if exist('desc','var') + ncol=length(desc); + for i=1:ncol + cdesc(i)={[desc{i} ' cdf']}; + end +end + + +% draw the line plot + +% (see "Using Multiple X- and Y-Axes" and "Overlaying Other +% Plots on Bar Graphs", or search on "YAxisLocation right") + +%hold all +%hold on +%plot(edges,cdf) +%plotyy([],[],edges,cdf) + +% ticks from the bar plot will show through, so make ticks +% from the line plot match + +nytick=length(get(ax1,'YTick')); +ytmax =ceil(max(max(cdf))/0.1-0.1)*0.1; +ytinc =ytmax/(nytick-1-1); + +ax2=axes('Position',get(ax1,'Position'),... + 'XLim',get(ax1,'XLim'),... + 'YLim',[0 ytinc*(nytick-1)],... + 'YTick',[0:ytinc:ytinc*(nytick-1)],... + 'XAxisLocation','bottom','YAxisLocation','right',... + 'Color','none','Layer','top'); +hl2=line(edges(1:end-1),cdf(1:end-1,:),'Parent',ax2); + +% add the annotation + +ylabel('Cumulative Percent') + +% legend doesn't combine with bar chart above +if exist('cdesc','var') + hleg2=legend(ax2,cdesc,'Location','NorthEast',... + 'Interpreter','none'); + set(hleg2,'Color','white') +else + hleg2=legend(ax2); + set(hleg2,'Color','white') +end + +set(gcf,'PaperPositionMode','auto') +%hold off + +end Index: /issm/trunk/src/m/classes/public/plot/plot_qmuimportancefactors.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_qmuimportancefactors.m (revision 345) +++ /issm/trunk/src/m/classes/public/plot/plot_qmuimportancefactors.m (revision 345) @@ -0,0 +1,87 @@ +function plot_qmuimportancefactors(md,options_structure,width,ii); +%PLOT_IMPORTANCEFACTORS - plot importance factors +% +% Usage: +% plot_importancefactors(md,options_structure,width,i); +% +% See also: PLOTMODEL + +%first recover design variable descriptor +if ~isnan(options_structure.designvariable), + descriptor=options_structure.designvariable; +else + error('plot_importancefactors error message: Need to supply design variable descriptor'); +end +descriptorlength=length(descriptor); + +%then recover responsfunction name +if ~isnan(options_structure.responsefunction), + responsefunctiondescriptor=options_structure.responsefunction; +else + error('plot_importancefactors error message: Need to supply response function descriptor'); +end + +%go through all response functions and find the one corresponding to the correct responsefunctiondescriptor +responsefunctions=md.dakotaresults.dresp_out; +found=0; +for i=1:length(responsefunctions), + if strcmpi(responsefunctions(i).descriptor,responsefunctiondescriptor), + found=i; + break; + end +end +if ~found, + error('plot_importancefactors error message: could not find correct response function'); +end +responsefunctions=responsefunctions(found); +nfun=size(responsefunctions.desvar,1); + +%Now recover response to the correct design variable +importancefactors=zeros(md.npart,1); +count=0; +for i=1:nfun, + desvar=responsefunctions.desvar{i}; + if strncmpi(desvar,descriptor,descriptorlength), + count=count+1; + importancefactors(count)=responsefunctions.impfac(i); + end +end +if count==0, + error('plot_importancefactors error message: could not find to response functions with corresponding design variable'); +end + +%log? +if ~isnan(options_structure.log), + logvalue=options_structure.log; + importancefactors=log(importancefactors)/log(logvalue); +end + +%Ok, get partitioning. +[epart npart]=MeshPartition(md,md.npart); + +%distribute importance factor +gridimportance=importancefactors(npart); + +%process data and model +[x y z elements is2d]=processmesh(md,options_structure); + +%edgecolor? +if ~isnan(options_structure.edgecolor), + edgecolor=options_structure.edgecolor; +else + edgecolor='none'; +end + +%standard plot: +subplot(width,width,ii); + +%ok, plot gridimportance now. +if is2d, + A=elements(:,1); B=elements(:,2); C=elements(:,3); + patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', gridimportance,'FaceColor','interp','EdgeColor',edgecolor); +else + error('plot_importancefactors error message: 3d meshes not supported yet'); +end + +%apply options +applyoptions(md,[],options_structure); Index: /issm/trunk/src/m/classes/public/plot/plot_qmunormplot.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_qmunormplot.m (revision 345) +++ /issm/trunk/src/m/classes/public/plot/plot_qmunormplot.m (revision 345) @@ -0,0 +1,45 @@ +% +% plot a normal probability plot of the response functions. +% +% []=plot_normplot(rfunc) +% +function []=plot_qmunormplot(rfunc,width,ii) + +if ~nargin + help plot_normplot + return +end + +%% assemble the data into a matrix + +desc=cell (1,length(rfunc)); +for i=1:length(rfunc) + ldata(i)=length(rfunc(i).sample); +end +data=zeros(max(ldata),length(rfunc)); + +for i=1:length(rfunc) + desc(i)=cellstr(rfunc(i).descriptor); + data(1:ldata(i),i)=rfunc(i).sample; +end + +%standard plot: +subplot(width,width,ii); + +%% draw the plot + +% draw normal probability plot + +normplot(data) +ax1=gca; + +% add the annotation + +title('Normal Probability Plot of Design Variables and/or Response Functions') +xlabel('Value') +ylabel('Probability') + +hleg1=legend(ax1,desc,'Location','EastOutside',... + 'Orientation','vertical','Interpreter','none'); + +end Index: /issm/trunk/src/m/classes/public/plot/plotdoc.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plotdoc.m (revision 344) +++ /issm/trunk/src/m/classes/public/plot/plotdoc.m (revision 345) @@ -42,4 +42,6 @@ disp(' - ''thermaltransient_results'': this will display all the time steps of a thermal transient run'); disp(' - ''importancefactors'': qmu results'); +disp(' - ''normplot'': norm plot'); +disp(' - ''histnorm'': histogram normal distribution'); answer = lower(input(['more?(y/n) \n'],'s'));