Index: /issm/trunk/src/m/solutions/dakota/dakota_cdfs.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/dakota_cdfs.m (revision 3092) +++ /issm/trunk/src/m/solutions/dakota/dakota_cdfs.m (revision 3092) @@ -0,0 +1,330 @@ +% +% calculate the same cumulative distribution functions that dakota +% calculates for given responses, probabilities, reliabilities, and/or +% general reliabilities. +% +% [dresp]=dakota_cdfs(method,dresp ,resp,prob,rel,grel) +% [cdf ]=dakota_cdfs(method,samp ,resp,prob,rel,grel) +% [cdf ]=dakota_cdfs(method,mean,stddev,resp,prob,rel,grel) +% +% the required input is: +% method (char, 'nond_sampling' or 'nond_local_reliability') +% dresp (structure array, responses) +% or +% samp (double array, lists of samples) +% or +% mean (double vector, means) +% stddev (double vector, standard deviations) +% resp (double vector, list of responses) +% prob (double vector, list of probabilities) +% rel (double vector, list of reliabilities) +% grel (double vector, list of general reliabilities) +% +% and the optional input is: +% alpha (numeric, confidence interval of 100(1-alpha)%) +% +% the required field of dresp is (for nond_sampling): +% sample (double vector, list of samples) +% or (for nond_local_reliability): +% mean (double, mean of samples) +% stddev (double, standard deviation of samples) +% +% the required output is: +% dresp (structure array, responses) +% or +% cdf(:,4) (double, array of resp/prob/rel/grel) +% +% and the output fields of dresp are: +% cdf(:,4) (double, array of resp/prob/rel/grel) +% +% for each response (or column of data) in the input array, this +% function calculates the responses, probabilities, reliabilities +% and general reliabilities for the cumulative distribution function, +% the same way as dakota would. if the input is a structure, the +% output is a field in the structure; if the input is arrays, the +% output is an array. +% +% dresp data would typically be contained in the dakota tabular +% output file from a sampling analysis, or in the dakota output file +% from a local reliability analysis, either read by dakota_out_parse. +% +% "Copyright 2009, by the California Institute of Technology. +% ALL RIGHTS RESERVED. United States Government Sponsorship +% acknowledged. Any commercial use must be negotiated with +% the Office of Technology Transfer at the California Institute +% of Technology. (NTR 47078) +% +% This software may be subject to U.S. export control laws. +% By accepting this software, the user agrees to comply with +% all applicable U.S. export laws and regulations. User has the +% responsibility to obtain export licenses, or other export +% authority as may be required before exporting such information +% to foreign countries or providing access to foreign persons." +% +function [varargout]=dakota_cdfs(varargin) + +if ~nargin + help dakota_cdfs + return +end + +%% process input data + +iarg=1; +if ischar(varargin{iarg}) + method=varargin{iarg}; + iarg=iarg+1; + if ~strncmpi(method,'nond_s',6) && ~strncmpi(method,'nond_l',6) + error(['Method ''' method ''' is unrecognized.']); + end +else + method=''; +end + +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; +else + if strncmpi(method,'nond_s',6) + samp=varargin{iarg}; + iarg=iarg+1; + elseif strncmpi(method,'nond_l',6) + mean =varargin{iarg}; + iarg=iarg+1; + stddev=varargin{iarg}; + iarg=iarg+1; + end +end + +if iarg <= nargin && isnumeric(varargin{iarg}) + resp=varargin{iarg}; + iarg=iarg+1; +else + resp=[]; +end +if iarg <= nargin && isnumeric(varargin{iarg}) + prob=varargin{iarg}; + iarg=iarg+1; +else + prob=[]; +end +if iarg <= nargin && isnumeric(varargin{iarg}) + rel =varargin{iarg}; + iarg=iarg+1; +else + rel =[]; +end +if iarg <= nargin && isnumeric(varargin{iarg}) + grel=varargin{iarg}; + iarg=iarg+1; +else + grel=[]; +end + +%% calculate the cumulative distribution functions by input type + +if exist('dresp','var') && ~isempty(dresp) + if strncmpi(method,'nond_s',6) + for i=1:length(dresp) + [dresp(i).cdf]=cdfs_samp_calc(dresp(i).sample,... + resp,prob,rel,grel); + end + elseif strncmpi(method,'nond_l',6) + for i=1:length(dresp) + [dresp(i).cdf]=cdfs_lr_calc(dresp(i).mean,dresp(i).stddev,... + resp,prob,rel,grel); + end + end + + varargout{1}=dresp; + +elseif exist('samp','var') && ~isempty(samp) + cdf=zeros(length(resp)+length(prob)+length(rel)+length(grel),... + 4,size(samp,2)); + + for i=1:size(samp,2) + [cdf(:,:,i)]=cdfs_samp_calc(samp(:,i),... + resp,prob,rel,grel); + end + + varargout{1}=cdf; + +elseif exist('mean','var' ) && ~isempty(mean ) && ... + exist('stddev','var') && ~isempty(stddev) + cdf=zeros(length(resp)+length(prob)+length(rel)+length(grel),... + 4,length(mean)); + + for i=1:length(mean) + [cdf(:,:,i)]=cdfs_lr_calc(mean(i),stddev(i),... + resp,prob,rel,grel); + end + + varargout{1}=cdf; +else + error(['Empty data ''' inputname(2) ''' of type ''' class(varargin{2}) '''.']); +end + +end + +%% function to calculate the results for a sampling analysis + +function [cdf]=cdfs_samp_calc(samp,resp,prob,rel,grel) + +% sort the samples and remove any NaN padding (should only occur at end) + + samp=sort(samp(~isnan(samp(:))),'ascend'); + nsamp=length(samp); + + mu =mean(samp); + sigma=std(samp); + + cdf=zeros(length(resp)+length(prob)+length(rel)+length(grel),4); + cdf(:,:)=NaN; + irow=0; + +% compute quantities, given response levels + + for i=1:length(resp) + irow=irow+1; + indx=bin_search_val(resp(i),samp); + cdf(irow,1)=resp(i); + cdf(irow,2)=indx/nsamp; + cdf(irow,3)=(mu-resp(i))/sigma; +% cdf(irow,4)=-sqrt(2)*erfinv((indx-nsamp/2)/(nsamp/2)); + cdf(irow,4)=sqrt(2)*erfcinv(indx/(nsamp/2)); + end + +% compute response levels, given probabilities + + for i=1:length(prob) + irow=irow+1; +% why not round(prob(i)*(nsamp-1)+1)? + indx=ceil(prob(i)*(nsamp)); + if (indx < 1) + indx=1; + elseif (indx > nsamp) + indx=nsamp; + end + cdf(irow,1)=samp(indx); + cdf(irow,2)=prob(i); + end + +% compute response levels, given reliabilities + + for i=1:length(rel) + irow=irow+1; + cdf(irow,1)=mu-sigma*rel(i); + cdf(irow,3)=rel(i); + end + +% compute response levels, given general reliabilities + + for i=1:length(grel) + irow=irow+1; +% indx=ceil(nsamp/2+nsamp/2*erf(-grel(i)/sqrt(2))); + indx=ceil((nsamp/2)*erfc(grel(i)/sqrt(2))); + if (indx < 1) + indx=1; + elseif (indx > nsamp) + indx=nsamp; + end + cdf(irow,1)=samp(indx); + cdf(irow,4)=grel(i); + end + +end + +%% function to calculate the results for a local reliability analysis + +function [cdf]=cdfs_lr_calc(mu,sigma,resp,prob,rel,grel) + + cdf=zeros(length(resp)+length(prob)+length(rel)+length(grel),4); + irow=0; + +% compute quantities, given response levels + + for i=1:length(resp) + irow=irow+1; + cdf(irow,1)=resp(i); + cdf(irow,2)=normcdf(resp(i),mu,sigma); + cdf(irow,3)=(mu-resp(i))/sigma; + cdf(irow,4)=(mu-resp(i))/sigma; + end + +% compute quantities, given probabilities + + for i=1:length(prob) + irow=irow+1; + cdf(irow,1)=norminv(prob(i),mu,sigma); + cdf(irow,2)=prob(i); + cdf(irow,3)=-norminv(prob(i),0,1); + cdf(irow,4)=-norminv(prob(i),0,1); + end + +% compute quantities, given reliabilities + + for i=1:length(rel) + irow=irow+1; + cdf(irow,1)=mu-sigma*rel(i); + cdf(irow,2)=normcdf(-rel(i),0,1); + cdf(irow,3)=rel(i); + cdf(irow,4)=rel(i); + end + +% compute quantities, given general reliabilities + + for i=1:length(grel) + irow=irow+1; + cdf(irow,1)=mu-sigma*grel(i); + cdf(irow,2)=normcdf(-grel(i),0,1); + cdf(irow,3)=grel(i); + cdf(irow,4)=grel(i); + end + +end +%% +% function to perform a recursive binary search for a matrix of values +% in an ordered vector (loop separately outside of recursion for +% efficiency purposes) +% +% function [ind]=bin_search(val,vect) +% +function [ind]=bin_search(val,vect) + +ind=zeros(size(val)); + +for i=1:numel(val) + ind(i)=bin_search_val(val(i),vect); +end + +end +%% +% function to perform a recursive binary search in an ordered vector, +% returning low if value does not exist (more efficient than find or +% ismember, which must use linear searches and/or sort) +% +% function [ind]=bin_search_val(val,vect) +% +function [ind]=bin_search_val(val,vect) + +imid=floor((1+length(vect))/2); + +if (val == vect(imid)) + ind=imid; +elseif (val < vect(imid)) + if (imid > 1) + ind= bin_search(val,vect(1:imid-1)); + else + ind=0; + end +elseif (val > vect(imid)) + if (imid < length(vect)) + ind=imid+bin_search(val,vect(imid+1:length(vect))); + else + ind=length(vect); + end +else + ind=NaN; +end + +end Index: /issm/trunk/src/m/solutions/dakota/dakota_in_params.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/dakota_in_params.m (revision 3091) +++ /issm/trunk/src/m/solutions/dakota/dakota_in_params.m (revision 3092) @@ -26,5 +26,5 @@ % acknowledged. Any commercial use must be negotiated with % the Office of Technology Transfer at the California Institute -% of Technology. (NTR 47078) +% of Technology. (J. Schiermeier, NTR 47078) % % This software may be subject to U.S. export control laws. @@ -66,4 +66,7 @@ % or result from the response level lists +if ~isfield(params,'compute') + params.compute='probabilities'; +end if ~isfield(params,'distribution') params.distribution='cumulative'; @@ -172,5 +175,5 @@ if ~isfield(params,'numerical_gradients') - params.numerical_gradients=true; + params.numerical_gradients=false; end if ~isfield(params,'method_source') @@ -198,5 +201,5 @@ % hessians not fully implemented if ~isfield(params,'numerical_hessians') - params.numerical_gradients=true; + params.numerical_hessians=true; end if ~isfield(params,'hessian_gradient_step_size') Index: /issm/trunk/src/m/solutions/dakota/dakota_in_write.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/dakota_in_write.m (revision 3091) +++ /issm/trunk/src/m/solutions/dakota/dakota_in_write.m (revision 3092) @@ -2,5 +2,6 @@ % write a Dakota .in input file. % -% []=dakota_in_write(method,dmeth,dvar,dresp,params,filei,varargin) +% []=dakota_in_write(method,dvar,dresp,params,filei,varargin) +% []=dakota_in_write(dmeth ,dvar,dresp,params,filei,varargin) % % where the required input is: @@ -12,6 +13,6 @@ % filei (character, name of .in file) % -% the method, dmeth, and filei will be prompted if empty. -% params may be empty, in which case defaults will be used. +% the method and filei will be prompted if empty. params +% may be empty, in which case defaults will be used. % % the optional varargin are not yet used. @@ -30,5 +31,5 @@ % acknowledged. Any commercial use must be negotiated with % the Office of Technology Transfer at the California Institute -% of Technology. (NTR 47078) +% of Technology. (J. Schiermeier, NTR 47078) % % This software may be subject to U.S. export control laws. @@ -39,5 +40,5 @@ % to foreign countries or providing access to foreign persons." % -function []=dakota_in_write(method,dmeth,dvar,dresp,params,filei,varargin) +function []=dakota_in_write(method,dvar,dresp,params,filei,varargin) if ~nargin @@ -51,7 +52,10 @@ method=input('Method? ','s'); end - -if ~exist('dmeth' ,'var') || isempty(dmeth) +if ischar(method) dmeth=dakota_method(method); +elseif isa(method,'dakota_method') + dmeth=method; +else + error(['Method ''' inputname(1) ''' is unrecognized class ''' class(method) '''.']); end @@ -134,42 +138,11 @@ % write response levels -if strcmp(dmeth.type,'nond') && isfield(dresp,'rf') - param_write(fidi,'\t ','distribution',' ','\n',params); - [respl,probl,rell,grell]=prop_levels(dresp.rf); - if ~isempty(respl) - rlev_write(fidi,'response_levels',respl); - param_write(fidi,'\t ','compute',' ','\n',params); - end - if ~isempty(probl) - rlev_write(fidi,'probability_levels',probl); - end - if ~isempty(rell) - rlev_write(fidi,'reliability_levels',rell); - end - if ~isempty(grell) - rlev_write(fidi,'gen_reliability_levels',grell); - end -end -fprintf(fidi,'\n'); - -end - -%% function to write response levels - -function []=rlev_write(fidi,ltype,levels) - -fprintf(fidi,'\t num_%s =',ltype); -for i=1:length(levels) - fprintf(fidi,' %d',length(levels{i})); -end -fprintf(fidi,'\n'); - -fprintf(fidi,'\t %s =\n',ltype); - -for i=1:length(levels) - if ~isempty(levels{i}) - vector_write(fidi,sprintf('\t '),levels{i},8,76); - end -end +if strcmp(dmeth.type,'nond') + for i=1:length(dmeth.responses) + fhresp=str2func([dmeth.responses{i} '.dakota_rlev_write']); + fhresp(fidi,dresp,params); + end +end +fprintf(fidi,'\n'); end @@ -196,165 +169,17 @@ % variables vary by method -vsets_write(fidi,dvar,dmeth.variables); +for i=1:length(dmeth.variables) + fhvar=str2func([dmeth.variables{i} '.dakota_write']); + fhvar(fidi,dvar); +end % linear constraints vary by method -lcsets_write(fidi,dvar,dmeth.lcspec); - -fprintf(fidi,'\n'); - -end - -%% function to write variable sets - -function []=vsets_write(fidi,dvar,variables) - -for i=1:length(variables) - switch(variables{i}) - case 'cdv' - cstring='continuous_design'; - case 'nuv' - cstring='normal_uncertain'; - case 'csv' - cstring='continuous_state'; - otherwise - warning('vsets_write:unrec_var',... - 'Unrecognized variable ''%s''.',variables{i}); - end - - if isfield(dvar,variables{i}) - vlist_write(fidi,cstring,variables{i},dvar.(variables{i})); - end -end - -end - -%% function to write variable list - -function []=vlist_write(fidi,cstring,cstring2,dvar) - -% put variables into lists for writing - -pinitpt=prop_initpt(dvar); -plower =prop_lower (dvar); -pupper =prop_upper (dvar); -pmean =prop_mean (dvar); -pstddev=prop_stddev(dvar); -pinitst=prop_initst(dvar); -pstype =prop_stype (dvar); -pscale =prop_scale (dvar); -pdesc =prop_desc (dvar); - -% write variables -% (using Dakota 4.1 syntax for backward compatability) - -disp(sprintf(' Writing %d %s variables.',length(dvar),class(dvar))); - -fprintf(fidi,'\t%s = %d\n',cstring,length(dvar)); -if ~isempty(pinitpt) - fprintf(fidi,'\t %s_initial_point =\n',cstring2); - vector_write(fidi,sprintf('\t '),pinitpt,6,76); -end -if ~isempty(plower) - fprintf(fidi,'\t %s_lower_bounds =\n',cstring2); - vector_write(fidi,sprintf('\t '),plower ,6,76); -end -if ~isempty(pupper) - fprintf(fidi,'\t %s_upper_bounds =\n',cstring2); - vector_write(fidi,sprintf('\t '),pupper ,6,76); -end -if ~isempty(pmean) - fprintf(fidi,'\t %s_means =\n',cstring2); - vector_write(fidi,sprintf('\t '),pmean ,6,76); -end -if ~isempty(pstddev) - fprintf(fidi,'\t %s_std_deviations =\n',cstring2); - vector_write(fidi,sprintf('\t '),pstddev,6,76); -end -if ~isempty(pinitst) - fprintf(fidi,'\t %s_initial_state =\n',cstring2); - vector_write(fidi,sprintf('\t '),pinitst,6,76); -end -if ~isempty(pstype) - fprintf(fidi,'\t %s_scale_types =\n',cstring2); - vector_write(fidi,sprintf('\t '),pstype ,6,76); -end -if ~isempty(pscale) - fprintf(fidi,'\t %s_scales =\n',cstring2); - vector_write(fidi,sprintf('\t '),pscale ,6,76); -end -if ~isempty(pdesc) - fprintf(fidi,'\t %s_descriptors =\n',cstring2); - vector_write(fidi,sprintf('\t '),pdesc ,6,76); -end - -end - -%% function to write linear constraint sets - -function []=lcsets_write(fidi,dvar,lcspec) - -for i=1:length(lcspec) - switch(lcspec{i}) - case 'lic' - cstring ='linear_inequality_constraints'; - cstring2='linear_inequality'; - case 'lec' - cstring ='linear_equality_constraints'; - cstring2='linear_equality'; - otherwise - warning('lcspec_write:unrec_lcspec',... - 'Unrecognized linear constraint ''%s''.',lcspec{i}); - end - - if isfield(dvar,lcspec{i}) - lclist_write(fidi,cstring,cstring2,dvar.(lcspec{i})); - end -end - -end - -%% function to write linear constraint list - -function []=lclist_write(fidi,cstring,cstring2,dvar) - -% put linear constraints into lists for writing - -pmatrix=prop_matrix(dvar); -plower =prop_lower (dvar); -pupper =prop_upper (dvar); -ptarget=prop_target(dvar); -pstype =prop_stype (dvar); -pscale =prop_scale (dvar); - -% write linear constraints - -disp(sprintf(' Writing %d %s linear constraints.',... - length(dvar),class(dvar))); - -if ~isempty(pmatrix) - fprintf(fidi,'\t %s_matrix =\n',cstring2); - vector_write(fidi,sprintf('\t '),pmatrix,6,76); -end -if ~isempty(plower) - fprintf(fidi,'\t %s_lower_bounds =\n',cstring2); - vector_write(fidi,sprintf('\t '),plower ,6,76); -end -if ~isempty(pupper) - fprintf(fidi,'\t %s_upper_bounds =\n',cstring2); - vector_write(fidi,sprintf('\t '),pupper ,6,76); -end -if ~isempty(ptarget) - fprintf(fidi,'\t %s_targets =\n',cstring2); - vector_write(fidi,sprintf('\t '),ptarget,6,76); -end -if ~isempty(pstype) - fprintf(fidi,'\t %s_scale_types =\n',cstring2); - vector_write(fidi,sprintf('\t '),pstype ,6,76); -end -if ~isempty(pscale) - fprintf(fidi,'\t %s_scales =\n',cstring2); - vector_write(fidi,sprintf('\t '),pscale ,6,76); -end +for i=1:length(dmeth.lcspec) + fhvar=str2func([dmeth.lcspec{i} '.dakota_write']); + fhvar(fidi,dvar); +end + +fprintf(fidi,'\n'); end @@ -407,10 +232,12 @@ param_write(fidi,'\t','direct','','\n',params); param_write(fidi,'\t ','analysis_driver',' = ''','''\n',params); - [pathstr,name,ext,versn] = fileparts(params.analysis_components); - if isempty(ext) - ext='.m'; - end - params.analysis_components=fullfile(pathstr,[name ext versn]); - param_write(fidi,'\t ','analysis_components',' = ''','''\n',params); + if ~isempty(params.analysis_components) + [pathstr,name,ext,versn] = fileparts(params.analysis_components); + if isempty(ext) + ext='.m'; + end + params.analysis_components=fullfile(pathstr,[name ext versn]); + param_write(fidi,'\t ','analysis_components',' = ''','''\n',params); + end if ~isempty(params.input_filter) param_write(fidi,'\t ','input_filter',' = ''','''\n',params); @@ -421,5 +248,4 @@ param_write(fidi,'\t ','failure_capture',' ','\n',params); param_write(fidi,'\t ','deactivate',' ','\n',params); - param_write(fidi,'\t ','output',' ','\n',params); param_write(fidi,'\t ','processors_per_analysis',' = ''','''\n',params); end @@ -439,45 +265,12 @@ % functions, gradients, and hessians vary by method -rsets_write(fidi,dresp,dmeth.responses); -ghspec_write(fidi,params,dmeth.ghspec); - -fprintf(fidi,'\n'); - -end - -%% function to write response sets - -function []=rsets_write(fidi,dresp,responses) - rdesc={}; -for i=1:length(responses) - switch(responses{i}) - case 'of' - cstring ='objective_functions'; - cstring2='objective_function'; - case 'lst' - cstring ='least_squares_terms'; - cstring2='least_squares_term'; - case 'nic' - cstring ='nonlinear_inequality_constraints'; - cstring2='nonlinear_inequality'; - case 'nec' - cstring ='nonlinear_equality_constraints'; - cstring2='nonlinear_equality'; - case 'rf' - cstring ='response_functions'; - cstring2='response_function'; - otherwise - warning('rsets_write:unrec_resp',... - 'Unrecognized response ''%s''.',responses{i}); - end - - if isfield(dresp,responses{i}) - [rdesc]=rlist_write(fidi,cstring,cstring2,dresp.(responses{i}),rdesc); - end -end - -% write response descriptors +for i=1:length(dmeth.responses) + fhresp=str2func([dmeth.responses{i} '.dakota_write']); + [rdesc]=fhresp(fidi,dresp,rdesc); +end + +% write accumulated response descriptors for all response classes if ~isempty(rdesc) @@ -486,58 +279,7 @@ end -end - -%% function to write response list - -function [rdesc]=rlist_write(fidi,cstring,cstring2,dresp,rdesc) - -% put responses into lists for writing - -pstype =prop_stype (dresp); -pscale =prop_scale (dresp); -pweight=prop_weight(dresp); -plower =prop_lower (dresp); -pupper =prop_upper (dresp); -ptarget=prop_target(dresp); - -% accumulate descriptors into list for subsequent writing - -rdesc(end+1:end+numel(dresp))=prop_desc (dresp); - -% write responses - -disp(sprintf(' Writing %d %s responses.',length(dresp),class(dresp))); - -fprintf(fidi,'\tnum_%s = %d\n',cstring,length(dresp)); -if ~isempty(pstype) - fprintf(fidi,'\t %s_scale_types =\n',cstring2); - vector_write(fidi,sprintf('\t '),pstype ,6,76); -end -if ~isempty(pscale) - fprintf(fidi,'\t %s_scales =\n',cstring2); - vector_write(fidi,sprintf('\t '),pscale ,6,76); -end -if ~isempty(pweight) - switch cstring2 - case 'objective_function' - fprintf(fidi,'\t %s_weights =\n','multi_objective'); - vector_write(fidi,sprintf('\t '),pweight,6,76); - case 'least_squares_term' - fprintf(fidi,'\t %s_weights =\n','least_squares'); - vector_write(fidi,sprintf('\t '),pweight,6,76); - end -end -if ~isempty(plower) - fprintf(fidi,'\t %s_lower_bounds =\n',cstring2); - vector_write(fidi,sprintf('\t '),plower ,6,76); -end -if ~isempty(pupper) - fprintf(fidi,'\t %s_upper_bounds =\n',cstring2); - vector_write(fidi,sprintf('\t '),pupper ,6,76); -end -if ~isempty(ptarget) - fprintf(fidi,'\t %s_targets =\n',cstring2); - vector_write(fidi,sprintf('\t '),ptarget,6,76); -end +ghspec_write(fidi,params,dmeth.ghspec); + +fprintf(fidi,'\n'); end @@ -550,4 +292,10 @@ if find_string(ghspec,'grad') + if ~params.numerical_gradients && ~params.analytic_gradients + params.numerical_gradients=true; + elseif (params.numerical_gradients+params.analytic_gradients > 1) + error('Too many gradients selected.') + end + if params.numerical_gradients param_write(fidi,'\t','numerical_gradients','','\n',params); @@ -572,78 +320,2 @@ end - -%% function to write a parameter - -function []=param_write(fidi,sbeg,pname,smid,send,params) - -if ~isfield(params,pname) - warning('param_write:param_not_found',... - 'Parameter ''%s'' not found in structure.',pname); - return -end - -if islogical(params.(pname)) && ~params.(pname) - return -end - -if islogical(params.(pname)) - fprintf(fidi,[sbeg '%s' send],pname); -elseif ischar (params.(pname)) - fprintf(fidi,[sbeg '%s' smid '%s' send],pname,params.(pname)); -elseif isnumeric(params.(pname)) - fprintf(fidi,[sbeg '%s' smid '%g' send],pname,params.(pname)); -end - -end - -%% function to write a vector on multiple lines - -function []=vector_write(fidi,sbeg,vec,nmax,cmax) - -if ~exist('nmax','var') || isempty(nmax) - nmax=Inf; -end -if ~exist('cmax','var') || isempty(cmax) - cmax=Inf; -end - -% set up first iteration - -svec =[]; -nitem=nmax; -lsvec=cmax; - -% transpose vector from column-wise to row-wise - -vec=vec'; - -% assemble each line, flushing when necessary - -for i=1:numel(vec) - if isnumeric(vec(i)) - sitem=sprintf('%g' ,vec(i)); - else - sitem=sprintf('''%s''',char(vec(i))); - end - nitem=nitem+1; - lsvec=lsvec+1+length(sitem); - - if (nitem <= nmax) && (lsvec <= cmax) - svec=[svec ' ' sitem]; - else - if ~isempty(svec) - fprintf(fidi,'%s\n',svec); - end - svec=[sbeg sitem]; - nitem=1; - lsvec=length(svec); - end -end - -% flush buffer at end, if necessary - -if ~isempty(svec) - fprintf(fidi,'%s\n',svec); -end - -end Index: /issm/trunk/src/m/solutions/dakota/dakota_m_write.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/dakota_m_write.m (revision 3092) +++ /issm/trunk/src/m/solutions/dakota/dakota_m_write.m (revision 3092) @@ -0,0 +1,296 @@ +% +% write a Matlab .m function file to be called by Dakota for +% the Matlab direct or external driver. +% +% []=dakota_m_write(method,dmeth,dvar,dresp,params,filem,package,varargin) +% +% where the required input is: +% method (character, dakota method name) +% dmeth (dakota_method, method class object) +% dvar (structure array, variable class objects) +% dresp (structure array, response class objects) +% params (structure array, method-independent parameters) +% filem (character, name of .m file) +% package (character, analysis package) +% +% the method, dmeth, and filem will be prompted if empty. +% params may be empty, in which case defaults will be used. +% +% the optional varargin are passed directly through to the +% QmuUpdateFunctions brancher to be used by the analysis +% package. for example, this could be model information. +% +% this function writes a matlab .m function file to be called +% by dakota for the matlab direct or external driver. for +% the direct driver, dakota is linked with matlab and +% automatically starts a matlab session, passing the variables +% and responses through the function; for the external driver, +% dakota calls a shell script to start the matlab session, +% passing the variables and responses through text files. +% this function must be tailored to the particular analysis +% package. +% +% this data would typically be generated by a matlab script +% for a specific model, using the method, variable, and +% response class objects. this function may be called by +% dakota_in_data. +% +% "Copyright 2009, by the California Institute of Technology. +% ALL RIGHTS RESERVED. United States Government Sponsorship +% acknowledged. Any commercial use must be negotiated with +% the Office of Technology Transfer at the California Institute +% of Technology. (NTR 47078) +% +% This software may be subject to U.S. export control laws. +% By accepting this software, the user agrees to comply with +% all applicable U.S. export laws and regulations. User has the +% responsibility to obtain export licenses, or other export +% authority as may be required before exporting such information +% to foreign countries or providing access to foreign persons." +% +function []=dakota_m_write(method,dmeth,dvar,dresp,params,filem,package,varargin) + +if ~nargin + help dakota_m_write + return +end + +% process the input parameters + +if ~exist('method','var') || isempty(method) + method=input('Method? ','s'); +end + +if ~exist('dmeth' ,'var') || isempty(dmeth) + dmeth=dakota_method(method); +end + +if ~exist('params','var') + params=[]; +end + +if ~exist('filem' ,'var') || isempty(filem) + filem=input('Matlab m-file to write? ','s'); +end +[pathstr,name,ext,versn] = fileparts(filem); +if isempty(ext) + ext='.m'; +end +filem2=fullfile(pathstr,[name ext versn]); + +display(sprintf('Opening Matlab m-file ''%s''.',filem2)); +fidm=fopen(sprintf('%s',filem2),'w'); +if (fidm < 0) + error('''%s'' could not be opened.',filem2); +end + +% write beginning of the Matlab m-file + +begin_write(fidm,name,params); + +% write variables into the Matlab m-file + +variables_write(fidm,dmeth,dvar,params,varargin{:}); + +% write solution into the Matlab m-file + +solution_write(fidm,package); + +% write responses into the Matlab m-file + +responses_write(fidm,dmeth,params,dresp); + +% write end of the Matlab m-file + +end_write(fidm,name,params); + +fclose(fidm); +display('End of file successfully written.'); + +end + +%% function to write beginning of the Matlab m-file + +function []=begin_write(fidm,name,params) + +display('Writing beginning of Matlab m-file.'); + +fprintf(fidm,'%%\n'); +if strcmpi(params.analysis_driver,'matlab') + fprintf(fidm,'%% Run the specified input variables and return the output responses.\n'); + fprintf(fidm,'%%\n'); + fprintf(fidm,'function Dakota=%s(Dakota)\n\n',name); + fprintf(fidm,'clk=clock;\n'); + fprintf(fidm,'cpu=cputime;\n\n'); + fprintf(fidm,'%% Dakota\n\n'); + fprintf(fidm,'infile=''%s'';\n','Qmu.model'); +else + fprintf(fidm,'%% Run the specified input file and return the output file.\n'); + fprintf(fidm,'%%\n'); + fprintf(fidm,'function %s(infile,outfile)\n\n',name); + fprintf(fidm,'clk=clock;\n'); + fprintf(fidm,'cpu=cputime;\n\n'); +end +fprintf(fidm,'try\n'); +fprintf(fidm,'\tloadmodel(infile);\n\n'); + +if strcmpi(params.analysis_driver,'matlab') + fprintf(fidm,'\tmd=qmuname(md,Dakota.fnEvalId);\n\n'); +else + fprintf(fidm,'\tmd=qmuname(md);\n\n'); +end + +end + +%% function to write variables into the Matlab m-file + +function []=variables_write(fidm,dmeth,dvar,params,varargin) + +display('Writing variables for Matlab m-file.'); + +fprintf(fidm,'%% Apply the variables.\n\n'); +ixc=0; + +% variables vary by method + +ixc=vsets_write(fidm,ixc,dvar,dmeth.variables,params,varargin{:}); + +end + +%% function to write variable sets into the Matlab m-file + +function [ixc]=vsets_write(fidm,ixc,dvar,variables,params,varargin) + +for i=1:length(variables) + if isfield(dvar,variables{i}) + ixc=vlist_write(fidm,ixc,variables{i},dvar.(variables{i}),params,varargin{:}); + end +end + +end + +%% function to write variable list into the Matlab m-file + +function [ixc]=vlist_write(fidm,ixc,vtype,dvar,params,varargin) + +disp(sprintf(' Writing %d %s variables.',length(dvar),class(dvar))); + +for i=1:length(dvar) + + %first find descriptor, without sample number on it. + descriptor=discardnum(dvar(i).descriptor); + + %is there a lock on this variable? We don't want to update the variable twice! + if exist([descriptor '_lock'],'var'), + %lock is in place, do not update! continue; + continue; + else + %first things first, put lock in place + eval([descriptor '_lock=1;']); + + %now, we need a string to put in the matlab file, which will update all the variables + %for this descriptor. + [string,ixc]=QmuUpdateFunctions(ixc,descriptor,dvar,params,i,varargin{:}); + + %dump this string in the matlab file. + fprintf(fidm,'%s',string); + end +end + +end + +%% function to write solution into the Matlab m-file + +function []=solution_write(fidm,package) + +display('Writing solution for Matlab m-file.'); +fprintf(fidm,'%% Run the solution.\n\n'); + +fprintf(fidm,['\tmd=solve(md,''diagnostic'',''' package ''');\n\n']); + +end + +%% function to write responses into the Matlab m-file + +function []=responses_write(fidm,dmeth,params,dresp) + +display('Writing responses for Matlab m-file.'); + +fprintf(fidm,'%% Calculate the responses.\n\n'); +ifnvals=0; + +if ~strcmpi(params.analysis_driver,'matlab') + fprintf(fidm,'\tfid=fopen(outfile,''w'');\n\n'); +end + +% responses vary by method + +ifnvals=rsets_write(fidm,ifnvals,params,dresp,dmeth.responses); + +fprintf(fidm,'\n'); +if ~strcmpi(params.analysis_driver,'matlab') + fprintf(fidm,'\tstatus=fclose(fid);\n\n'); +end + +end + +%% function to write response sets into the Matlab m-file + +function [ifnvals]=rsets_write(fidm,ifnvals,params,dresp,responses) + +for i=1:length(responses) + if isfield(dresp,responses{i}) + ifnvals=rlist_write(fidm,ifnvals,params,responses{i},dresp.(responses{i})); + end +end + +end + +%% function to write response list into the Matlab m-file + +function [ifnvals]=rlist_write(fidm,ifnvals,params,rtype,dresp) + +disp(sprintf(' Writing %d %s responses.',length(dresp),class(dresp))); +for i=1:length(dresp) + ifnvals=ifnvals+1; + if strcmpi(params.analysis_driver,'matlab') + fprintf(fidm,'\tDakota.fnVals(%d)=QmuResponseValue(md,''%s'');\n',ifnvals,dresp(i).descriptor); + else + fprintf(fidm,'\tfprintf(fid,''%%f\\n'',QmuResponseValue(md,''%s''));\n',dresp(i).descriptor); + end +end + +end + +%% function to write end of the Matlab m-file + +function []=end_write(fidm,name,params) + +display('Writing end of Matlab m-file.'); + +fprintf(fidm,'%% Error condition.\n\n'); + +fprintf(fidm,'catch ME\n'); +fprintf(fidm,'\tME\n'); +fprintf(fidm,'\tfor i=1:length(ME.stack)\n'); +fprintf(fidm,'\t\tdisplay(sprintf('' file(%%d): %%s'', i,ME.stack(i).file));\n'); +fprintf(fidm,'\t\tdisplay(sprintf('' name(%%d): %%s'', i,ME.stack(i).name));\n'); +fprintf(fidm,'\t\tdisplay(sprintf('' line(%%d): %%d\\n'',i,ME.stack(i).line));\n'); +fprintf(fidm,'\tend\n'); +if strcmpi(params.analysis_driver,'matlab') + fprintf(fidm,'\tDakota.failure=1;\n'); +else + fprintf(fidm,'\tif exist(''fid'',''var'')\n'); + fprintf(fidm,'\t\tstatus=fclose(fid);\n'); + fprintf(fidm,'\tend\n'); + fprintf(fidm,'\tfid=fopen(outfile,''w'');\n'); + fprintf(fidm,'\tfprintf(fid,''fail\\n'');\n'); + fprintf(fidm,'\tstatus=fclose(fid);\n'); +end +fprintf(fidm,'end\n\n'); + +fprintf(fidm,'disp(sprintf(''%s -- %%f CPU seconds; %%f clock seconds\\n'',...\n',name); +fprintf(fidm,' cputime-cpu,etime(clock,clk)))\n\n'); +fprintf(fidm,'end\n\n'); + +end Index: /issm/trunk/src/m/solutions/dakota/dakota_moments.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/dakota_moments.m (revision 3092) +++ /issm/trunk/src/m/solutions/dakota/dakota_moments.m (revision 3092) @@ -0,0 +1,140 @@ +% +% calculate the same moments and confidence intervals that dakota +% calculates for a sample. +% +% [dresp ]=dakota_moments(dresp,alpha) +% [mean,stddev,meanci,stddevci]=dakota_moments(samp ,alpha) +% +% the required input is: +% dresp (structure array, responses) +% or +% samp (double array, lists of samples) +% +% and the optional input is: +% alpha (numeric, confidence interval of 100(1-alpha)%) +% +% the required field of dresp is: +% sample (double vector, list of samples) +% +% the required output is: +% dresp (structure array, responses) +% or +% mean (double, mean of sample) +% stddev (double, standard deviation of sample) +% meanci(2) (double, confidence interval of mean) +% stddevci(2) (double, confidence interval of standard deviation) +% +% and the output fields of dresp are: +% mean (double, mean of sample) +% stddev (double, standard deviation of sample) +% meanci(2) (double, confidence interval of mean) +% stddevci(2) (double, confidence interval of standard deviation) +% +% for each response (or column of data) in the input array, this +% function calculates the mean, standard deviation, and their +% confidence intervals for a normal distribution, the same way as +% dakota would. if the input is a structure, the output is fields +% in the structure; if the input is an array, the output is arrays. +% +% dresp data would typically be contained in the dakota tabular +% output file from a sampling analysis, read by dakota_out_parse. +% +% "Copyright 2009, by the California Institute of Technology. +% ALL RIGHTS RESERVED. United States Government Sponsorship +% acknowledged. Any commercial use must be negotiated with +% the Office of Technology Transfer at the California Institute +% of Technology. (NTR 47078) +% +% This software may be subject to U.S. export control laws. +% By accepting this software, the user agrees to comply with +% all applicable U.S. export laws and regulations. User has the +% responsibility to obtain export licenses, or other export +% authority as may be required before exporting such information +% to foreign countries or providing access to foreign persons." +% +function [varargout]=dakota_moments(varargin) + +if ~nargin + help dakota_moments + return +end + +%% process input data + +iarg=1; +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; +elseif isnumeric(varargin{iarg}) + samp=varargin{iarg}; + iarg=iarg+1; +else + error(['Unknown data ''' inputname(1) ''' of type ''' class(varargin{1}) '''.']); +end + +if iarg <= nargin && isnumeric(varargin{iarg}) + alpha=varargin{iarg}; + iarg=iarg+1; +else +% use dakota default of 95% + alpha=0.05; +end + +%% calculate the moments and confidence intervals by input type + +if exist('dresp','var') && ~isempty(dresp) + for i=1:length(dresp) + [dresp(i).mean,dresp(i).stddev,... + dresp(i).meanci,dresp(i).stddevci]=... + moments_calc(dresp(i).sample,alpha); + end + + varargout{1}=dresp; + +elseif exist('samp','var') && ~isempty(samp) + mean =zeros(1,size(samp,2)); + stddev =zeros(1,size(samp,2)); + meanci =zeros(2,size(samp,2)); + stddevci=zeros(2,size(samp,2)); + +% could do this using vector math rather than loop, but want to allow +% the case of differently sized samples padded by NaN's + + for i=1:size(samp,2) + [mean(i),stddev(i),... + meanci(:,i),stddevci(:,i)]=... + moments_calc(samp(:,i),alpha); + end + + varargout{1}=mean; + varargout{2}=stddev; + varargout{3}=meanci; + varargout{4}=stddevci; +else + error(['Empty data ''' inputname(1) ''' of type ''' class(varargin{1}) '''.']); +end + +end + +%% function to calculate the results + +function [mu,sigma,muci,sigmaci]=moments_calc(samp,alpha) + +% remove any NaN padding (should only occur at end) + + samp=samp(~isnan(samp(:))); + nsamp=length(samp); + prob=1-alpha/2; + +% could use Matlab normfit, but make calculations explicit +% [mu,sigma,muci,sigmaci]=normfit(samp,alpha); + + mu =mean(samp); + sigma=std(samp); + + muci(1) =mu-tinv(prob,nsamp-1)*sigma/sqrt(nsamp); + muci(2) =mu+tinv(prob,nsamp-1)*sigma/sqrt(nsamp); + sigmaci(1)=sigma*sqrt((nsamp-1)/chi2inv(prob ,nsamp-1)); + sigmaci(2)=sigma*sqrt((nsamp-1)/chi2inv(1-prob,nsamp-1)); + +end Index: /issm/trunk/src/m/solutions/dakota/dakota_out_parse.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/dakota_out_parse.m (revision 3091) +++ /issm/trunk/src/m/solutions/dakota/dakota_out_parse.m (revision 3092) @@ -88,5 +88,5 @@ [ntokens,tokens]=fltokens(fline); method=tokens{1}{3}; -display(sprintf('Dakota method=%s.',method)); +display(sprintf('Dakota methodName=''%s''.',method)); dresp=struct([]); @@ -96,12 +96,8 @@ prcm=struct([]); -%% loop through the file to find the iterator completion - -[fline]=findline(fidi,['<<<<< Iterator ' method ' completed.']); -if ~ischar(fline) - return -end -% display([' ' deblank(fline)]); - +%% loop through the file to find the function evaluation summary + +fline=''; +[nfeval]=nfeval_read(fidi,fline); fline=fgetl(fidi); @@ -110,5 +106,6 @@ while ischar(fline) % ipos=ftell(fidi); - if strncmp(fline,'<<<<< Function evaluation summary',33) + if isempty(fline) + elseif strncmp(fline,'<<<<< Function evaluation summary',33) [nfeval]=nfeval_read(fidi,fline); elseif strncmp(fline,'Statistics based on ',20) @@ -134,6 +131,9 @@ elseif strncmp(fline,'The following lists volumetric uniformity measures',50) [dresp]=vum_read(fidi,dresp,fline); + elseif strncmp(fline,'<<<<< Iterator ',15) && ... + (length(fline) > 26) && ... + ~isempty(strfind(fline(16:end),' completed.')) + [method]=itcomp_read(fidi,fline); elseif strncmp(fline,'-----',5) - elseif isempty(fline) else display(['Unexpected line: ' deblank(fline)]); @@ -678,4 +678,19 @@ end +% read and add best residual norms + + elseif strncmpi(cellstr(tokens{1}{3}),'residual', 8) && ... + strncmpi(cellstr(tokens{1}{4}),'norm' , 4) + display([' ' deblank(fline)]); + dresp.best.norm = tokens{1}{ 6}; + dresp.best.hnormsq= tokens{1}{11}; + + fline=fgetl(fidi); + + while ischar(fline) && ~isempty(fline) && ... + ~strncmpi(fline,'<<<<< Best ',11) + fline=fgetl(fidi); + end + % read and add best residual term(s) @@ -715,21 +730,22 @@ strncmpi(cellstr(tokens{1}{4}),'captured', 8) display([' ' deblank(fline)]); - [ntokens,tokens]=fltokens(fline); dresp.best.eval= tokens{1}{8}; fline=fgetl(fidi); - while ischar(fline) && ~strncmpi(fline,'<<<<< Best ',11) - fline=fgetl(fidi); - end - -% read until next best or end + while ischar(fline) && ~isempty(fline) && ... + ~strncmpi(fline,'<<<<< Best ',11) + fline=fgetl(fidi); + end + +% read until next best or blank or end else - display([' ' deblank(fline) ' (not processed)']); - - fline=fgetl(fidi); - - while ischar(fline) && ~strncmpi(fline,'<<<<< Best ',11) + display([' ' deblank(fline) ' (ignored)']); + + fline=fgetl(fidi); + + while ischar(fline) && ~isempty(fline) && ... + ~strncmpi(fline,'<<<<< Best ',11) fline=fgetl(fidi); end @@ -775,4 +791,27 @@ end +%% function to find and read the iterator completion + +function [method]=itcomp_read(fidi,fline) + +if ~exist('fline','var') || isempty(fline) || ~ischar(fline) + while 1 + [fline]=findline(fidi,'<<<<< Iterator '); + if ~ischar(fline) + return + end + if (length(fline) > 26) && ... + ~isempty(strfind(fline(16:end),' completed.')) + break + end + end +end + +[ntokens,tokens]=fltokens(fline); +method=tokens{1}{3}; +display(sprintf('Dakota iterator ''%s'' completed.',method)); + +end + %% function to find a file line starting with a specified string Index: /issm/trunk/src/m/solutions/dakota/plot_boxplot.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/plot_boxplot.m (revision 3091) +++ /issm/trunk/src/m/solutions/dakota/plot_boxplot.m (revision 3092) @@ -2,8 +2,16 @@ % plot a box plot of the responses. % -% []=plot_boxplot(dresp) +% []=plot_boxplot(dresp ,params) +% []=plot_boxplot(dresp,descr,params) +% []=plot_boxplot(sampr,descr,params) % % where the required input is: % dresp (structure array, responses) +% or +% dresp (structure array, responses) +% descr (cell array, list of response descriptions desired) +% or +% sampr (double array, lists of response samples) +% descr (cell array, list of response descriptions) % % the required fields of dresp are: @@ -11,5 +19,12 @@ % sample (double vector, list of samples) % -% for each response in the input array, this functionplots a +% the optional input is: +% params (string/numeric, parameter names and values) +% +% where the optional parameters are: +% ymin (numeric, minimum of y-axis) +% ymax (numeric, maximum of y-axis) +% +% for each response in the input array, this function plots a % matlab box plot of the list of samples and annotates it % with the description. the lists of samples need not all be @@ -23,5 +38,5 @@ % acknowledged. Any commercial use must be negotiated with % the Office of Technology Transfer at the California Institute -% of Technology. (NTR 47078) +% of Technology. (J. Schiermeier, NTR 47078) % % This software may be subject to U.S. export control laws. @@ -32,5 +47,5 @@ % to foreign countries or providing access to foreign persons." % -function []=plot_boxplot(dresp) +function []=plot_boxplot(varargin) if ~nargin @@ -39,29 +54,92 @@ end -%% assemble the data into a matrix +%% process input data and assemble into matrices as needed -desc=cell (1,length(dresp)); -for i=1:length(dresp) - ldata(i)=length(dresp(i).sample); -end -data=zeros(max(ldata),length(dresp)); +% responses -for i=1:length(dresp) - desc(i)=cellstr(dresp(i).descriptor); - data(1:ldata(i),i)=dresp(i).sample; +iarg=1; +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp=struc_desc(dresp,varargin{iarg}); + iarg=iarg+1; + end + + descr=cell (1,length(dresp)); + lsamp=zeros(1,length(dresp)); + for i=1:length(dresp) + lsamp(i)=length(dresp(i).sample); + end + sampr=zeros(max(lsamp),length(dresp)); + sampr(:,:)=NaN; + + for i=1:length(dresp) + descr(i)=cellstr(dresp(i).descriptor); + sampr(1:lsamp(i),i)=dresp(i).sample; + end +else + sampr=varargin{iarg}; + iarg=iarg+1; + + if iarg <= nargin && iscell(varargin{iarg}) + descr=varargin{iarg}; + iarg=iarg+1; +% elseif iarg <= nargin && ischar(varargin{iarg}) +% descr=cellstr(varargin{iarg}); +% iarg=iarg+1; + else + descr=cell(1,size(sampr,2)); + end end +for i=1:length(descr) + if isempty(descr{i}) + descr(i)={['resp_' num2str(i)]}; + end +end + +% parameters + +while (iarg <= nargin-1) + if ischar(varargin{iarg}) + eval([varargin{iarg} '=varargin{iarg+1};']); + disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']); + else + error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']); + end + iarg=iarg+2; +end + %% draw the plot % draw box plot -boxplot(data,'labels',desc,'notch','on') +figure +boxplot(sampr,'labels',descr,'notch','on') ax1=gca; % add the annotation -title('Box Plot of Variables and/or Responses') -xlabel('Variable or Response') -ylabel('Value') +ylim('auto') +[ylims]=ylim; +if exist('ymin','var') + ylims(1)=ymin; +end +if exist('ymax','var') + ylims(2)=ymax; +end +ylim(ylims) + +if (size(sampr,2) == 1) + tlabc=descr{1}; +else + tlabc='Responses'; +end +title(['Box Plot of ' tlabc],'Interpreter','none'); +xlabel('Response','Interpreter','none'); +ylabel('Value' ,'Interpreter','none'); end Index: /issm/trunk/src/m/solutions/dakota/plot_hist_norm.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/plot_hist_norm.m (revision 3092) +++ /issm/trunk/src/m/solutions/dakota/plot_hist_norm.m (revision 3092) @@ -0,0 +1,367 @@ +% +% plot a relative histogram and cdf optionally along with +% a normal distribution. +% +% []=plot_hist_norm(dresp1 ,dresp2 ,params) +% []=plot_hist_norm(dresp1,desc1,dresp2,desc2,params) +% []=plot_hist_norm(sampr ,descr,mu,sigma ,params) +% +% where the required input is: +% dresp1 (structure array, responses) +% or +% dresp1 (structure array, responses) +% desc1 (cell array, list of response descriptions desired) +% or +% sampr (double array, lists of samples) +% descr (cell array, list of descriptions) +% +% and the optional input is: +% dresp2 (structure array, responses) +% or +% dresp2 (structure array, responses) +% desc2 (cell array, list of response descriptions desired) +% or +% mu (double vector, means) +% sigma (double vector, standard deviations) +% +% the required fields of dresp1 are: +% descriptor (char, description) +% sample (double vector, list of samples) +% +% and the required fields of dresp2 are: +% mean (double, mean of sample) +% stddev (double, standard deviation of sample) +% +% the optional input is: +% params (string/numeric, parameter names and values) +% +% where the optional parameters are: +% hmin (numeric, minimum for histogram) +% hmax (numeric, maximum for histogram) +% hnint (numeric, number of intervals for histogram) +% ymin1 (numeric, minimum of histogram y-axis) +% ymax1 (numeric, maximum of histogram y-axis) +% ymin2 (numeric, minimum of cdf y-axis) +% ymax2 (numeric, maximum of cdf y-axis) +% cdfplt (char, 'off' to turn off cdf line plots) +% cdfleg (char, 'off' to turn off cdf legends) +% +% for each response in the input array, this function +% calculates and plots a relative histogram and CDF of the list +% of samples, and annotates it with the description. in +% addition, a mean and standard deviation may be supplied or, +% if empty, calculated so that a normal distribution and CDF may +% be plotted. +% +% dresp1 data would typically be contained in the dakota tabular +% output file from a sampling analysis, and dresp2 data would +% typically be contained in the dakota output file from a local +% sensitivity analysis, both read by dakota_out_parse. +% +% "Copyright 2009, by the California Institute of Technology. +% ALL RIGHTS RESERVED. United States Government Sponsorship +% acknowledged. Any commercial use must be negotiated with +% the Office of Technology Transfer at the California Institute +% of Technology. (J. Schiermeier, NTR 47078) +% +% This software may be subject to U.S. export control laws. +% By accepting this software, the user agrees to comply with +% all applicable U.S. export laws and regulations. User has the +% responsibility to obtain export licenses, or other export +% authority as may be required before exporting such information +% to foreign countries or providing access to foreign persons." +% +function []=plot_hist_norm(varargin) + +if ~nargin + help plot_hist_norm + return +end + +%% process input data and assemble into matrices as needed + +% responses for histograms + +iarg=1; +if isstruct(varargin{iarg}) + dresp1=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp1=struc_desc(dresp1,varargin{iarg}); + iarg=iarg+1; + end + + descr=cell (1,length(dresp1)); + lsamp=zeros(1,length(dresp1)); + for i=1:length(dresp1) + lsamp(i)=length(dresp1(i).sample); + end + sampr=zeros(max(lsamp),length(dresp1)); + sampr(:,:)=NaN; + + for i=1:length(dresp1) + descr(i)=cellstr(dresp1(i).descriptor); + sampr(1:lsamp(i),i)=dresp1(i).sample; + end +else + sampr=varargin{iarg}; + iarg=iarg+1; + + lsamp(1:size(sampr,2))=size(sampr,1); + + if iarg <= nargin && iscell(varargin{iarg}) + descr=varargin{iarg}; + iarg=iarg+1; +% elseif iarg <= nargin && ischar(varargin{iarg}) +% descr=cellstr(varargin{iarg}); +% iarg=iarg+1; + else + descr=cell(1:size(sampr,2)); + end +end + +for i=1:length(descr) + if isempty(descr{i}) + descr(i)={['resp_' num2str(i)]}; + end +end + +% responses for normal distributions + +if iarg <= nargin && isstruct(varargin{iarg}) + dresp2=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp2=struc_desc(dresp2,varargin{iarg}); + iarg=iarg+1; + end + + mu =zeros(1,length(dresp2)); + sigma=zeros(1,length(dresp2)); + + for i=1:length(dresp2) + mu (i)=dresp2(i).mean; + sigma(i)=dresp2(i).stddev; + end +elseif iarg+1 <= nargin && ... + isnumeric(varargin{iarg}) && isnumeric(varargin{iarg+1}) + if ~isempty(varargin{iarg}) + mu =varargin{iarg}; + else + mu =mean(sampr); + display('Using calculated means.') + end + iarg=iarg+1; + if ~isempty(varargin{iarg}) + sigma=varargin{iarg}; + else + sigma=std(sampr); + display('Using calculated standard deviations.') + end + iarg=iarg+1; +end + +% parameters + +while (iarg <= nargin-1) + if ischar(varargin{iarg}) + eval([varargin{iarg} '=varargin{iarg+1};']); + disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']); + else + error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']); + end + iarg=iarg+2; +end + +%% generate the intervals + +if ~exist('hmin','var') + hmin=min(min(sampr)); +end +if ~exist('hmax','var') + hmax=max(max(sampr)); +end +if ~exist('hnint','var') + hnint=50; +end +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(sampr,edges); +for i=1:size(sampr,2) + dbelow(i) =length(find(sampr(:,i)edges(end)))/lsamp(i); +end + +if exist('mu','var') && exist('sigma','var') + ncol=size(sampr,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('descr','var') + descr(ncol+i)={[descr{i} ' norm']}; + end + end +end + +% draw the bar plot + +figure +hl1=bar(edges(1:end-1),dhistc(1:end-1,:)); +ax1=gca; + +% set barseries properties for clarity + +if (length(hl1) > 1) + for i=1:length(hl1) + set(hl1(i),'BarWidth',1,'EdgeColor','none'); + end +end + +xlim('auto') +[xlims]=xlim; +if exist('hmin','var') + xlims(1)=edges(1); +end +if exist('hmax','var') + xlims(2)=edges(end-1); +end +xlim(xlims) + +ylim('auto') +[ylims]=ylim; +if exist('ymin1','var') + ylims(1)=ymin1; +end +if exist('ymax1','var') + ylims(2)=ymax1; +end +ylim(ylims) + +% add the annotation + +if exist('cdfplt','var') && strcmpi(cdfplt,'off') + title('Relative Frequency Histogram') +else + title('Relative Frequency Histogram with CDF') +end +xlabel('Interval Edge Value'); +ylabel('Relative Frequency'); + +if exist('descr','var') + hleg1=legend(ax1,descr,'Location','NorthWest',... + 'Color','none','Interpreter','none'); +else + hleg1=legend(ax1); +end + +%% generate the cumulative distribution functions + +if ~exist('cdfplt','var') || ~strcmpi(cdfplt,'off') +% 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('descr','var') + ncol=length(descr); + for i=1:ncol + cdescr(i)={[descr{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 on the right side, +% so make equal number of ticks for the line plot on right side + + nytick=length(get(ax1,'YTick')); +% ylim('auto') +% [ylims]=ylim; + [ylims]=[0 ceil(max(max(cdf))/0.1-0.1)*0.1]; + if exist('ymin2','var') + ylims(1)=ymin2; + end + if exist('ymax2','var') + ylims(2)=ymax2; + else + ylims(2)=ylims(1)+(nytick-1)/(nytick-1-1)*(ylims(2)-ylims(1)); + end +% ylim(ylims) + ytinc =(ylims(2)-ylims(1))/(nytick-1); + + ax2=axes('Position',get(ax1,'Position'),... + 'XLim',get(ax1,'XLim'),... + 'XTick',get(ax1,'XTick'),... + 'YLim',ylims,... + 'YTick',[ylims(1):ytinc:ylims(2)],... + 'XAxisLocation','bottom','YAxisLocation','right',... + 'Color','none','Layer','top'); + hl2=line(edges(1:end-1),cdf(1:end-1,:),'Parent',ax2); + +% set line property colors to match barseries property +% (if barseries is "flat", must interpolate and round to colormap) + + cmap=colormap; + for i=1:length(hl2) + if ischar(get(hl1(i),'FaceColor')) && ... + strcmpi(get(hl1(i),'FaceColor'),'flat') + if (length(hl2) > 1) + imap=round((i-1)/(length(hl2)-1)*(size(cmap,1)-1))+1; + else + imap=1; + end + set(hl2(i),'Color',cmap(imap,:)) + else + set(hl2(i),'Color',get(hl1(i),'FaceColor')) + end + end + +% add the annotation + + ylabel('Cumulative Percent'); + + if ~exist('cdfleg','var') || ~strcmpi(cdfleg,'off') +% legend doesn't combine with bar chart above + if exist('cdescr','var') + hleg2=legend(ax2,cdescr,'Location','NorthEast',... + 'Color','none','Interpreter','none'); +% set(hleg2,'Color','white') + else + hleg2=legend(ax2); +% set(hleg2,'Color','white') + end + end + + set(gcf,'PaperPositionMode','auto') +% hold off +end + +end Index: /issm/trunk/src/m/solutions/dakota/plot_hist_norm_ci.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/plot_hist_norm_ci.m (revision 3092) +++ /issm/trunk/src/m/solutions/dakota/plot_hist_norm_ci.m (revision 3092) @@ -0,0 +1,455 @@ +% +% plot a relative histogram and cdf optionally along with a +% normal distribution for the sample and confidence intervals. +% +% []=plot_hist_norm_ci(dresp ,params) +% []=plot_hist_norm_ci(dresp,descr,params) +% []=plot_hist_norm_ci(sampr,descr,params) +% +% where the required input is: +% dresp (structure array, responses) +% or +% dresp (structure array, responses) +% descr (cell array, list of response descriptions desired) +% or +% sampr (double array, lists of samples) +% descr (cell array, list of descriptions) +% +% the required fields of dresp are: +% descriptor (char, description) +% sample (double vector, list of samples) +% +% and the optional fields of dresp are: +% mean (double, mean of sample) +% stddev (double, standard deviation of sample) +% meanci(2) (double, confidence interval of mean) +% stddevci(2) (double, confidence interval of standard deviation) +% +% the optional input is: +% params (string/numeric, parameter names and values) +% +% and the optional input is: +% hmin (numeric, minimum for histogram) +% hmax (numeric, maximum for histogram) +% hnint (numeric, number of intervals for histogram) +% ymin1 (numeric, minimum of histogram y-axis) +% ymax1 (numeric, maximum of histogram y-axis) +% ymin2 (numeric, minimum of cdf y-axis) +% ymax2 (numeric, maximum of cdf y-axis) +% nrmplt (char, 'line' or 'off' to change nrm plots from 'bar') +% ciplt (char, 'line' or 'off' to change ci plots from 'bar') +% cdfplt (char, 'off' to turn off cdf line plots) +% cdfleg (char, 'off' to turn off cdf legends) +% cmap (char or numeric, colormap definition) +% +% for each response in the input array, this function +% calculates and plots a relative histogram and CDF of the list +% of samples, and annotates it with the description. in +% addition, the normal distribution and CDF are plotted, and +% four CDF's are plotted for the confidence intervals. +% +% dresp data would typically be contained in the dakota tabular +% output file from a sampling analysis, read by dakota_out_parse. +% +% "Copyright 2009, by the California Institute of Technology. +% ALL RIGHTS RESERVED. United States Government Sponsorship +% acknowledged. Any commercial use must be negotiated with +% the Office of Technology Transfer at the California Institute +% of Technology. (J. Schiermeier, NTR 47078) +% +% This software may be subject to U.S. export control laws. +% By accepting this software, the user agrees to comply with +% all applicable U.S. export laws and regulations. User has the +% responsibility to obtain export licenses, or other export +% authority as may be required before exporting such information +% to foreign countries or providing access to foreign persons." +% +function []=plot_hist_norm_ci(varargin) + +if ~nargin + help plot_hist_norm_ci + return +end + +%% process input data and assemble into matrices as needed + +iarg=1; +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp=struc_desc(dresp,varargin{iarg}); + iarg=iarg+1; + end + + descr=cell (1,length(dresp)); + lsamp=zeros(1,length(dresp)); + for i=1:length(dresp) + lsamp(i)=length(dresp(i).sample); + end + sampr=zeros(max(lsamp),length(dresp)); + sampr(:,:)=NaN; + + mu =zeros(1,length(dresp)); + sigma =zeros(1,length(dresp)); + muci =zeros(2,length(dresp)); + sigmaci=zeros(2,length(dresp)); + + for i=1:length(dresp) + descr(i)=cellstr(dresp(i).descriptor); + sampr(1:lsamp(i),i)=dresp(i).sample; + mu (i) =dresp(i).mean; + sigma (i) =dresp(i).stddev; + muci (:,i)=dresp(i).meanci; + sigmaci(:,i)=dresp(i).stddevci; + end +else + sampr=varargin{iarg}; + iarg=iarg+1; + + lsamp(1:size(sampr,2))=size(sampr,1); + + if iarg <= nargin && iscell(varargin{iarg}) + descr=varargin{iarg}; + iarg=iarg+1; +% elseif iarg <= nargin && ischar(varargin{iarg}) +% descr=cellstr(varargin{iarg}); +% iarg=iarg+1; + else + descr=cell(1,size(sampr,2)); + end + + mu =zeros(1,size(sampr,2)); + sigma =zeros(1,size(sampr,2)); + muci =zeros(2,size(sampr,2)); + sigmaci=zeros(2,size(sampr,2)); + for i=1:size(sampr,2) +% calculate 95% confidence intervals (same as dakota) + [mu(i),sigma(i),muci(:,i),sigmaci(:,i)]=... + normfit(sampr(:,i),0.05); + end + display('Using calculated normal fits.') +end + +for i=1:length(descr) + if isempty(descr{i}) + descr(i)={['resp_' num2str(i)]}; + end +end + +% parameters + +while (iarg <= nargin-1) + if ischar(varargin{iarg}) + eval([varargin{iarg} '=varargin{iarg+1};']); + disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']); + else + error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']); + end + iarg=iarg+2; +end + +%% generate the intervals + +if ~exist('hmin','var') + hmin=min(min(sampr)); +end +if ~exist('hmax','var') + hmax=max(max(sampr)); +end +if ~exist('hnint','var') + hnint=50; +end +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(sampr,edges); +for i=1:size(sampr,2) + dbelow(i) =length(find(sampr(:,i)edges(end)))/lsamp(i); +end + +if exist('mu','var') && exist('sigma','var') + ncol=size(sampr,2); + for i=1:ncol + dbelow(end+1)=normcdf(edges( 1),mu(i),sigma(i)); + dhistc(1:size(dhistc,1)-1,end+1)=... + normcdf(edges(2:end ),mu(i),sigma(i))-... + normcdf(edges(1:end-1),mu(i),sigma(i)); + dabove(end+1)=norminv(edges(end),mu(i),sigma(i)); + if exist('descr','var') + descr(end+1)={[descr{i} ' norm']}; + end + end +end + +if exist('muci','var') && exist('sigmaci','var') + for i=1:ncol + dbelow(end+1)=normcdf(edges( 1),muci(1,i),sigmaci(2,i)); + dhistc(1:size(dhistc,1)-1,end+1)=... + normcdf(edges(2:end ),muci(1,i),sigmaci(2,i))-... + normcdf(edges(1:end-1),muci(1,i),sigmaci(2,i)); + dabove(end+1)=norminv(edges(end),muci(1,i),sigmaci(2,i)); + if exist('descr','var') + descr(end+1)={[descr{i} ' norm-+']}; + end + end + for i=1:ncol + dbelow(end+1)=normcdf(edges( 1),muci(1,i),sigmaci(1,i)); + dhistc(1:size(dhistc,1)-1,end+1)=... + normcdf(edges(2:end ),muci(1,i),sigmaci(1,i))-... + normcdf(edges(1:end-1),muci(1,i),sigmaci(1,i)); + dabove(end+1)=norminv(edges(end),muci(1,i),sigmaci(1,i)); + if exist('descr','var') + descr(end+1)={[descr{i} ' norm--']}; + end + end + for i=1:ncol + dbelow(end+1)=normcdf(edges( 1),muci(2,i),sigmaci(1,i)); + dhistc(1:size(dhistc,1)-1,end+1)=... + normcdf(edges(2:end ),muci(2,i),sigmaci(1,i))-... + normcdf(edges(1:end-1),muci(2,i),sigmaci(1,i)); + dabove(end+1)=norminv(edges(end),muci(2,i),sigmaci(1,i)); + if exist('descr','var') + descr(end+1)={[descr{i} ' norm+-']}; + end + end + for i=1:ncol + dbelow(end+1)=normcdf(edges( 1),muci(2,i),sigmaci(2,i)); + dhistc(1:size(dhistc,1)-1,end+1)=... + normcdf(edges(2:end ),muci(2,i),sigmaci(2,i))-... + normcdf(edges(1:end-1),muci(2,i),sigmaci(2,i)); + dabove(end+1)=norminv(edges(end),muci(2,i),sigmaci(2,i)); + if exist('descr','var') + descr(end+1)={[descr{i} ' norm++']}; + end + end +end + +% draw the bar plot + +figure +if ~exist('nrmplt','var') + nrmplt='bar'; +end +if ~exist('ciplt','var') + ciplt='bar'; +end + +hold all +% hl1=bar(edges(1:end-1),dhistc(1:end-1,1:6*ncol)); +% hl1=line(edges(1:end-1),dhistc(1:end-1,1:6*ncol)); +if strcmpi(nrmplt,'bar') + if strcmpi(ciplt,'bar') + hl1=bar (edges(1:end-1),dhistc(1:end-1,1:6*ncol)); + elseif strcmpi(ciplt,'line') + hl1=bar (edges(1:end-1),dhistc(1:end-1,1:2*ncol)); + hl1(2*ncol+1:6*ncol)=... + line(edges(1:end-1),dhistc(1:end-1,2*ncol+1:6*ncol),... + 'LineWidth',2); + elseif strcmpi(ciplt,'off') + hl1=bar (edges(1:end-1),dhistc(1:end-1,1:2*ncol)); + end +elseif strcmpi(nrmplt,'line') + if strcmpi(ciplt,'bar') || strcmpi(ciplt,'line') + hl1=bar (edges(1:end-1),dhistc(1:end-1,1:1*ncol)); + hl1(1*ncol+1:2*ncol)=... + line(edges(1:end-1),dhistc(1:end-1,1*ncol+1:2*ncol),... + 'LineWidth',2); + hl1(2*ncol+1:6*ncol)=... + line(edges(1:end-1),dhistc(1:end-1,2*ncol+1:6*ncol),... + 'LineWidth',1); + elseif strcmpi(ciplt,'off') + hl1=bar (edges(1:end-1),dhistc(1:end-1,1:1*ncol)); + hl1(1*ncol+1:2*ncol)=... + line(edges(1:end-1),dhistc(1:end-1,1*ncol+1:2*ncol),... + 'LineWidth',2); + end +elseif strcmpi(nrmplt,'off') + hl1=bar (edges(1:end-1),dhistc(1:end-1,1:1*ncol)); +end +ax1=gca; +hold off + +% set barseries properties for clarity + +if (length(hl1) > 1) + for i=1:length(hl1) + if strcmpi(get(hl1(i),'Type'),'hggroup') + set(hl1(i),'BarWidth',1,'EdgeColor','none'); + end + end +end + +% set bars and lines to have a continuous colormap +% (if barseries is "flat", must interpolate and round to colormap) + +if exist('cmap','var') + colormap(cmap) +end + +cmap=colormap; +for i=1:length(hl1) + if (length(hl1) > 1) + imap=round((i-1)/(length(hl1)-1)*(size(cmap,1)-1))+1; + else + imap=1; + end + if strcmpi(get(hl1(i),'Type'),'hggroup') + if ischar(get(hl1(i),'FaceColor')) && ... + strcmpi(get(hl1(i),'FaceColor'),'flat') + set(hl1(i),'FaceColor',cmap(imap,:)) + else + set(hl1(i),'FaceColor',get(hl1(i),'FaceColor')) + end + elseif strcmpi(get(hl1(i),'Type'),'line') + set(hl1(i),'Color',cmap(imap,:)) + end +end + +xlim('auto') +[xlims]=xlim; +if exist('hmin','var') + xlims(1)=edges(1); +end +if exist('hmax','var') + xlims(2)=edges(end-1); +end +xlim(xlims) + +ylim('auto') +[ylims]=ylim; +if exist('ymin1','var') + ylims(1)=ymin1; +end +if exist('ymax1','var') + ylims(2)=ymax1; +end +ylim(ylims) + +% add the annotation + +if exist('cdfplt','var') && strcmpi(cdfplt,'off') + title('Relative Frequency Histogram') +else + title('Relative Frequency Histogram with CDF') +end +xlabel('Interval Edge Value'); +ylabel('Relative Frequency'); + +if exist('descr','var') + hleg1=legend(ax1,descr(1:length(hl1)),'Location','NorthWest',... + 'Color','none','Interpreter','none'); +else + hleg1=legend(ax1); +end + +%% generate the cumulative distribution functions + +if ~exist('cdfplt','var') || ~strcmpi(cdfplt,'off') +% 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('descr','var') + ncol=length(descr); + for i=1:ncol + cdescr(i)={[descr{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 on the right side, +% so make equal number of ticks for the line plot on right side + + nytick=length(get(ax1,'YTick')); +% ylim('auto') +% [ylims]=ylim; + [ylims]=[0 ceil(max(max(cdf))/0.1-0.1)*0.1]; + if exist('ymin2','var') + ylims(1)=ymin2; + end + if exist('ymax2','var') + ylims(2)=ymax2; + else + ylims(2)=ylims(1)+(nytick-1)/(nytick-1-1)*(ylims(2)-ylims(1)); + end +% ylim(ylims) + ytinc =(ylims(2)-ylims(1))/(nytick-1); + + ax2=axes('Position',get(ax1,'Position'),... + 'XLim',get(ax1,'XLim'),... + 'XTick',get(ax1,'XTick'),... + 'YLim',ylims,... + 'YTick',[ylims(1):ytinc:ylims(2)],... + 'XAxisLocation','bottom','YAxisLocation','right',... + 'Color','none','Layer','top'); + hl2=line(edges(1:end-1),cdf(1:end-1,1:length(hl1)),'Parent',ax2); + +% set line property colors to match barseries or line property +% (if barseries is "flat", must interpolate and round to colormap) + + cmap=colormap; + for i=1:length(hl2) + if (length(hl2) > 1) + imap=round((i-1)/(length(hl2)-1)*(size(cmap,1)-1))+1; + else + imap=1; + end + if strcmpi(get(hl1(i),'Type'),'hggroup') + if ischar(get(hl1(i),'FaceColor')) && ... + strcmpi(get(hl1(i),'FaceColor'),'flat') + set(hl2(i),'Color',cmap(imap,:)) + else + set(hl2(i),'Color',get(hl1(i),'FaceColor')) + end + elseif strcmpi(get(hl1(i),'Type'),'line') + set(hl2(i),'Color',get(hl1(i),'Color')) + end + end + +% add the annotation + + ylabel('Cumulative Percent'); + + if ~exist('cdfleg','var') || ~strcmpi(cdfleg,'off') +% legend doesn't combine with bar chart above + if exist('cdescr','var') + hleg2=legend(ax2,cdescr(1:length(hl2)),'Location','NorthEast',... + 'Color','none','Interpreter','none'); +% set(hleg2,'Color','white') + else + hleg2=legend(ax2); +% set(hleg2,'Color','white') + end + end + + set(gcf,'PaperPositionMode','auto') +% hold off +end + +end Index: /issm/trunk/src/m/solutions/dakota/plot_if_bars.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/plot_if_bars.m (revision 3091) +++ /issm/trunk/src/m/solutions/dakota/plot_if_bars.m (revision 3092) @@ -2,19 +2,29 @@ % plot a stacked bar chart of the importance factors. % -% []=plot_if_bars(dresp,ifmin,isort) +% []=plot_if_bars(dresp ,params) +% []=plot_if_bars(dresp,descr,params) % % where the required input is: % dresp (structure array, responses) -% -% and the optional input is: +% or +% dresp (structure array, responses) +% descr (cell array, list of response descriptions desired) +% +% the required fields of dresp are: +% descriptor (char, description) +% var (cell array, variables) +% impfac (double array, importance factors) +% +% the optional input is: +% params (string/numeric, parameter names and values) +% +% where the optional parameters are: +% ymin (numeric, minimum of y-axis) +% ymax (numeric, maximum of y-axis) % ifmin (double, minimum importance factor) % isort (numeric, sort flag: 0, no sorting; % 1, highest at bottom; % -1, lowest at bottom) -% -% the required fields of dresp are: -% descriptor (char, description) -% var (cell array, variables) -% impfac (double array, importance factors) +% xtlrot (numeric, rotation in degrees of x-tick labels) % % for each response in the input array, this function plots @@ -32,5 +42,5 @@ % acknowledged. Any commercial use must be negotiated with % the Office of Technology Transfer at the California Institute -% of Technology. (NTR 47078) +% of Technology. (J. Schiermeier, NTR 47078) % % This software may be subject to U.S. export control laws. @@ -41,5 +51,5 @@ % to foreign countries or providing access to foreign persons." % -function []=plot_if_bars(dresp,ifmin,isort) +function []=plot_if_bars(varargin) if ~nargin @@ -48,4 +58,47 @@ end +%% process input data and assemble into matrices + +% responses + +iarg=1; +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp=struc_desc(dresp,varargin{iarg}); + iarg=iarg+1; + end + + descr=cell (1,length(dresp)); + lifr =zeros(1,length(dresp)); + for i=1:length(dresp) + lifr(i)=length(dresp(i).impfac); + end + ifr =zeros(length(dresp),max(lifr)); + dvar=dresp(find(lifr == max(lifr),1,'first')).var; + + for i=1:length(dresp) + descr(i)=cellstr(dresp(i).descriptor); + ifr(i,1:lifr(i))=dresp(i).impfac; + end +else + error(['''' inputname(iarg) ''' is not a structure.']); +end + +% parameters + +while (iarg <= nargin-1) + if ischar(varargin{iarg}) + eval([varargin{iarg} '=varargin{iarg+1};']); + disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']); + else + error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']); + end + iarg=iarg+2; +end + if ~exist('ifmin','var') || isempty(ifmin) ifmin=0; @@ -56,22 +109,8 @@ end -%% assemble the data into a matrix - -desc=cell (1,length(dresp)); -for i=1:length(dresp) - ldata(i)=length(dresp(i).impfac); -end -data=zeros(length(dresp),max(ldata)); -dvar=dresp(find(ldata == max(ldata),1,'first')).var; - -for i=1:length(dresp) - desc(i)=cellstr(dresp(i).descriptor); - data(i,1:ldata(i))=dresp(i).impfac; -end - %% sort the data, if necessary if (isort) - ifmean=mean(data,1); + ifmean=mean(ifr,1); if (isort > 0) [ifmean,index]=sort(ifmean,'descend'); @@ -82,5 +121,5 @@ dvar=dvar(index); - data=data(:,index); + ifr =ifr (:,index); end @@ -90,19 +129,19 @@ nif=length(dvar); dvar(nif+1,1)=cellstr(sprintf('others < %f',ifmin)); - data(:,nif+1)=0.; + ifr (:,nif+1)=0.; nif2=0; dvar2=cell (size(dvar)); - data2=zeros(size(data)); + ifr2 =zeros(size(ifr )); % sum filtered rows and copy unfiltered rows for i=1:nif - if (max(data(:,i)) < ifmin) - data(:,nif+1)=data(:,nif+1)+data(:,i); + if (max(ifr(:,i)) < ifmin) + ifr(:,nif+1)=ifr(:,nif+1)+ifr(:,i); else nif2=nif2+1; dvar2(nif2) =dvar(i); - data2(:,nif2)=data(:,i); + ifr2 (:,nif2)=ifr (:,i); end end @@ -111,19 +150,19 @@ dvar2(nif2+1) =dvar(nif+1); - data2(:,nif2+1)=data(:,nif+1); + ifr2 (:,nif2+1)=ifr (:,nif+1); % copy back and truncate filtered rows - clear dvar data + clear dvar ifr if (isort >= 0) dvar(1:nif2+1) =dvar2(1:nif2+1); - data(:,1:nif2+1)=data2(:,1:nif2+1); + ifr (:,1:nif2+1)=ifr2 (:,1:nif2+1); else dvar(1 ) =dvar2( nif2+1); dvar(2:nif2+1) =dvar2(1:nif2 ); - data(:,1 )=data2(:, nif2+1); - data(:,2:nif2+1)=data2(:,1:nif2 ); - end - clear nif nif2 dvar2 data2 + ifr (:,1 )=ifr2 (:, nif2+1); + ifr (:,2:nif2+1)=ifr2 (:,1:nif2 ); + end + clear nif nif2 dvar2 ifr2 end @@ -134,15 +173,36 @@ if length(dresp) == 1 - data=[data; data]; -end - -hl1=bar(data,'stacked'); + ifr=[ifr; ifr]; +end + +figure +hl1=bar(ifr,'stacked'); ax1=gca; -set(ax1,'ylim',[0 1.2]) if length(dresp) == 1 set(ax1,'xlim',[0.5 1.5]) end -set(ax1,'xticklabel',desc) + +% set(ax1,'ylim',[0 1.2]) +% ylim('auto') +% [ylims]=ylim; +[ylims]=[0 1.2]; +if exist('ymin','var') + ylims(1)=ymin; +end +if exist('ymax','var') + ylims(2)=ymax; +end +ylim(ylims) + +set(ax1,'xtick',1:length(descr)) +set(ax1,'xticklabel',descr) +if exist('xtlrot','var') + htl=rotateticklabel(ax1,xtlrot); + tlext=zeros(length(htl),4); + for i=1:length(htl) + tlext(i,:)=get(htl(i),'Extent'); + end +end % add the annotation @@ -150,4 +210,10 @@ title('Importance Factors') xlabel('Response') +if exist('xtlrot','var') + xlext=get(get(ax1,'xlabel'),'Extent'); + nskip=ceil(max(tlext(:,4))/xlext(4)); + xlabel(cellstr([repmat(' ',nskip,1);'Response'])); + clear nskip xlext tlext +end ylabel('Importance Factor Value') Index: /issm/trunk/src/m/solutions/dakota/plot_normdist_bars.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/plot_normdist_bars.m (revision 3091) +++ /issm/trunk/src/m/solutions/dakota/plot_normdist_bars.m (revision 3092) @@ -2,12 +2,16 @@ % plot a stacked bar chart of the sample distributions. % -% []=plot_normdist_bars(dresp,plev,lstr) +% []=plot_normdist_bars(dresp ,params) +% []=plot_normdist_bars(dresp,descr,params) +% []=plot_normdist_bars(sampr,descr,params) % % where the required input is: % dresp (structure array, responses) -% -% and the optional input is: -% plev (double vector, probability levels) -% lstr (cell array, legend labels) +% or +% dresp (structure array, responses) +% descr (cell array, list of response descriptions desired) +% or +% sampr (double array, lists of response samples) +% descr (cell array, list of response descriptions) % % the required fields of dresp are: @@ -18,4 +22,14 @@ % mean (double, mean of sample) % stddev (double, standard deviation of sample) +% +% the optional input is: +% params (string/numeric, parameter names and values) +% +% where the optional parameters are: +% prob (double vector, probability levels) +% ymin (numeric, minimum of y-axis) +% ymax (numeric, maximum of y-axis) +% xtlrot (numeric, rotation in degrees of x-tick labels) +% lstr (cell array, legend labels) % % for each response in the input array, this function plots @@ -35,5 +49,5 @@ % acknowledged. Any commercial use must be negotiated with % the Office of Technology Transfer at the California Institute -% of Technology. (NTR 47078) +% of Technology. (J. Schiermeier, NTR 47078) % % This software may be subject to U.S. export control laws. @@ -44,5 +58,5 @@ % to foreign countries or providing access to foreign persons." % -function []=plot_normdist_bars(dresp,plev,lstr) +function []=plot_normdist_bars(varargin) if ~nargin @@ -51,28 +65,78 @@ end -if ~exist('plev','var') || isempty(plev) - plev=[0.01 0.25 0.50 0.75 0.99]; -end - -%% assemble the data into a matrix and calculate the increments - -desc=cell (1,length(dresp)); -data=zeros(length(dresp),5); +%% process input data and assemble into dresp as needed + +% responses + +iarg=1; +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp=struc_desc(dresp,varargin{iarg}); + iarg=iarg+1; + end +else + sampr=varargin{iarg}; + iarg=iarg+1; + + if iarg <= nargin && iscell(varargin{iarg}) + descr=varargin{iarg}; + iarg=iarg+1; +% elseif iarg <= nargin && ischar(varargin{iarg}) +% descr=cellstr(varargin{iarg}); +% iarg=iarg+1; + else + descr=cell(1:size(sampr,2)); + end + + dresp=struct([]); + for i=1:size(sampr,2) + dresp(end+1).sample=samp(:,i); + if ~isempty(descr) + dresp(i).descriptor=descr{i}; + else + dresp(i).descriptor=['dresp_' num2str(i)]; + end + end +end + +% parameters + +while (iarg <= nargin-1) + if ischar(varargin{iarg}) + eval([varargin{iarg} '=varargin{iarg+1};']); + disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']); + else + error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']); + end + iarg=iarg+2; +end + +%% calculate any missing information (noting that dresp is local) if ~isfield(dresp,'mean') || ~isfield(dresp,'stddev') for i=1:length(dresp) -% calculate 95% confidence intervals (same as dakota) - [dresp(i).mean,dresp(i).stddev,... - dresp(i).meanci,dresp(i).stddevci]=... - normfit(dresp(i).sample,0.05); - end -end + [dresp(i).mean,dresp(i).stddev]=normfit(dresp(i).sample); + end +end + +%% assemble the data into a matrix and calculate the increments + +if ~exist('prob','var') || isempty(prob) + prob=[0.01 0.25 0.50 0.75 0.99]; +end + +descr=cell (1,length(dresp)); +data =zeros(length(dresp),length(prob)); for i=1:length(dresp) - desc(i)=cellstr(dresp(i).descriptor); - data(i,:)=norminv(plev,dresp(i).mean,dresp(i).stddev); -end - -for j=length(plev):-1:2 + descr(i)=cellstr(dresp(i).descriptor); + data(i,:)=norminv(prob,dresp(i).mean,dresp(i).stddev); +end + +for j=length(prob):-1:2 data(:,j)=data(:,j)-data(:,j-1); end @@ -87,4 +151,5 @@ end +figure hl1=bar(data,'stacked'); % set barseries properties for lowest value @@ -97,17 +162,41 @@ set(ax1,'xlim',[0.5 1.5]) end -set(ax1,'xtick',1:1:max(length(dresp),2)); -set(ax1,'xticklabel',desc) + +ylim('auto') +[ylims]=ylim; +if exist('ymin','var') + ylims(1)=ymin; +end +if exist('ymax','var') + ylims(2)=ymax; +end +ylim(ylims) + +set(ax1,'xtick',1:length(descr)) +set(ax1,'xticklabel',descr) +if exist('xtlrot','var') + htl=rotateticklabel(ax1,xtlrot); + tlext=zeros(length(htl),4); + for i=1:length(htl) + tlext(i,:)=get(htl(i),'Extent'); + end +end % add the annotation -title('Normal Distributions of Variables and/or Responses') -xlabel('Variable or Response') +title('Normal Distributions of Responses') +xlabel('Response') +if exist('xtlrot','var') + xlext=get(get(ax1,'xlabel'),'Extent'); + nskip=ceil(max(tlext(:,4))/xlext(4)); + xlabel(cellstr([repmat(' ',nskip,1);'Response'])); + clear nskip xlext tlext +end ylabel('Value') if ~exist('lstr','var') || isempty(lstr) - lstr=cell(1,length(plev)); - for i=1:length(plev) - lstr(i)=cellstr(sprintf('%g%%',100*plev(i))); + lstr=cell(1,length(prob)); + for i=1:length(prob) + lstr(i)=cellstr(sprintf('%g%%',100*prob(i))); end end Index: /issm/trunk/src/m/solutions/dakota/plot_normplot.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/plot_normplot.m (revision 3091) +++ /issm/trunk/src/m/solutions/dakota/plot_normplot.m (revision 3092) @@ -2,12 +2,27 @@ % plot a normal probability plot of the responses. % -% []=plot_normplot(dresp) +% []=plot_normplot(dresp ,params) +% []=plot_normplot(dresp,descr,params) +% []=plot_normplot(sampr,descr,params) % % where the required input is: % dresp (structure array, responses) +% or +% dresp (structure array, responses) +% descr (cell array, list of response descriptions desired) +% or +% sampr (double array, lists of response samples) +% descr (cell array, list of response descriptions) % % the required fields of dresp are: % descriptor (char, description) % sample (double vector, list of samples) +% +% the optional input is: +% params (string/numeric, parameter names and values) +% +% where the optional parameters are: +% xmin (numeric, minimum of x-axis) +% xmax (numeric, maximum of x-axis) % % for each response in the input array, this function plots @@ -23,5 +38,5 @@ % acknowledged. Any commercial use must be negotiated with % the Office of Technology Transfer at the California Institute -% of Technology. (NTR 47078) +% of Technology. (J. Schiermeier, NTR 47078) % % This software may be subject to U.S. export control laws. @@ -32,5 +47,5 @@ % to foreign countries or providing access to foreign persons." % -function []=plot_normplot(dresp) +function []=plot_normplot(varargin) if ~nargin @@ -39,31 +54,94 @@ end -%% assemble the data into a matrix +%% process input data and assemble into matrices as needed -desc=cell (1,length(dresp)); -for i=1:length(dresp) - ldata(i)=length(dresp(i).sample); -end -data=zeros(max(ldata),length(dresp)); +% responses -for i=1:length(dresp) - desc(i)=cellstr(dresp(i).descriptor); - data(1:ldata(i),i)=dresp(i).sample; +iarg=1; +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp=struc_desc(dresp,varargin{iarg}); + iarg=iarg+1; + end + + descr=cell (1,length(dresp)); + lsamp=zeros(1,length(dresp)); + for i=1:length(dresp) + lsamp(i)=length(dresp(i).sample); + end + sampr=zeros(max(lsamp),length(dresp)); + sampr(:,:)=NaN; + + for i=1:length(dresp) + descr(i)=cellstr(dresp(i).descriptor); + sampr(1:lsamp(i),i)=dresp(i).sample; + end +else + sampr=varargin{iarg}; + iarg=iarg+1; + + if iarg <= nargin && iscell(varargin{iarg}) + descr=varargin{iarg}; + iarg=iarg+1; +% elseif iarg <= nargin && ischar(varargin{iarg}) +% descr=cellstr(varargin{iarg}); +% iarg=iarg+1; + else + descr=cell(1,size(sampr,2)); + end end +for i=1:length(descr) + if isempty(descr{i}) + descr(i)={['resp_' num2str(i)]}; + end +end + +% parameters + +while (iarg <= nargin-1) + if ischar(varargin{iarg}) + eval([varargin{iarg} '=varargin{iarg+1};']); + disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']); + else + error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']); + end + iarg=iarg+2; +end + %% draw the plot % draw normal probability plot -normplot(data) +figure +normplot(sampr) ax1=gca; % add the annotation -title('Normal Probability Plot of Variables and/or Responses') -xlabel('Value') -ylabel('Probability') +xlim('auto') +[xlims]=xlim; +if exist('xmin','var') + xlims(1)=xmin; +end +if exist('xmax','var') + xlims(2)=xmax; +end +xlim(xlims) -hleg1=legend(ax1,desc,'Location','EastOutside',... +if (size(sampr,2) == 1) + tlabc=descr{1}; +else + tlabc='Responses'; +end +title(['Normal Probability Plot of ' tlabc],'Interpreter','none'); +xlabel('Value' ,'Interpreter','none'); +ylabel('Probability','Interpreter','none'); + +hleg1=legend(ax1,descr,'Location','EastOutside',... 'Orientation','vertical','Interpreter','none'); Index: /issm/trunk/src/m/solutions/dakota/plot_prob_bars.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/plot_prob_bars.m (revision 3091) +++ /issm/trunk/src/m/solutions/dakota/plot_prob_bars.m (revision 3092) @@ -2,15 +2,26 @@ % plot a stacked bar chart of the probabilities in the CDF. % -% []=plot_prob_bars(dresp,lstr) +% []=plot_prob_bars(dresp ,params) +% []=plot_prob_bars(dresp,descr,params) % % where the required input is: % dresp (structure array, responses) -% -% and the optional input is: -% lstr (cell array, legend labels) +% or +% dresp (structure array, responses) +% descr (cell array, list of response descriptions desired) % % the required fields of dresp are: % descriptor (char, description) % cdf(:,4) (double matrix, CDF table) +% +% the optional input is: +% params (string/numeric, parameter names and values) +% +% where the optional parameters are: +% ymin (numeric, minimum of y-axis) +% ymax (numeric, maximum of y-axis) +% xtlrot (numeric, rotation in degrees of x-tick labels) +% cmap (char, 'stoplight' for 6-color stoplight colormap) +% lstr (cell array, legend labels) % % for each response in the input array, this function plots @@ -28,5 +39,5 @@ % acknowledged. Any commercial use must be negotiated with % the Office of Technology Transfer at the California Institute -% of Technology. (NTR 47078) +% of Technology. (J. Schiermeier, NTR 47078) % % This software may be subject to U.S. export control laws. @@ -37,5 +48,5 @@ % to foreign countries or providing access to foreign persons." % -function []=plot_prob_bars(dresp,lstr) +function []=plot_prob_bars(varargin) if ~nargin @@ -46,25 +57,57 @@ %% assemble the data into a matrix and calculate the increments -desc=cell (1,length(dresp)); -for i=1:length(dresp) - ldata(i)=size(dresp(i).cdf,1); -end -data=zeros(length(dresp),max(ldata)); -for i=1:length(dresp) - desc(i)=cellstr(dresp(i).descriptor); - if ~isempty(dresp(i).cdf) - data(i,1)=dresp(i).cdf(1,2); - for j=2:size(dresp(i).cdf,1) - if (dresp(i).cdf(j,2) > dresp(i).cdf(j-1,2)) - data(i,j)=dresp(i).cdf(j,2)-dresp(i).cdf(j-1,2); +%% process input data and assemble into matrices and increments + +% responses + +iarg=1; +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp=struc_desc(dresp,varargin{iarg}); + iarg=iarg+1; + end + + descr=cell (1,length(dresp)); + lcdfr=zeros(1,length(dresp)); + for i=1:length(dresp) + lcdfr(i)=size(dresp(i).cdf,1); + end + cdfr=zeros(length(dresp),max(lcdfr)); + + for i=1:length(dresp) + descr(i)=cellstr(dresp(i).descriptor); + if ~isempty(dresp(i).cdf) + cdfr(i,1)=dresp(i).cdf(1,2); + for j=2:size(dresp(i).cdf,1) + if (dresp(i).cdf(j,2) > dresp(i).cdf(j-1,2)) + cdfr(i,j)=dresp(i).cdf(j,2)-dresp(i).cdf(j-1,2); + end end end end +else + error(['''' inputname(iarg) ''' is not a structure.']); end % convert to percentage -data=data*100.; +cdfr=cdfr*100.; + +% parameters + +while (iarg <= nargin-1) + if ischar(varargin{iarg}) + eval([varargin{iarg} '=varargin{iarg+1};']); + disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']); + else + error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']); + end + iarg=iarg+2; +end %% draw the stacked bar plot @@ -74,21 +117,49 @@ if length(dresp) == 1 - data=[data; data]; + cdfr=[cdfr; cdfr]; end -hl1=bar(data,'stacked'); -% set(hl1(1),'FaceColor','green') -% set(hl1(2),'FaceColor','blue') -% set(hl1(3),'FaceColor','yellow') -% set(hl1(4),'FaceColor','cyan') -% set(hl1(5),'FaceColor','magenta') -% set(hl1(6),'FaceColor','red') +figure +hl1=bar(cdfr,'stacked'); +if exist('cmap','var') + if strncmpi(cmap,'stop',4) +% set(hl1(1),'FaceColor','green') +% set(hl1(2),'FaceColor',[.7 1 0]) +% set(hl1(3),'FaceColor','yellow') +% set(hl1(4),'FaceColor',[1 .7 0]) +% set(hl1(5),'FaceColor',[1 .5 0]) +% set(hl1(6),'FaceColor','red') + colormap([0 1 0;.7 1 0;1 1 0;1 .7 0;1 .5 0;1 0 0]); + else + colormap(cmap); + end +end ax1=gca; -set(ax1,'ylim',[0 120]) if length(dresp) == 1 set(ax1,'xlim',[0.5 1.5]) end -set(ax1,'xticklabel',desc) + +% set(ax1,'ylim',[0 120]) +% ylim('auto') +% [ylims]=ylim; +[ylims]=[0 120]; +if exist('ymin','var') + ylims(1)=ymin; +end +if exist('ymax','var') + ylims(2)=ymax; +end +ylim(ylims) + +set(ax1,'xtick',1:length(descr)) +set(ax1,'xticklabel',descr) +if exist('xtlrot','var') + htl=rotateticklabel(ax1,xtlrot); + tlext=zeros(length(htl),4); + for i=1:length(htl) + tlext(i,:)=get(htl(i),'Extent'); + end +end % add the annotation @@ -96,13 +167,19 @@ title('Probabilities for Specified Response Levels (RIA)') xlabel('Response') +if exist('xtlrot','var') + xlext=get(get(ax1,'xlabel'),'Extent'); + nskip=ceil(max(tlext(:,4))/xlext(4)); + xlabel(cellstr([repmat(' ',nskip,1);'Response'])); + clear nskip xlext tlext +end ylabel('Percent Below Level') if ~exist('lstr','var') || isempty(lstr) - lstr=cell(1,max(ldata)); - for i=1:max(ldata) + lstr=cell(1,max(lcdfr)); + for i=1:max(lcdfr) lstr(i)=cellstr(sprintf('%g',... - dresp(find(ldata == max(ldata),1,'first')).cdf(i,1))); + dresp(find(lcdfr == max(lcdfr),1,'first')).cdf(i,1))); end - if ~isempty(find(ldata < max(ldata))) + if ~isempty(find(lcdfr < max(lcdfr))) warning('Variable number of levels for responses.'); end Index: /issm/trunk/src/m/solutions/dakota/plot_rlev_bars.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/plot_rlev_bars.m (revision 3091) +++ /issm/trunk/src/m/solutions/dakota/plot_rlev_bars.m (revision 3092) @@ -2,15 +2,25 @@ % plot a stacked bar chart of the response levels in the cdf. % -% []=plot_rlev_bars(dresp,lstr) +% []=plot_rlev_bars(dresp ,params) +% []=plot_rlev_bars(dresp,descr,params) % % where the required input is: % dresp (structure array, responses) -% -% and the optional input is: -% lstr (cell array, legend labels) +% or +% dresp (structure array, responses) +% descr (cell array, list of response descriptions desired) % % the required fields of dresp are: % descriptor (char, description) % cdf(:,4) (double matrix, CDF table) +% +% the optional input is: +% params (string/numeric, parameter names and values) +% +% where the optional parameters are: +% ymin (numeric, minimum of y-axis) +% ymax (numeric, maximum of y-axis) +% xtlrot (numeric, rotation in degrees of x-tick labels) +% lstr (cell array, legend labels) % % for each response in the input array, this function plots @@ -28,5 +38,5 @@ % acknowledged. Any commercial use must be negotiated with % the Office of Technology Transfer at the California Institute -% of Technology. (NTR 47078) +% of Technology. (J. Schiermeier, NTR 47078) % % This software may be subject to U.S. export control laws. @@ -37,5 +47,5 @@ % to foreign countries or providing access to foreign persons." % -function []=plot_rlev_bars(dresp,lstr) +function []=plot_rlev_bars(varargin) if ~nargin @@ -44,22 +54,51 @@ end -%% assemble the data into a matrix and calculate the increments +%% process input data and assemble into matrices and increments -desc=cell (1,length(dresp)); -for i=1:length(dresp) - ldata(i)=size(dresp(i).cdf,1); -end -data=zeros(length(dresp),max(ldata)); +% responses -for i=1:length(dresp) - desc(i)=cellstr(dresp(i).descriptor); - if ~isempty(dresp(i).cdf) - data(i,1)=dresp(i).cdf(1,1); - for j=2:size(dresp(i).cdf,1) - if (dresp(i).cdf(j,1) > dresp(i).cdf(j-1,1)) - data(i,j)=dresp(i).cdf(j,1)-dresp(i).cdf(j-1,1); +iarg=1; +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp=struc_desc(dresp,varargin{iarg}); + iarg=iarg+1; + end + + descr=cell (1,length(dresp)); + lcdfr=zeros(1,length(dresp)); + for i=1:length(dresp) + lcdfr(i)=size(dresp(i).cdf,1); + end + cdfr=zeros(length(dresp),max(lcdfr)); + + for i=1:length(dresp) + descr(i)=cellstr(dresp(i).descriptor); + if ~isempty(dresp(i).cdf) + cdfr(i,1)=dresp(i).cdf(1,1); + for j=2:size(dresp(i).cdf,1) + if (dresp(i).cdf(j,1) > dresp(i).cdf(j-1,1)) + cdfr(i,j)=dresp(i).cdf(j,1)-dresp(i).cdf(j-1,1); + end end end end +else + error(['''' inputname(iarg) ''' is not a structure.']); +end + +% parameters + +while (iarg <= nargin-1) + if ischar(varargin{iarg}) + eval([varargin{iarg} '=varargin{iarg+1};']); + disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']); + else + error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']); + end + iarg=iarg+2; end @@ -70,8 +109,9 @@ if length(dresp) == 1 - data=[data; data]; + cdfr=[cdfr; cdfr]; end -hl1=bar(data,'stacked'); +figure +hl1=bar(cdfr,'stacked'); % set barseries properties for lowest value whitebg('white') @@ -83,19 +123,44 @@ set(ax1,'xlim',[0.5 1.5]) end -set(ax1,'xticklabel',desc) + +ylim('auto') +[ylims]=ylim; +if exist('ymin','var') + ylims(1)=ymin; +end +if exist('ymax','var') + ylims(2)=ymax; +end +ylim(ylims) + +set(ax1,'xtick',1:length(descr)) +set(ax1,'xticklabel',descr) +if exist('xtlrot','var') + htl=rotateticklabel(ax1,xtlrot); + tlext=zeros(length(htl),4); + for i=1:length(htl) + tlext(i,:)=get(htl(i),'Extent'); + end +end % add the annotation title('Response Levels for Specified Probabilities (PMA)') -xlabel('Response') +xlabel('Response'); +if exist('xtlrot','var') + xlext=get(get(ax1,'xlabel'),'Extent'); + nskip=ceil(max(tlext(:,4))/xlext(4)); + xlabel(cellstr([repmat(' ',nskip,1);'Response'])); + clear nskip xlext tlext +end ylabel('Response Level') if ~exist('lstr','var') || isempty(lstr) - lstr=cell(1,max(ldata)); - for i=1:max(ldata) + lstr=cell(1,max(lcdfr)); + for i=1:max(lcdfr) lstr(i)=cellstr(sprintf('%g%%',... - 100*dresp(find(ldata == max(ldata),1,'first')).cdf(i,2))); + 100*dresp(find(lcdfr == max(lcdfr),1,'first')).cdf(i,2))); end - if ~isempty(find(ldata < max(ldata))) + if ~isempty(find(lcdfr < max(lcdfr),1)) warning('Variable number of probabilities for responses.'); end Index: /issm/trunk/src/m/solutions/dakota/plot_rlev_bars_ci.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/plot_rlev_bars_ci.m (revision 3092) +++ /issm/trunk/src/m/solutions/dakota/plot_rlev_bars_ci.m (revision 3092) @@ -0,0 +1,263 @@ +% +% plot a stacked bar chart of the response levels in the cdf +% for the sample and confidence intervals. +% +% []=plot_rlev_bars_ci(dresp ,params) +% []=plot_rlev_bars_ci(dresp,descr,params) +% []=plot_rlev_bars_ci(sampr,descr,params) +% +% where the required input is: +% dresp (structure array, responses) +% or +% dresp (structure array, responses) +% descr (cell array, list of response descriptions desired) +% or +% sampr (double array, lists of response samples) +% descr (cell array, list of response descriptions) +% +% the required fields of dresp are: +% descriptor (char, description) +% cdf(:,4) (double matrix, CDF table) +% +% and the optional fields of dresp are: +% mean (double, mean of sample) +% stddev (double, standard deviation of sample) +% meanci(2) (double, confidence interval of mean) +% stddevci(2) (double, confidence interval of standard deviation) +% +% the optional input is: +% params (string/numeric, parameter names and values) +% +% where the optional parameters are: +% ymin (numeric, minimum of y-axis) +% ymax (numeric, maximum of y-axis) +% xtlrot (numeric, rotation in degrees of x-tick labels) +% lstr (cell array, legend labels) +% +% for each response in the input array, this function plots +% a stacked bar plot of the responses, where the bars are +% stacked by the response levels corresponding to the given +% probabilities in the CDF, and annotates it with the +% description. the response levels for the normal distribution +% and the confidence intervals are also plotted. the legend +% labels can be given or constructed from the probabilities. +% +% dresp data would typically be contained in the dakota tabular +% output file from a sampling analysis, read by dakota_out_parse. +% +% "Copyright 2009, by the California Institute of Technology. +% ALL RIGHTS RESERVED. United States Government Sponsorship +% acknowledged. Any commercial use must be negotiated with +% the Office of Technology Transfer at the California Institute +% of Technology. (J. Schiermeier, NTR 47078) +% +% This software may be subject to U.S. export control laws. +% By accepting this software, the user agrees to comply with +% all applicable U.S. export laws and regulations. User has the +% responsibility to obtain export licenses, or other export +% authority as may be required before exporting such information +% to foreign countries or providing access to foreign persons." +% +function []=plot_rlev_bars_ci(varargin) + +if ~nargin + help plot_rlev_bars_ci + return +end + +%% process input data and assemble into dresp as needed + +% responses + +iarg=1; +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp=struc_desc(dresp,varargin{iarg}); + iarg=iarg+1; + end +else + sampr=varargin{iarg}; + iarg=iarg+1; + + if iarg <= nargin && iscell(varargin{iarg}) + descr=varargin{iarg}; + iarg=iarg+1; +% elseif iarg <= nargin && ischar(varargin{iarg}) +% descr=cellstr(varargin{iarg}); +% iarg=iarg+1; + else + descr=cell(1:size(sampr,2)); + end + + dresp=struct([]); + for i=1:size(sampr,2) + dresp(end+1).sample=sampr(:,i); + if ~isempty(descr) + dresp(i).descriptor=descr{i}; + else + dresp(i).descriptor=['dresp_' num2str(i)]; + end + end +end + +% parameters + +while (iarg <= nargin-1) + if ischar(varargin{iarg}) + eval([varargin{iarg} '=varargin{iarg+1};']); + disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']); + else + error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']); + end + iarg=iarg+2; +end + +%% calculate any missing information (noting that dresp is local) + +for i=1:length(dresp) + if ~isfield(dresp(i),'mean') || isempty(dresp(i).mean) || ... + ~isfield(dresp(i),'stddev') || isempty(dresp(i).stddev) || ... + ~isfield(dresp(i),'meanci') || isempty(dresp(i).meanci) || ... + ~isfield(dresp(i),'stddevci') || isempty(dresp(i).stddevci) +% calculate 95% confidence intervals (same as dakota) + [dresp(i).mean,dresp(i).stddev,... + dresp(i).meanci,dresp(i).stddevci]=... + normfit(sampr(:,i),0.05); + display('Using calculated normal fits from sample data.') + end + + if ~isfield(dresp(i),'cdf') || isempty(dresp(i).cdf) +% use minus/plus integer standard deviations + sdvect=[-4 -3 -2 -1 0 1 2 3 4]; + dresp(i).cdf(:,2)=normcdf(sdvect,0,1); + dresp(i).cdf(:,1)=norminv(dresp(i).cdf(:,2),... + dresp(i).mean,dresp(i).stddev); + display('Using integer standard deviations for percentages.') + + if ~exist('lstr','var') || isempty(lstr) + lstr=cell(1,size(dresp(i).cdf,1)); + for j=1:size(dresp(i).cdf,1) + if sdvect(j) + lstr{j}=sprintf('mu %+d sigma',sdvect(j)); + else + lstr{j}='mu'; + end + end + end + end +end + +%% assemble the data into a matrix and calculate the increments + +descr=cell (1,0); +lcdfr=zeros(1,length(dresp)); +for i=1:length(dresp) + lcdfr(i)=size(dresp(i).cdf,1); +end +cdfr=zeros(0,max(lcdfr)); + +% fill in the cdf data + +for i=1:length(dresp) + if ~isempty(dresp(i).cdf) + descr(end+1)=cellstr([dresp(i).descriptor]); + cdfr(end+1,:)=dresp(i).cdf(:,1); + if isfield(dresp(i),'mean' ) && ~isempty(dresp(i).mean ) && ... + isfield(dresp(i),'stddev') && ~isempty(dresp(i).stddev) + descr(end+1)=cellstr([dresp(i).descriptor ' norm']); + cdfr(end+1,:)=norminv(dresp(i).cdf(:,2),dresp(i).mean,dresp(i).stddev); + end + if isfield(dresp(i),'meanci' ) && ~isempty(dresp(i).meanci ) && ... + isfield(dresp(i),'stddevci') && ~isempty(dresp(i).stddevci) + descr(end+1)=cellstr([dresp(i).descriptor ' norm-+']); + descr(end+1)=cellstr([dresp(i).descriptor ' norm--']); + descr(end+1)=cellstr([dresp(i).descriptor ' norm+-']); + descr(end+1)=cellstr([dresp(i).descriptor ' norm++']); + cdfr(end+1,:)=norminv(dresp(i).cdf(:,2),dresp(i).meanci(1),dresp(i).stddevci(2)); + cdfr(end+1,:)=norminv(dresp(i).cdf(:,2),dresp(i).meanci(1),dresp(i).stddevci(1)); + cdfr(end+1,:)=norminv(dresp(i).cdf(:,2),dresp(i).meanci(2),dresp(i).stddevci(1)); + cdfr(end+1,:)=norminv(dresp(i).cdf(:,2),dresp(i).meanci(2),dresp(i).stddevci(2)); + end + end +end + +% calculate the increments + +for i=1:size(cdfr,1) + for j=find(cdfr(i,:),1,'last'):-1:2 + cdfr(i,j)=cdfr(i,j)-cdfr(i,j-1); + end +end + +%% draw the stacked bar plot + +% if there's only one row, Matlab 7.5 interprets it as a column, +% so add an extra row, then reduce xlim + +if length(descr) == 1 + cdfr=[cdfr; cdfr]; +end + +figure +hl1=bar(cdfr,'stacked'); +% set barseries properties for lowest value +whitebg('white') +set(hl1(1),'FaceColor','white') +set(hl1(1),'Visible','off') + +ax1=gca; +if length(descr) == 1 + set(ax1,'xlim',[0.5 1.5]) +end + +ylim('auto') +[ylims]=ylim; +if exist('ymin','var') + ylims(1)=ymin; +end +if exist('ymax','var') + ylims(2)=ymax; +end +ylim(ylims) + +set(ax1,'xtick',1:length(descr)) +set(ax1,'xticklabel',descr) +if exist('xtlrot','var') + htl=rotateticklabel(ax1,xtlrot); + tlext=zeros(length(htl),4); + for i=1:length(htl) + tlext(i,:)=get(htl(i),'Extent'); + end +end + +% add the annotation + +title('Response Levels for Specified Probabilities (PMA)'); +xlabel('Response'); +if exist('xtlrot','var') + xlext=get(get(ax1,'xlabel'),'Extent'); + nskip=ceil(max(tlext(:,4))/xlext(4)); + xlabel(cellstr([repmat(' ',nskip,1);'Response'])); + clear nskip xlext tlext +end +ylabel('Response Level'); + +if ~exist('lstr','var') || isempty(lstr) + lstr=cell(1,max(lcdfr)); + for i=1:max(lcdfr) + lstr(i)=cellstr(sprintf('%g%%',... + 100*dresp(find(lcdfr == max(lcdfr),1,'first')).cdf(i,2))); + end + if ~isempty(find(lcdfr < max(lcdfr),1,'first')) + warning('Variable number of probabilities for responses.'); + end +end + +hleg1=legend(ax1,lstr,'Location','EastOutside',... + 'Orientation','vertical','Interpreter','none'); + +end Index: /issm/trunk/src/m/solutions/dakota/plot_rvsv_line.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/plot_rvsv_line.m (revision 3092) +++ /issm/trunk/src/m/solutions/dakota/plot_rvsv_line.m (revision 3092) @@ -0,0 +1,257 @@ +% +% plot line plots of responses vs. variables. +% +% []=plot_rvsv_line(dvar ,dresp ,params) +% []=plot_rvsv_line(dvar ,descv,dresp,descr,params) +% []=plot_rvsv_line(sampv,descv,sampr,descr,params) +% +% where the required input is: +% dvar (structure array, variables) +% or +% dvar (structure array, variables) +% descv (cell array, list of variable descriptions desired) +% or +% sampv (double array, lists of variable samples) +% descv (cell array, list of variable descriptions) +% +% dresp (structure array, responses) +% or +% dresp (structure array, responses) +% descr (cell array, list of response descriptions desired) +% or +% sampr (double array, lists of response samples) +% descr (cell array, list of response descriptions) +% +% the required fields of dvar and dresp are: +% descriptor (char, description) +% sample (double vector, list of samples) +% +% the optional input is: +% params (string/numeric, parameter names and values) +% +% where the optional parameters are: +% nplotr (numeric, number of plot rows) +% nplotc (numeric, number of plot columns) +% ymin (numeric, minimum of y-axis) +% ymax (numeric, maximum of y-axis) +% yscat (char, 'off' to turn off y-axis scattergram) +% +% for each variable/response combination in the input array, this +% function plots a line plot. all of the variables and responses +% are plotted on the same axes, if nplotr and nplotc are not +% specified, so some scaling might otherwise be desired. +% +% dvar and dresp data would typically be contained in the dakota +% tabular output file from a sampling or parametric analysis, and +% read by dakota_out_parse. +% +% "Copyright 2009, by the California Institute of Technology. +% ALL RIGHTS RESERVED. United States Government Sponsorship +% acknowledged. Any commercial use must be negotiated with +% the Office of Technology Transfer at the California Institute +% of Technology. (J. Schiermeier, NTR 47078) +% +% This software may be subject to U.S. export control laws. +% By accepting this software, the user agrees to comply with +% all applicable U.S. export laws and regulations. User has the +% responsibility to obtain export licenses, or other export +% authority as may be required before exporting such information +% to foreign countries or providing access to foreign persons." +% +function []=plot_rvsv_line(varargin) + +if ~nargin + help plot_rvsv_line + return +end + +%% process input data and assemble into matrices as needed + +% variables + +iarg=1; +if isstruct(varargin{iarg}) + dvar=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dvar=struc_desc(dvar,varargin{iarg}); + iarg=iarg+1; + end + + descv=cell (1,length(dvar)); + lsamp=zeros(1,length(dvar)); + for i=1:length(dvar) + lsamp(i)=length(dvar(i).sample); + end + sampv=zeros(max(lsamp),length(dvar)); + sampv(:,:)=NaN; + + for i=1:length(dvar) + descv(i)=cellstr(dvar(i).descriptor); + sampv(1:lsamp(i),i)=dvar(i).sample; + end +else + sampv=varargin{iarg}; + iarg=iarg+1; + + if iarg <= nargin && iscell(varargin{iarg}) + descv=varargin{iarg}; + iarg=iarg+1; +% elseif iarg <= nargin && ischar(varargin{iarg}) +% descv=cellstr(varargin{iarg}); +% iarg=iarg+1; + else + descv=cell(1,size(sampv,2)); + end +end + +for i=1:length(descv) + if isempty(descv{i}) + descv(i)={['var_' i]}; + end +end + +% responses + +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp=struc_desc(dresp,varargin{iarg}); + iarg=iarg+1; + end + + descr=cell (1,length(dresp)); + lsamp=zeros(1,length(dresp)); + for i=1:length(dresp) + lsamp(i)=length(dresp(i).sample); + end + sampr=zeros(max(lsamp),length(dresp)); + sampr(:,:)=NaN; + + for i=1:length(dresp) + descr(i)=cellstr(dresp(i).descriptor); + sampr(1:lsamp(i),i)=dresp(i).sample; + end +else + sampr=varargin{iarg}; + iarg=iarg+1; + + if iarg <= nargin && iscell(varargin{iarg}) + descr=varargin{iarg}; + iarg=iarg+1; +% elseif iarg <= nargin && ischar(varargin{iarg}) +% descr=cellstr(varargin{iarg}); +% iarg=iarg+1; + else + descr=cell(1,size(sampr,2)); + end +end + +for i=1:length(descr) + if isempty(descr{i}) + descr(i)={['resp_' num2str(i)]}; + end +end + +% parameters + +while (iarg <= nargin-1) + if ischar(varargin{iarg}) + eval([varargin{iarg} '=varargin{iarg+1};']); + disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']); + else + error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']); + end + iarg=iarg+2; +end + +if ~exist('nplotr','var') && ~exist('nplotc','var') + nplotr=1; + nplotc=1; +elseif ~exist('nplotr','var') + nplotr=ceil(size(sampr,2)*size(sampv,2)/nplotc); +elseif ~exist('nplotc','var') + nplotc=ceil(size(sampr,2)*size(sampv,2)/nplotr); +end + +%% filter, sort, and plot the data + +% while it would be preferable for the outer loop to be responses, +% it is more efficient for the outer loop to be variables + +figure +haxes=[]; +hplot=[]; +cdesc={}; +hscat=[]; + +iplot=0; + +for ivar=1:size(sampv,2) + [vval,indxv,indxvi]=unique(sampv(:,ivar),'first'); + indxv2=setdiff(1:size(sampv,1),indxv); + + for iresp=1:size(sampr,2) + +% initialize the subplot + + if (ivar*iresp == 1) || ... + (nplotr*nplotc > 1) + iplot=iplot+1; + haxes(end+1)=subplot(nplotr,nplotc,iplot); + hold all + end + + hplot(end+1)=plot (sampv(indxv ,ivar),sampr(indxv ,iresp),'-x'); + cdesc(end+1)={[descr{iresp} ' wrt ' descv{ivar}]}; + if ~exist('yscat','var') || ~strcmpi(yscat,'off') + hscat(end+1)=scatter(sampv(indxv2,ivar),sampr(indxv2,iresp),'+k'); +% see "controlling legends" in Matlab on-line docs +% cdesc(end+1)={['constant ' descv{ivar}]}; + set(get(get(hscat(end),'Annotation'),'LegendInformation'),... + 'IconDisplayStyle','off'); % Exclude line from legend + end + +% add the annotation + + if (ivar*iresp == size(sampv,2)*size(sampr,2)) || ... + (nplotr*nplotc > 1) + hold off + + ylim('auto') + [ylims]=ylim; + if exist('ymin','var') + ylims(1)=ymin; + end + if exist('ymax','var') + ylims(2)=ymax; + end + ylim(ylims) + + if (size(sampv,2) == 1) || (nplotr*nplotc > 1) + xlabc=descv{ivar}; + else + xlabc='Variables'; + end + if (size(sampr,2) == 1) || (nplotr*nplotc > 1) + ylabc=descr{iresp}; + else + ylabc='Responses'; + end + title([ylabc ' vs. ' xlabc],'Interpreter','none'); + xlabel(xlabc,'Interpreter','none'); + ylabel(ylabc,'Interpreter','none'); + + if (nplotr*nplotc == 1) && (size(sampv,2)*size(sampr,2) > 1) + legend(cdesc,'Location','EastOutside','Interpreter','none'); + end + end + end +end + +end Index: /issm/trunk/src/m/solutions/dakota/plot_rvsv_scat3.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/plot_rvsv_scat3.m (revision 3092) +++ /issm/trunk/src/m/solutions/dakota/plot_rvsv_scat3.m (revision 3092) @@ -0,0 +1,273 @@ +% +% plot 3-d scatter plots of variables vs. responses. +% +% []=plot_rvsv_scat3(dvar ,dresp ,params) +% []=plot_rvsv_scat3(dvar ,descv,dresp,descr,params) +% []=plot_rvsv_scat3(sampv,descv,sampr,descr,params) +% +% where the required input is: +% dvar (structure array, variables) +% or +% dvar (structure array, variables) +% descv (cell array, list of variable descriptions desired) +% or +% sampv (double array, lists of variable samples) +% descv (cell array, list of variable descriptions) +% +% dresp (structure array, responses) +% or +% dresp (structure array, responses) +% descr (cell array, list of response descriptions desired) +% or +% sampr (double array, lists of response samples) +% descr (cell array, list of response descriptions) +% +% the required fields of dvar and dresp are: +% descriptor (char, description) +% sample (double vector, list of samples) +% +% the optional input is: +% params (string/numeric, parameter names and values) +% +% where the optional parameters are: +% nplotr (numeric, number of plot rows) +% nplotc (numeric, number of plot columns) +% ymin (numeric, minimum of y-axis) +% ymax (numeric, maximum of y-axis) +% zmin (numeric, minimum of z-axis) +% zmax (numeric, maximum of z-axis) +% cmin (numeric, minimum of colorbar) +% cmax (numeric, maximum of colorbar) +% smark (numeric, size of markers) +% +% for each response in the input array, this function plots a 3-d +% scatter plot. there should be no more than three variables. +% each response will be in a separate scatter plot; hence the +% need for nplotr and nplotc. +% +% dvar and dresp data would typically be contained in the dakota +% tabular output file from a sampling or parametric analysis, and +% read by dakota_out_parse. +% +% "Copyright 2009, by the California Institute of Technology. +% ALL RIGHTS RESERVED. United States Government Sponsorship +% acknowledged. Any commercial use must be negotiated with +% the Office of Technology Transfer at the California Institute +% of Technology. (J. Schiermeier, NTR 47078) +% +% This software may be subject to U.S. export control laws. +% By accepting this software, the user agrees to comply with +% all applicable U.S. export laws and regulations. User has the +% responsibility to obtain export licenses, or other export +% authority as may be required before exporting such information +% to foreign countries or providing access to foreign persons." +% +function []=plot_rvsv_scat3(varargin) + +if ~nargin + help plot_rvsv_scat3 + return +end + +%% process input data and assemble into matrices as needed + +iarg=1; + +% variables + +if isstruct(varargin{iarg}) + dvar=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dvar=struc_desc(dvar,varargin{iarg}); + iarg=iarg+1; + end + + descv=cell (1,length(dvar)); + lsamp=zeros(1,length(dvar)); + for i=1:length(dvar) + lsamp(i)=length(dvar(i).sample); + end + sampv=zeros(max(lsamp),length(dvar)); + sampv(:,:)=NaN; + + for i=1:length(dvar) + descv(i)=cellstr(dvar(i).descriptor); + sampv(1:lsamp(i),i)=dvar(i).sample; + end +else + sampv=varargin{iarg}; + iarg=iarg+1; + + if iarg <= nargin && iscell(varargin{iarg}) + descv=varargin{iarg}; + iarg=iarg+1; +% elseif iarg <= nargin && ischar(varargin{iarg}) +% descv=cellstr(varargin{iarg}); +% iarg=iarg+1; + else + descv=cell(1,size(sampv,2)); + end +end + +for i=1:length(descv) + if isempty(descv{i}) + descv(i)={['var_' i]}; + end +end + +if (size(sampv,2) > 3) + error('No more than three variables required for 3-d scatter plot.'); +end + +% responses + +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp=struc_desc(dresp,varargin{iarg}); + iarg=iarg+1; + end + + descr=cell (1,length(dresp)); + lsamp=zeros(1,length(dresp)); + for i=1:length(dresp) + lsamp(i)=length(dresp(i).sample); + end + sampr=zeros(max(lsamp),length(dresp)); + sampr(:,:)=NaN; + + for i=1:length(dresp) + descr(i)=cellstr(dresp(i).descriptor); + sampr(1:lsamp(i),i)=dresp(i).sample; + end +else + sampr=varargin{iarg}; + iarg=iarg+1; + + if iarg <= nargin && iscell(varargin{iarg}) + descr=varargin{iarg}; + iarg=iarg+1; +% elseif iarg <= nargin && ischar(varargin{iarg}) +% descr=cellstr(varargin{iarg}); +% iarg=iarg+1; + else + descr=cell(1,size(sampr,2)); + end +end + +for i=1:length(descr) + if isempty(descr{i}) + descr(i)={['resp_' num2str(i)]}; + end +end + +% parameters + +while (iarg <= nargin-1) + if ischar(varargin{iarg}) + eval([varargin{iarg} '=varargin{iarg+1};']); + disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']); + else + error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']); + end + iarg=iarg+2; +end + +if ~exist('nplotr','var') && ~exist('nplotc','var') + nplotr=ceil(sqrt(size(sampr,2))); + nplotc=ceil(size(sampr,2)/nplotr); +elseif ~exist('nplotr','var') + nplotr=ceil(size(sampr,2)/nplotc); +elseif ~exist('nplotc','var') + nplotc=ceil(size(sampr,2)/nplotr); +end + +if ~exist('smark','var') + smark=100; +end + +%% filter, sort, and plot the data + +figure +haxes =[]; +hscat3=[]; + +for iresp=1:size(sampr,2) + +% initialize the subplot + + haxes(end+1)=subplot(nplotr,nplotc,iresp); + switch size(sampv,2) + case 1 + hscat3(end+1)=scatter (sampv(:,1),sampr(:,iresp),... + smark,sampr(:,iresp),'filled'); + case 2 + hscat3(end+1)=scatter3(sampv(:,1),sampv(:,2),sampr(:,iresp),... + smark,sampr(:,iresp),'filled'); + case 3 + hscat3(end+1)=scatter3(sampv(:,1),sampv(:,2),sampv(:,3),... + smark,sampr(:,iresp),'filled'); + end + + ylim('auto') + [ylims]=ylim; + if exist('ymin','var') + ylims(1)=ymin; + end + if exist('ymax','var') + ylims(2)=ymax; + end + ylim(ylims) + + zlim('auto') + [zlims]=zlim; + if exist('zmin','var') + zlims(1)=zmin; + end + if exist('zmax','var') + zlims(2)=zmax; + end + zlim(zlims) + +% add the annotation + + switch size(sampv,2) + case 1 + title([descr{iresp} ' wrt ' descv{1}],... + 'Interpreter','none'); + xlabel(descv{1},'Interpreter','none'); + ylabel(descr{iresp},'Interpreter','none'); + case 2 + title([descr{iresp} ' wrt ' descv{1} ' and ' descv{2}],... + 'Interpreter','none'); + xlabel(descv{1},'Interpreter','none'); + ylabel(descv{2},'Interpreter','none'); + zlabel(descr{iresp},'Interpreter','none'); + case 3 + title([descr{iresp} ' wrt ' descv{1} ' and ' descv{2} ' and ' descv{3}],... + 'Interpreter','none'); + xlabel(descv{1},'Interpreter','none'); + ylabel(descv{2},'Interpreter','none'); + zlabel(descv{3},'Interpreter','none'); + end + + caxis('auto') + [cmini,cmaxi]=caxis; + if exist('cmin','var') + cmini=cmin; + end + if exist('cmax','var') + cmaxi=cmax; + end + caxis([cmini cmaxi]) + + colorbar +end + +end Index: /issm/trunk/src/m/solutions/dakota/plot_rvsv_surf.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/plot_rvsv_surf.m (revision 3092) +++ /issm/trunk/src/m/solutions/dakota/plot_rvsv_surf.m (revision 3092) @@ -0,0 +1,240 @@ +% +% plot surface plots of variables vs. responses. +% +% []=plot_rvsv_surf(dvar ,dresp ,params) +% []=plot_rvsv_surf(dvar ,descv,dresp,descr,params) +% []=plot_rvsv_surf(sampv,descv,sampr,descr,params) +% +% where the required input is: +% dvar (structure array, variables) +% or +% dvar (structure array, variables) +% descv (cell array, list of variable descriptions desired) +% or +% sampv (double array, lists of variable samples) +% descv (cell array, list of variable descriptions) +% +% dresp (structure array, responses) +% or +% dresp (structure array, responses) +% descr (cell array, list of response descriptions desired) +% or +% sampr (double array, lists of response samples) +% descr (cell array, list of response descriptions) +% +% the required fields of dvar and dresp are: +% descriptor (char, description) +% sample (double vector, list of samples) +% +% the optional input is: +% params (string/numeric, parameter names and values) +% +% where the optional parameters are: +% nplotr (numeric, number of plot rows) +% nplotc (numeric, number of plot columns) +% zmin (numeric, minimum of z-axis) +% zmax (numeric, maximum of z-axis) +% cmin (numeric, minimum of colorbar) +% cmax (numeric, maximum of colorbar) +% +% for each response in the input array, this function plots a +% surface plot. there should be two and only two variables. +% each response will be in a separate surface plot; hence the +% need for nplotr and nplotc. +% +% dvar and dresp data would typically be contained in the dakota +% tabular output file from a sampling or parametric analysis, and +% read by dakota_out_parse. +% +% "Copyright 2009, by the California Institute of Technology. +% ALL RIGHTS RESERVED. United States Government Sponsorship +% acknowledged. Any commercial use must be negotiated with +% the Office of Technology Transfer at the California Institute +% of Technology. (J. Schiermeier, NTR 47078) +% +% This software may be subject to U.S. export control laws. +% By accepting this software, the user agrees to comply with +% all applicable U.S. export laws and regulations. User has the +% responsibility to obtain export licenses, or other export +% authority as may be required before exporting such information +% to foreign countries or providing access to foreign persons." +% +function []=plot_rvsv_surf(varargin) + +if ~nargin + help plot_rvsv_surf + return +end + +%% process input data and assemble into matrices as needed + +iarg=1; + +% variables + +if isstruct(varargin{iarg}) + dvar=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dvar=struc_desc(dvar,varargin{iarg}); + iarg=iarg+1; + end + + descv=cell (1,length(dvar)); + lsamp=zeros(1,length(dvar)); + for i=1:length(dvar) + lsamp(i)=length(dvar(i).sample); + end + sampv=zeros(max(lsamp),length(dvar)); + sampv(:,:)=NaN; + + for i=1:length(dvar) + descv(i)=cellstr(dvar(i).descriptor); + sampv(1:lsamp(i),i)=dvar(i).sample; + end +else + sampv=varargin{iarg}; + iarg=iarg+1; + + if iarg <= nargin && iscell(varargin{iarg}) + descv=varargin{iarg}; + iarg=iarg+1; +% elseif iarg <= nargin && ischar(varargin{iarg}) +% descv=cellstr(varargin{iarg}); +% iarg=iarg+1; + else + descv=cell(1,size(sampv,2)); + end +end + +for i=1:length(descv) + if isempty(descv{i}) + descv(i)={['var_' i]}; + end +end + +if (size(sampv,2) ~= 2) + error('Two and only two variables required for surface plot.'); +end + +% responses + +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp=struc_desc(dresp,varargin{iarg}); + iarg=iarg+1; + end + + descr=cell (1,length(dresp)); + lsamp=zeros(1,length(dresp)); + for i=1:length(dresp) + lsamp(i)=length(dresp(i).sample); + end + sampr=zeros(max(lsamp),length(dresp)); + sampr(:,:)=NaN; + + for i=1:length(dresp) + descr(i)=cellstr(dresp(i).descriptor); + sampr(1:lsamp(i),i)=dresp(i).sample; + end +else + sampr=varargin{iarg}; + iarg=iarg+1; + + if iarg <= nargin && iscell(varargin{iarg}) + descr=varargin{iarg}; + iarg=iarg+1; +% elseif iarg <= nargin && ischar(varargin{iarg}) +% descr=cellstr(varargin{iarg}); +% iarg=iarg+1; + else + descr=cell(1,size(sampr,2)); + end +end + +for i=1:length(descr) + if isempty(descr{i}) + descr(i)={['resp_' num2str(i)]}; + end +end + +% parameters + +while (iarg <= nargin-1) + if ischar(varargin{iarg}) + eval([varargin{iarg} '=varargin{iarg+1};']); + disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']); + else + error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']); + end + iarg=iarg+2; +end + +if ~exist('nplotr','var') && ~exist('nplotc','var') + nplotr=ceil(sqrt(size(sampr,2))); + nplotc=ceil(size(sampr,2)/nplotr); +elseif ~exist('nplotr','var') + nplotr=ceil(size(sampr,2)/nplotc); +elseif ~exist('nplotc','var') + nplotc=ceil(size(sampr,2)/nplotr); +end + +%% filter, sort, and plot the data + +figure +haxes=[]; +hsurf=[]; + +[x,ix,ixi]=unique(sampv(:,1),'first'); +[y,iy,iyi]=unique(sampv(:,2),'first'); + +for iresp=1:size(sampr,2) + z=zeros(length(x),length(y)); + for i=1:size(sampr,1) + z(ixi(i),iyi(i))=sampr(i,iresp); + end + +% initialize the subplot + + haxes(iresp)=subplot(nplotr,nplotc,iresp); +% hsurf(iresp)=surfc(x,y,z); + surfc(x,y,z); + + zlim('auto') + [zlims]=zlim; + if exist('zmin','var') + zlims(1)=zmin; + end + if exist('zmax','var') + zlims(2)=zmax; + end + zlim(zlims) + +% add the annotation + + title([descr{iresp} ' wrt ' descv{1} ' and ' descv{2}],... + 'Interpreter','none'); + xlabel(descv{1},'Interpreter','none'); + ylabel(descv{2},'Interpreter','none'); + zlabel(descr{iresp},'Interpreter','none'); + + caxis('auto') + [cmini,cmaxi]=caxis; + if exist('cmin','var') + cmini=cmin; + end + if exist('cmax','var') + cmaxi=cmax; + end + caxis([cmini cmaxi]) + + colorbar +end + +end Index: /issm/trunk/src/m/solutions/dakota/plot_sampdist_bars.m =================================================================== --- /issm/trunk/src/m/solutions/dakota/plot_sampdist_bars.m (revision 3091) +++ /issm/trunk/src/m/solutions/dakota/plot_sampdist_bars.m (revision 3092) @@ -2,11 +2,16 @@ % plot a stacked bar chart of the sample distributions. % -% []=plot_sampdist_bars(dresp,lstr) +% []=plot_sampdist_bars(dresp ,params) +% []=plot_sampdist_bars(dresp,descr,params) +% []=plot_sampdist_bars(sampr,descr,params) % % where the required input is: % dresp (structure array, responses) -% -% and the optional input is: -% lstr (cell array, legend labels) +% or +% dresp (structure array, responses) +% descr (cell array, list of response descriptions desired) +% or +% sampr (double array, lists of response samples) +% descr (cell array, list of response descriptions) % % the required fields of dresp are: @@ -21,4 +26,13 @@ % max (double, maximum of sample) % +% the optional input is: +% params (string/numeric, parameter names and values) +% +% where the optional parameters are: +% ymin (numeric, minimum of y-axis) +% ymax (numeric, maximum of y-axis) +% xtlrot (numeric, rotation in degrees of x-tick labels) +% lstr (cell array, legend labels) +% % for each response in the input array, this function plots % a stacked bar plot of the list of samples, where the bars @@ -34,5 +48,5 @@ % acknowledged. Any commercial use must be negotiated with % the Office of Technology Transfer at the California Institute -% of Technology. (NTR 47078) +% of Technology. (J. Schiermeier, NTR 47078) % % This software may be subject to U.S. export control laws. @@ -43,5 +57,5 @@ % to foreign countries or providing access to foreign persons." % -function []=plot_sampdist_bars(dresp,lstr) +function []=plot_sampdist_bars(varargin) if ~nargin @@ -50,8 +64,56 @@ end -%% assemble the data into a matrix and calculate the increments - -desc=cell (1,length(dresp)); -data=zeros(length(dresp),5); +%% process input data and assemble into dresp as needed + +% responses + +iarg=1; +if isstruct(varargin{iarg}) + dresp=varargin{iarg}; + iarg=iarg+1; + +% if iarg <= nargin && (iscell(varargin{iarg}) || ischar(varargin{iarg})) + if iarg <= nargin && iscell(varargin{iarg}) + dresp=struc_desc(dresp,varargin{iarg}); + iarg=iarg+1; + end +else + sampr=varargin{iarg}; + iarg=iarg+1; + + if iarg <= nargin && iscell(varargin{iarg}) + descr=varargin{iarg}; + iarg=iarg+1; +% elseif iarg <= nargin && ischar(varargin{iarg}) +% descr=cellstr(varargin{iarg}); +% iarg=iarg+1; + else + descr=cell(1:size(sampr,2)); + end + + dresp=struct([]); + for i=1:size(sampr,2) + dresp(end+1).sample=sampr(:,i); + if ~isempty(descr) + dresp(i).descriptor=descr{i}; + else + dresp(i).descriptor=['dresp_' num2str(i)]; + end + end +end + +% parameters + +while (iarg <= nargin-1) + if ischar(varargin{iarg}) + eval([varargin{iarg} '=varargin{iarg+1};']); + disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']); + else + error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']); + end + iarg=iarg+2; +end + +%% calculate any missing information (noting that dresp is local) if ~isfield(dresp,'min') || ~isfield(dresp,'quart1') || ... @@ -67,6 +129,11 @@ end +%% assemble the data into a matrix and calculate the increments + +descr=cell (1,length(dresp)); +data =zeros(length(dresp),5); + for i=1:length(dresp) - desc(i)=cellstr(dresp(i).descriptor); + descr(i)=cellstr(dresp(i).descriptor); data(i,1)=dresp(i).min; data(i,2)=dresp(i).quart1-dresp(i).min; @@ -85,4 +152,5 @@ end +figure hl1=bar(data,'stacked'); % set barseries properties for lowest value @@ -95,11 +163,35 @@ set(ax1,'xlim',[0.5 1.5]) end -set(ax1,'xtick',1:1:max(length(dresp),2)); -set(ax1,'xticklabel',desc) + +ylim('auto') +[ylims]=ylim; +if exist('ymin','var') + ylims(1)=ymin; +end +if exist('ymax','var') + ylims(2)=ymax; +end +ylim(ylims) + +set(ax1,'xtick',1:length(descr)) +set(ax1,'xticklabel',descr) +if exist('xtlrot','var') + htl=rotateticklabel(ax1,xtlrot); + tlext=zeros(length(htl),4); + for i=1:length(htl) + tlext(i,:)=get(htl(i),'Extent'); + end +end % add the annotation -title('Sample Distributions of Variables and/or Responses') -xlabel('Variable or Response') +title('Sample Distributions of Responses') +xlabel('Response') +if exist('xtlrot','var') + xlext=get(get(ax1,'xlabel'),'Extent'); + nskip=ceil(max(tlext(:,4))/xlext(4)); + xlabel(cellstr([repmat(' ',nskip,1);'Response'])); + clear nskip xlext tlext +end ylabel('Value')