Index: /issm/trunk/src/m/dakota/dakota_cdfs.m =================================================================== --- /issm/trunk/src/m/dakota/dakota_cdfs.m (revision 4129) +++ /issm/trunk/src/m/dakota/dakota_cdfs.m (revision 4129) @@ -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/dakota/dakota_in_data.m =================================================================== --- /issm/trunk/src/m/dakota/dakota_in_data.m (revision 4129) +++ /issm/trunk/src/m/dakota/dakota_in_data.m (revision 4129) @@ -0,0 +1,137 @@ +% +% define the data to write the dakota .in and .m files. +% +% []=dakota_in_data(dmeth,variables,responses,dparams,filei,varargin) +% +% where the required input is: +% dmeth (dakota_method, method class object) +% variables (structure array, variable class objects) +% responses (structure array, response class objects) +% dparams (structure array, method-independent parameters) +% filei (character, name of .in and .m files) +% +% 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 defines the data to write the dakota .in and +% .m files. it is necessary for multiple reasons. first, +% it collects the parameters and applies some defaults that +% are unique to the environment. second, some analysis package +% variables and/or responses may be treated differently by +% dakota. for example, an analysis package variable may be +% defined as an array, so the QmuSetupDesign brancher will +% create dakota variables for each element of the array. +% finally it calls the functions to write the .in and .m files. +% this function is independent of 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. +% +% "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_in_data(dmeth,variables,responses,dparams,filei,varargin) + +if ~nargin + help dakota_in_data + return +end + +%% parameters + +% get default set of parameters + +params=dakota_in_params(struct()); + +% merge specified parameters into default set, whether or not +% they already exist + +fnames=fieldnames(dparams); + +for i=1:numel(fnames) + if ~isfield(params,fnames{i}) + warning('dakota_in_data:unknown_param',... + 'No parameter ''%s'' in default parameter set.',... + fnames{i}); + end + params.(fnames{i})=dparams.(fnames{i}); +end + +if params.direct && ... + isempty(params.analysis_driver) + params.analysis_driver='matlab'; +end + +if strcmpi(params.analysis_driver,'matlab') && ... + isempty(params.analysis_components) + [pathstr,name,ext,versn] = fileparts(filei); + params.analysis_components=fullfile(pathstr,[name '.m' versn]); +end + +% merge method parameters, though they shouldn't be in dparams + +% dmeth=dmeth_params_merge(dmeth,dparams) + + +%% variables + +fnames=fieldnames(variables); + +for i=1:length(fnames) + +% for linear constraints, just copy + + if isa(variables.(fnames{i}),'linear_inequality_constraint') || ... + isa(variables.(fnames{i}),'linear_equality_constraint' ) + dvar.(fnames{i})=variables.(fnames{i}); + +% for variables, call the setup function + + else + fhandle=str2func([class(variables.(fnames{i})) '.empty']); + dvar.(fnames{i})=fhandle(); + for j=1:length(variables.(fnames{i})) + %call setupdesign + dvar.(fnames{i})=QmuSetupDesign(dvar.(fnames{i}),variables.(fnames{i})(j),params,varargin{:}); + end + end +end + + +%% responses + +fnames=fieldnames(responses); + +for i=1:length(fnames) +% fhandle=str2func([class(responses.(fnames{i})) '.empty']); +% dresp.(fnames{i})=fhandle(); +% for j=1:length(responses.(fnames{i})) +% dresp.(fnames{i})(j)=responses.(fnames{i})(j); +% end + +% currently all response types can just be copied + + dresp.(fnames{i})=responses.(fnames{i}); +end + +%% write files + +%Write in file +dakota_in_write(dmeth,dvar,dresp,params,filei,varargin{:}); + +end Index: /issm/trunk/src/m/dakota/dakota_in_params.m =================================================================== --- /issm/trunk/src/m/dakota/dakota_in_params.m (revision 4129) +++ /issm/trunk/src/m/dakota/dakota_in_params.m (revision 4129) @@ -0,0 +1,210 @@ +% +% populate a Dakota parameter structure. +% +% [params]=dakota_in_params(params) +% +% where the optional input is: +% params (structure array, method-independent parameters) +% +% and the output is the same. +% +% this function takes a structure of method-independent dakota +% parameters, which may be empty, and adds default parameters +% for those parameters which do not exist. +% +% the field names of the structure are identical to the dakota +% parameter names (and are in fact used to write them to the +% files). logical values are used for parameters which have +% no associated data and are determined only by their presence +% or absence. +% +% note that the method-dependent parameters are contained in +% the dakota_method class object. +% +% "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 [params]=dakota_in_params(params) + +if ~nargin + help dakota_in_params + return +end + +%% process the input parameters + +if ~exist('params','var') + params=struct(); +end + +%% strategy section + +if ~isfield(params,'graphics') + params.graphics=false; +end +if ~isfield(params,'tabular_graphics_data') + params.tabular_graphics_data=false; +end +% could use unique file name rather than 'dakota_tabular.dat' +if ~isfield(params,'tabular_graphics_file') + params.tabular_graphics_file=false; +end + +%% method section + +% nearly all method parameters are in the dakota_method class +% or result from the response level lists + +if ~isfield(params,'compute') + params.compute='probabilities'; +end +if ~isfield(params,'distribution') + params.distribution='cumulative'; +end + +%% model section + +%% interface section + +if ~isfield(params,'system') + params.system=false; +end +if ~isfield(params,'fork') + params.fork=false; +end +if ~isfield(params,'direct') + params.direct=false; +end + +% interface parallelism controls + +if ~isfield(params,'asynchronous') + params.asynchronous=true; +end +if ~isfield(params,'evaluation_concurrency') + params.evaluation_concurrency=false; +end +if ~isfield(params,'analysis_concurrency') + params.analysis_concurrency=false; +end +if ~isfield(params,'evaluation_servers') + params.evaluation_servers=false; +end +if ~isfield(params,'evaluation_self_scheduling') + params.evaluation_self_scheduling=false; +end +if ~isfield(params,'evaluation_static_scheduling') + params.evaluation_static_scheduling=false; +end +if ~isfield(params,'analysis_servers') + params.analysis_servers=false; +end +if ~isfield(params,'analysis_self_scheduling') + params.analysis_self_scheduling=false; +end +if ~isfield(params,'analysis_static_scheduling') + params.analysis_static_scheduling=false; +end + +% algebraic mappings + +if ~isfield(params,'algebraic_mappings') + params.algebraic_mappings=false; +end + +% simulation interface controls + +if ~isfield(params,'analysis_driver') + params.analysis_driver=''; +end +if ~isfield(params,'analysis_components') + params.analysis_components=''; +end +if ~isfield(params,'input_filter') + params.input_filter=''; +end +if ~isfield(params,'output_filter') + params.output_filter=''; +end + +if ~isfield(params,'failure_capture') + params.failure_capture='abort'; +end +if ~isfield(params,'deactivate') + params.deactivate='evaluation_cache restart_file'; +end + +% system call or fork interface + +if ~isfield(params,'parameters_file') + params.parameters_file='params.in'; +end +if ~isfield(params,'results_file') + params.results_file='results.out'; +end +if ~isfield(params,'verbatim') + params.verbatim=false; +end +if ~isfield(params,'aprepro') + params.aprepro=false; +end +if ~isfield(params,'file_tag') + params.file_tag=true; +end +if ~isfield(params,'file_save') + params.file_save=true; +end + +% direct function interface + +if ~isfield(params,'processors_per_analysis') + params.processors_per_analysis=false; +end + +%% responses section + +if ~isfield(params,'numerical_gradients') + params.numerical_gradients=false; +end +if ~isfield(params,'method_source') + params.method_source='dakota'; +end +if ~isfield(params,'interval_type') + params.interval_type='forward'; +end +if ~isfield(params,'fd_gradient_step_size') + params.fd_gradient_step_size=0.001; +end +if ~isfield(params,'analytic_gradients') + params.analytic_gradients=false; +end +% mixed_gradients not fully implemented +if ~isfield(params,'mixed_gradients') + params.mixed_gradients=false; +end +if ~isfield(params,'id_analytic_gradients') + params.id_analytic_gradients=false; +end +if ~isfield(params,'id_numerical_gradients') + params.id_numerical_gradients=false; +end +% hessians not fully implemented +if ~isfield(params,'numerical_hessians') + params.numerical_hessians=true; +end +if ~isfield(params,'hessian_gradient_step_size') + params.hessian_gradient_step_size=0.001; +end + +end + Index: /issm/trunk/src/m/dakota/dakota_in_parse.m =================================================================== --- /issm/trunk/src/m/dakota/dakota_in_parse.m (revision 4129) +++ /issm/trunk/src/m/dakota/dakota_in_parse.m (revision 4129) @@ -0,0 +1,803 @@ +% +% read a Dakota .in input file and parse it. +% +% [method,dvar,dresp]=dakota_in_parse(filei) +% +% where the required input is: +% filei (character, name of .in file) +% +% the required output is: +% method (character, dakota method name) +% dvar (structure array, variables) +% dresp (structure array, responses) +% +% the filei will be prompted if empty. the fields of dvar and +% dresp are particular to the data contained within the file. +% +% this function reads a dakota .in input file and parses it +% into the matlab workspace. it operates in a content-driven +% fashion, where it parses whatever input data it encounters +% in the file, rather than searching for data based on the +% particular method. (this makes it independent of method.) +% +% as of now, parameters are generally not parsed. also, the +% variable and response classes are not used for output. +% +% this data would typically be used for modifying and submitting +% a subsequent dakota run. it could also be used with output +% data for post-processing or annotation purposes. +% +% "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 [method,dvar,dresp]=dakota_in_parse(filei) + +if ~nargin + help dakota_in_parse + return +end + +if ~exist('filei' ,'var') || isempty(filei) + filei=input('Input file? ','s'); +end +fidi=fopen(sprintf('%s',filei),'r'); +if (fidi < 0) + error('%s could not be opened.',filei); +end + +%% loop through the file to find the Dakota method + +method=[]; +fseek(fidi,0,'bof'); +[fline]=findline(fidi,'method'); +if ~ischar(fline) + return +end + +[ntokens,tokens]=fltokens(fline); +itoken=1; +[tokens,itoken]=nextkey(fidi,tokens,itoken); +method=tokens{1}{itoken}; +display(sprintf('Dakota method=%s.',method)); + +%% loop through the file to find the Dakota variables + +fseek(fidi,0,'bof'); +[fline]=findline(fidi,'variables'); +if ~ischar(fline) + error('No Dakota variables in file %s.',filei); +end + +[ntokens,tokens]=fltokens(fline); +itoken=1; +[dvar]=variables_parse(fidi,tokens,itoken); + +%% loop through the file to find the Dakota responses + +fseek(fidi,0,'bof'); +[fline]=findline(fidi,'responses'); +if ~ischar(fline) + error('No Dakota responses in file %s.',filei); +end + +[ntokens,tokens]=fltokens(fline); +itoken=1; +[dresp]=responses_parse(fidi,tokens,itoken); + +%% loop through the file to find the Dakota response and probability levels +% (even though they're in method section, process after responses) + +fseek(fidi,0,'bof'); +[fline]=findline(fidi,'method'); + +[ntokens,tokens]=fltokens(fline); +itoken=1; +[dresp]=resplevels(fidi,tokens,itoken,dresp); + +%% loop through the file to verify the end + +display('End of file successfully reached.'); +fclose(fidi); + +end + +%% function to parse the dakota variables + +function [dvar]=variables_parse(fidi,tokens,itoken) + +display('Reading Dakota variables.'); +dvar=[]; +ncdv=0; +nnuv=0; +ncsv=0; + +% read next keyword + +[tokens,itoken]=nextkey(fidi,tokens,itoken); +if ~itoken + warning('variables_parse:empty',... + 'Dakota variables section is empty.'); +end + +% process current keyword +% (note that this is using dakota 4.1 keywords. dakota 4.2 +% keywords are order-dependent.) + +while itoken + keyword=tokens{1}{itoken}; + display(sprintf(' Dakota keyword=%s.',keyword)); + +% switch according to the keyword + + switch lower(keyword) + case 'continuous_design' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + ncdv=tlist; + dvar.cdv=[]; + display(sprintf(' Number of Dakota %s variables=%d.',... + 'continuous_design',ncdv)); + case 'cdv_initial_point' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dvar.cdv(i).initpt =tlist(i); + end + case 'cdv_lower_bounds' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dvar.cdv(i).lower =tlist(i); + end + case 'cdv_upper_bounds' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dvar.cdv(i).upper =tlist(i); + end + case 'cdv_descriptors' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dvar.cdv(i).descriptor=char(tlist(i)); + end + + case 'normal_uncertain' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + nnuv=tlist; + dvar.nuv=[]; + display(sprintf(' Number of Dakota %s variables=%d.',... + 'normal_uncertain',nnuv)); + case 'nuv_means' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dvar.nuv(i).mean =tlist(i); + end + case 'nuv_std_deviations' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dvar.nuv(i).stddev =tlist(i); + end + case 'nuv_lower_bounds' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dvar.nuv(i).lower =tlist(i); + end + case 'nuv_upper_bounds' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dvar.nuv(i).upper =tlist(i); + end + case 'nuv_descriptors' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dvar.nuv(i).descriptor=char(tlist(i)); + end + + case 'continuous_state' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + ncsv=tlist; + dvar.csv=[]; + display(sprintf(' Number of Dakota %s variables=%d.',... + 'continuous_state',ncsv)); + case 'csv_initial_state' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dvar.csv(i).initst =tlist(i); + end + case 'csv_lower_bounds' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dvar.csv(i).lower =tlist(i); + end + case 'csv_upper_bounds' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dvar.csv(i).upper =tlist(i); + end + case 'csv_descriptors' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dvar.csv(i).descriptor=char(tlist(i)); + end + + otherwise + warning('variables_parse:unrec_key',... + 'Unrecognized keyword ''%s''.',keyword); + [tokens,itoken]=nextkey(fidi,tokens,itoken); + end + +% check for eof or start of new section + + if (~itoken) || ... + strncmpi(tokens{1}{itoken},'strategy' ,8) || ... + strncmpi(tokens{1}{itoken},'method' ,6) || ... + strncmpi(tokens{1}{itoken},'model' ,5) || ... + strncmpi(tokens{1}{itoken},'variables',9) || ... + strncmpi(tokens{1}{itoken},'interface',9) || ... + strncmpi(tokens{1}{itoken},'responses',9) + +% supply default descriptors if necessary + + if isfield(dvar,'cdv') && ~isfield(dvar.cdv,'descriptor') + for i=1:ncdv + dvar.cdv(i).descriptor=sprintf('cdv_%d',i); + end + end + if isfield(dvar,'nuv') && ~isfield(dvar.nuv,'descriptor') + for i=1:nnuv + dvar.nuv(i).descriptor=sprintf('nuv_%d',i); + end + end + if isfield(dvar,'csv') && ~isfield(dvar.csv,'descriptor') + for i=1:ncsv + dvar.csv(i).descriptor=sprintf('csv_%d',i); + end + end + return; + end +end + +end + +%% function to parse the dakota responses + +function [dresp]=responses_parse(fidi,tokens,itoken) + +display('Reading Dakota responses.'); +dresp=[]; +nof =0; +nlst=0; +nnic=0; +nnec=0; +nrf =0; + +% read next keyword + +[tokens,itoken]=nextkey(fidi,tokens,itoken); +if ~itoken + warning('responses_parse:empty',... + 'Dakota responses section is empty.'); +end + +% process current keyword + +while itoken + keyword=tokens{1}{itoken}; + display(sprintf(' Dakota keyword=%s.',keyword)); + +% switch according to the keyword + + switch lower(keyword) + case 'num_objective_functions' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + nof =tlist; + dresp.of =[]; + display(sprintf(' Number of Dakota %s=%d.',... + 'objective_functions',nof)); + case 'objective_function_scale_types' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.of(i).scale_type=char(tlist(i)); + end + case 'objective_function_scales' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.of(i).scale =tlist(i); + end + case 'multi_objective_weights' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.of(i).weight =tlist(i); + end + + case 'num_least_squares_terms' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + nlst=tlist; + dresp.lst=[]; + display(sprintf(' Number of Dakota %s=%d.',... + 'least_squares_terms',nlst)); + case 'least_squares_term_scale_types' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.lst(i).scale_type=char(tlist(i)); + end + case 'least_squares_term_scales' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.lst(i).scale =tlist(i); + end + case 'least_squares_weights' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.lst(i).weight =tlist(i); + end + + case 'num_nonlinear_inequality_constraints' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + nnic=tlist; + dresp.nic=[]; + display(sprintf(' Number of Dakota %s=%d.',... + 'nonlinear_inequality_constraints',nnic)); + case 'nonlinear_inequality_scale_types' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nic(i).scale_type=char(tlist(i)); + end + case 'nonlinear_inequality_scales' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nic(i).scale =tlist(i); + end + case 'nonlinear_inequality_lower_bounds' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nic(i).lower =tlist(i); + end + case 'nonlinear_inequality_upper_bounds' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nic(i).upper =tlist(i); + end + + case 'num_nonlinear_equality_constraints' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + nnec=tlist; + dresp.nec=[]; + display(sprintf(' Number of Dakota %s=%d.',... + 'nonlinear_equality_constraints',nnec)); + case 'nonlinear_equality_scale_types' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nec(i).scale_type=char(tlist(i)); + end + case 'nonlinear_equality_scales' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nec(i).scale =tlist(i); + end + case 'nonlinear_equality_targets' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + for i=1:length(tlist) + dresp.nec(i).target =tlist(i); + end + + case 'num_response_functions' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + nrf =tlist; + dresp.rf =[]; + display(sprintf(' Number of Dakota %s=%d.',... + 'response_functions',nrf)); + + case 'response_descriptors' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + desc=tlist; + otherwise + warning('responses_parse:unrec_key',... + 'Unrecognized keyword ''%s''.',keyword); + [tokens,itoken]=nextkey(fidi,tokens,itoken); + end + +% check for eof or start of new section + + if (~itoken) || ... + strncmpi(tokens{1}{itoken},'strategy' ,8) || ... + strncmpi(tokens{1}{itoken},'method' ,6) || ... + strncmpi(tokens{1}{itoken},'model' ,5) || ... + strncmpi(tokens{1}{itoken},'variables',9) || ... + strncmpi(tokens{1}{itoken},'interface',9) || ... + strncmpi(tokens{1}{itoken},'responses',9) + +% assign specified or supply default descriptors + + if exist('desc','var') + idesc=0; + if isfield(dresp,'of' ) + for i=1:nof + idesc=idesc+1; + dresp.of(i).descriptor=char(desc(idesc)); + end + end + if isfield(dresp,'lst') + for i=1:nlst + idesc=idesc+1; + dresp.lst(i).descriptor=char(desc(idesc)); + end + end + if isfield(dresp,'nic') + for i=1:nnic + idesc=idesc+1; + dresp.nic(i).descriptor=char(desc(idesc)); + end + end + if isfield(dresp,'nec') + for i=1:nnec + idesc=idesc+1; + dresp.nec(i).descriptor=char(desc(idesc)); + end + end + if isfield(dresp,'rf' ) + for i=1:nrf + idesc=idesc+1; + dresp.rf(i).descriptor=char(desc(idesc)); + end + end + + else + if isfield(dresp,'of' ) + for i=1:nof + dresp.of(i).descriptor=sprintf('obj_fn_%d',i); + end + end + if isfield(dresp,'lst') + for i=1:nlst + dresp.lst(i).descriptor=sprintf('least_sq_term_%d',i); + end + end + if isfield(dresp,'nic') + for i=1:nnic + dresp.nic(i).descriptor=sprintf('nln_ineq_con_%d',i); + end + end + if isfield(dresp,'nec') + for i=1:nnec + dresp.nec(i).descriptor=sprintf('nln_eq_con_%d',i); + end + end + if isfield(dresp,'rf' ) + for i=1:nrf + dresp.rf(i).descriptor=sprintf('response_fn_%d',i); + end + end + end + return; + end +end + +end + +%% function to read the number and levels of responses + +function [dresp]=resplevels(fidi,tokens,itoken,dresp) + +display('Reading Dakota response levels.'); + +% read next keyword + +[tokens,itoken]=nextkey(fidi,tokens,itoken); +if ~itoken + warning('resplevels:empty',... + 'Dakota method section is empty.'); +end + +% process current keyword + +while itoken + keyword=tokens{1}{itoken}; + display(sprintf(' Dakota keyword=%s.',keyword)); + +% switch according to the keyword + + switch lower(keyword) + case 'nond_sampling' + [tokens,itoken]=nextkey(fidi,tokens,itoken); + case 'nond_local_reliability' + [tokens,itoken]=nextkey(fidi,tokens,itoken); + case 'num_response_levels' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + nresp=tlist; + case 'response_levels' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + nrespl=tlist; + case 'num_probability_levels' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + nprob=tlist; + case 'probability_levels' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + nprobl=tlist; + case 'num_reliability_levels' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + nrel =tlist; + case 'reliability_levels' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + nrell =tlist; + case 'num_gen_reliability_levels' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + ngrel=tlist; + case 'gen_reliability_levels' + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); + ngrell=tlist; + case 'compute' + [tokens,itoken]=nexttoken(fidi,tokens,itoken); + [tokens,itoken]=nextkey(fidi,tokens,itoken); + otherwise + warning('resplevels:unrec_key',... + 'Unrecognized keyword ''%s''.',keyword); + [tokens,itoken]=nextkey(fidi,tokens,itoken); + end + +% check for eof or start of new section + + if (~itoken) || ... + strncmpi(tokens{1}{itoken},'strategy' ,8) || ... + strncmpi(tokens{1}{itoken},'method' ,6) || ... + strncmpi(tokens{1}{itoken},'model' ,5) || ... + strncmpi(tokens{1}{itoken},'variables',9) || ... + strncmpi(tokens{1}{itoken},'interface',9) || ... + strncmpi(tokens{1}{itoken},'responses',9) + +% assemble the lists by response + + if exist('nrespl','var') && isfield(dresp,'rf') + if ~exist('nresp','var') + nresp(1:length(dresp.rf))=floor(length(nrespl)/length(dresp.rf)); + end + ilist=1; + for i=1:length(dresp.rf) + dresp.rf(i).respl=nrespl(ilist:ilist+nresp(i)-1); + ilist=ilist+nresp(i); + end + end + + if exist('nprobl','var') && isfield(dresp,'rf') + if ~exist('nprob','var') + nprob(1:length(dresp.rf))=floor(length(nprobl)/length(dresp.rf)); + end + ilist=1; + for i=1:length(dresp.rf) + dresp.rf(i).probl=nprobl(ilist:ilist+nprob(i)-1); + ilist=ilist+nprob(i); + end + end + + if exist('nrell' ,'var') && isfield(dresp,'rf') + if ~exist('nrel' ,'var') + nrel (1:length(dresp.rf))=floor(length(nrell )/length(dresp.rf)); + end + ilist=1; + for i=1:length(dresp.rf) + dresp.rf(i).rell =nrell (ilist:ilist+nrel (i)-1); + ilist=ilist+nrel (i); + end + end + + if exist('ngrell','var') && isfield(dresp,'rf') + if ~exist('ngrel','var') + ngrel(1:length(dresp.rf))=floor(length(ngrell)/length(dresp.rf)); + end + ilist=1; + for i=1:length(dresp.rf) + dresp.rf(i).grell=ngrell(ilist:ilist+ngrel(i)-1); + ilist=ilist+ngrel(i); + end + end + + return; + end +end + +end + +%% function to find the next keyword + +function [tokens,itoken]=nextkey(fidi,tokens,itoken) + +% start with next token + +[tokens,itoken]=nexttoken(fidi,tokens,itoken); +if ~itoken + return; +end + +% check for equal sign and skip subsequent list + +if (itoken <= length(tokens{1})) && ... + strncmp(tokens{1}{itoken},'=',1) + [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken); +end + +end + +%% function to find the next token + +function [tokens,itoken]=nexttoken(fidi,tokens,itoken) + +% start with next token + +itoken=itoken+1; + +% read next line if necessary + +if (itoken > length(tokens{1})) + fline=readline(fidi); + if isempty(fline) + tokens={}; + itoken=0; + return; + end + [ntokens,tokens]=fltokens(fline); + itoken=1; +end + +end + +%% function to read a list of tokens + +function [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken) + +% start with next token (which should be equal sign, unless +% equal sign was already read to determine existence of list) + +itoken=itoken+1; + +% read next line if necessary + +if (itoken > length(tokens{1})) + fline=readline(fidi); + if isempty(fline) + tokens={}; + itoken=0; + return; + end + [ntokens,tokens]=fltokens(fline); + itoken=1; +end + +% check for equal sign and skip + +if strncmp(tokens{1}{itoken},'=',1) + itoken=itoken+1; +end + +ilist=0; + +% accumulate list until non-numeric and non-quoted-string (or eof) +% is encountered + +while 1 + for i=itoken:length(tokens{1}) + if isnumeric(tokens{1}{i}) + ilist=ilist+1; + tlist(ilist)=tokens{1}{i}; + elseif ischar(tokens{1}{i}) && ... + (strncmp(tokens{1}{i}(1) ,'''',1) && ... + strncmp(tokens{1}{i}(end),'''',1)) || ... + (strncmp(tokens{1}{i}(1) ,'"',1) && ... + strncmp(tokens{1}{i}(end),'"',1)) + ilist=ilist+1; + tlist(ilist)=cellstr(tokens{1}{i}(2:end-1)); + else + itoken=i; + return + end + end + fline=readline(fidi); + if isempty(fline) + tokens={}; + itoken=0; + return; + end + [ntokens,tokens]=fltokens(fline); + itoken=1; +end + +end + +%% function to find a file line starting with a specified string + +function [fline]=findline(fidi,string) + +ipos=ftell(fidi); + +while 1 + fline=readline(fidi); + if isempty(fline) + break; + else + if (strncmpi(fline,string,length(string))) + return; + end + end +end + +% issue warning and reset file position + +warning('findline:str_not_found',... + 'String ''%s'' not found in file.',string); +fseek(fidi,ipos,'bof'); + +end + +%% function to read a file line ignoring comments and blanks + +function [fline]=readline(fidi) + +while 1 + fline=fgetl(fidi); + if ~ischar(fline) + fline=[]; + return; + end + + for ichar=1:length(fline) + if ~strncmp(fline(ichar),' ',1) && ... + ~strncmp(fline(ichar),' ',1) + break; + end + end + if isempty(fline) || ... + (ichar > length(fline)) || ... + strncmp(fline(ichar),'#',1) + continue; + else + return; + end +end + +end + +%% function to parse a file line into tokens + +function [ntokens,tokens]=fltokens(fline) + +if ~ischar(fline) + ntokens=-1; + tokens={}; + return; +end +if isempty(fline) + ntokens=0; + tokens={}; + return; +end + +strings=textscan(fline,'%s','delimiter',' :,'); +%for i=1:length(strings{1}) +% display(sprintf('i=%d; strings{1}{%d}=%s',i,i,strings{1}{i})) +%end +ntokens=0; +tokens{1}{length(strings)}=''; + +for i=1:length(strings{1}) + if isempty(strings{1}{i}) + continue + end + ntokens=ntokens+1; + inum=sscanf(strings{1}{i},'%f'); + if isempty(inum) + tokens{1}{ntokens}=strings{1}{i}; +% display(sprintf('i=%d; tokens{1}{%d}=%s',... +% i,ntokens,tokens{1}{ntokens})) + else + tokens{1}{ntokens}=inum; +% display(sprintf('i=%d; tokens{1}{%d}=%f',... +% i,ntokens,tokens{1}{ntokens})) + end +end + +end Index: /issm/trunk/src/m/dakota/dakota_in_write.m =================================================================== --- /issm/trunk/src/m/dakota/dakota_in_write.m (revision 4129) +++ /issm/trunk/src/m/dakota/dakota_in_write.m (revision 4129) @@ -0,0 +1,321 @@ +% +% write a Dakota .in input file. +% +% []=dakota_in_write(method,dvar,dresp,params,filei,varargin) +% []=dakota_in_write(dmeth ,dvar,dresp,params,filei,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) +% filei (character, name of .in file) +% +% 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. +% +% this function writes a dakota .in input file to be used +% by dakota. this file is independent of 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. (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 []=dakota_in_write(method,dvar,dresp,params,filei,varargin) + +if ~nargin + help dakota_in_write + return +end + +% process the input parameters + +if ~exist('method','var') || isempty(method) + method=input('Method? ','s'); +end +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 + +if ~exist('filei' ,'var') || isempty(filei) + filei=input('Dakota input file to write? ','s'); +end +[pathstr,name,ext,versn] = fileparts(filei); +if isempty(ext) +% fileparts only considers '.in' to be the extension, not '.qmu.in' + ext='.qmu.in'; +end +filei2=fullfile(pathstr,[name ext versn]); + +display(sprintf('Opening Dakota input file ''%s''.',filei2)); +fidi=fopen(sprintf('%s',filei2),'w'); +if (fidi < 0) + error('''%s'' could not be opened.',filei2); +end + +if ~exist('params','var') + params=struct(); +end +params=dakota_in_params(params); + +% write the strategy section + +strategy_write(fidi,params); + +% write the method section + +method_write(fidi,dmeth,dresp,params); + +% write the model section + +model_write(fidi); + +% write the variables section + +variables_write(fidi,dmeth,dvar); + +% write the interface section + +interface_write(fidi,params); + +% write the responses section + +responses_write(fidi,dmeth,dresp,params); + +fclose(fidi); +display('End of file successfully written.'); + +end + +%% function to write the strategy section of the file + +function []=strategy_write(fidi,params) + +display('Writing strategy section of Dakota input file.'); + +fprintf(fidi,'strategy,\n'); +fprintf(fidi,'\tsingle_method\n'); +param_write(fidi,'\t ','graphics','','\n',params); +param_write(fidi,'\t ','tabular_graphics_data','','\n',params); +param_write(fidi,'\t ','tabular_graphics_file',' ''','''\n',params); +fprintf(fidi,'\n'); + +end + +%% function to write the method section of the file + +function []=method_write(fidi,dmeth,dresp,params) + +display('Writing method section of Dakota input file.'); + +fprintf(fidi,'method,\n'); +fprintf(fidi,'\t%s\n',dmeth.method); + +dmeth_params_write(dmeth,fidi); + +% write response levels + +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 + +%% function to write the model section of the file + +function []=model_write(fidi) + +display('Writing model section of Dakota input file.'); + +fprintf(fidi,'model,\n'); +fprintf(fidi,'\t%s\n\n','single'); + +end + +%% function to write the variables section of the file + +function []=variables_write(fidi,dmeth,dvar) + +display('Writing variables section of Dakota input file.'); + +fprintf(fidi,'variables,\n'); + +% variables vary by method + +for i=1:length(dmeth.variables) + fhvar=str2func([dmeth.variables{i} '.dakota_write']); + fhvar(fidi,dvar); +end + +% linear constraints vary by method + +for i=1:length(dmeth.lcspec) + fhvar=str2func([dmeth.lcspec{i} '.dakota_write']); + fhvar(fidi,dvar); +end + +fprintf(fidi,'\n'); + +end + +%% function to write the interface section of the file + +function []=interface_write(fidi,params) + +display('Writing interface section of Dakota input file.'); + +fprintf(fidi,'interface,\n'); + +if ~params.system && ~params.fork && ~params.direct + params.fork=true; +elseif (params.system+params.fork+params.direct > 1) + error('Too many interfaces selected.') +end + +if params.system || params.fork + param_write(fidi,'\t','asynchronous','','\n',params); + param_write(fidi,'\t ','evaluation_concurrency',' = ','\n',params); + param_write(fidi,'\t ','analysis_concurrency',' = ','\n',params); + param_write(fidi,'\t ','evaluation_servers',' = ','\n',params); + param_write(fidi,'\t ','evaluation_self_scheduling','','\n',params); + param_write(fidi,'\t ','evaluation_static_scheduling','','\n',params); + param_write(fidi,'\t ','analysis_servers',' = ','\n',params); + param_write(fidi,'\t ','analysis_self_scheduling','','\n',params); + param_write(fidi,'\t ','analysis_static_scheduling','','\n',params); + param_write(fidi,'\t','algebraic_mappings',' = ','\n',params); + param_write(fidi,'\t','system','','\n',params); + param_write(fidi,'\t','fork','','\n',params); + param_write(fidi,'\t ','analysis_driver',' = ''','''\n',params); + if ~isempty(params.input_filter) + param_write(fidi,'\t ','input_filter',' = ''','''\n',params); + end + if ~isempty(params.output_filter) + param_write(fidi,'\t ','output_filter',' = ''','''\n',params); + end + param_write(fidi,'\t ','failure_capture',' ','\n',params); + param_write(fidi,'\t ','deactivate',' ','\n',params); + param_write(fidi,'\t ','parameters_file',' = ''','''\n',params); + param_write(fidi,'\t ','results_file',' = ''','''\n',params); + param_write(fidi,'\t ','verbatim', '','\n',params); + param_write(fidi,'\t ','aprepro', '','\n',params); + param_write(fidi,'\t ','file_tag', '','\n',params); + param_write(fidi,'\t ','file_save','','\n',params); +elseif params.direct +% Error: asynchronous capability not yet supported in direct interfaces. + param_write(fidi,'\t','algebraic_mappings',' = ','\n',params); + param_write(fidi,'\t','direct','','\n',params); + param_write(fidi,'\t ','analysis_driver',' = ''','''\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); + end + if ~isempty(params.output_filter) + param_write(fidi,'\t ','output_filter',' = ''','''\n',params); + end + param_write(fidi,'\t ','failure_capture',' ','\n',params); + param_write(fidi,'\t ','deactivate',' ','\n',params); + param_write(fidi,'\t ','processors_per_analysis',' = ''','''\n',params); +end + +fprintf(fidi,'\n'); + +end + +%% function to write the responses section of the file + +function []=responses_write(fidi,dmeth,dresp,params) + +display('Writing responses section of Dakota input file.'); + +fprintf(fidi,'responses,\n'); + +% functions, gradients, and hessians vary by method + +rdesc={}; + +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) + fprintf(fidi,'\tresponse_descriptors =\n'); + vector_write(fidi,sprintf('\t '),rdesc,6,76); +end + +ghspec_write(fidi,params,dmeth.ghspec); + +fprintf(fidi,'\n'); + +end + +%% function to write gradient and hessian specifications + +function []=ghspec_write(fidi,params,ghspec) + +% gradients + +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); + param_write(fidi,'\t ','method_source',' ','\n',params); + param_write(fidi,'\t ','interval_type',' ','\n',params); + param_write(fidi,'\t ','fd_gradient_step_size',' = ','\n',params); + elseif params.analytic_gradients + param_write(fidi,'\t','analytic_gradients','','\n',params); +% elseif params.mixed_gradients + end +else + fprintf(fidi,'\tno_gradients\n'); +end + +% hessians (no implemented methods use them yet) + +if find_string(ghspec,'hess') + error('Hessians needed by method but not provided.'); +else + fprintf(fidi,'\tno_hessians\n'); +end + +end Index: /issm/trunk/src/m/dakota/dakota_m_write.m =================================================================== --- /issm/trunk/src/m/dakota/dakota_m_write.m (revision 4129) +++ /issm/trunk/src/m/dakota/dakota_m_write.m (revision 4129) @@ -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/dakota/dakota_moments.m =================================================================== --- /issm/trunk/src/m/dakota/dakota_moments.m (revision 4129) +++ /issm/trunk/src/m/dakota/dakota_moments.m (revision 4129) @@ -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/dakota/dakota_out_parse.m =================================================================== --- /issm/trunk/src/m/dakota/dakota_out_parse.m (revision 4129) +++ /issm/trunk/src/m/dakota/dakota_out_parse.m (revision 4129) @@ -0,0 +1,865 @@ +% +% read a Dakota .out or .dat output file and parse it. +% +% [method,dresp,scm,pcm,srcm,prcm]=dakota_out_parse(filei) +% +% where the required input is: +% filei (character, name of .out file) +% +% the required output is: +% method (character, dakota method name) +% dresp (structure array, responses) +% +% and the optional output is: +% scm (double array, simple correlation matrix) +% pcm (double array, partial correlation matrix) +% srcm (double array, simple rank correlation matrix) +% prcm (double array, partial rank correlation matrix) +% +% the filei will be prompted if empty. the fields of dresp +% are particular to the data contained within the file. the +% scm, pcm, srcm, and prcm are output by dakota only for the +% sampling methods. +% +% this function reads a dakota .out output file and parses it +% into the matlab workspace. it operates in a content-driven +% fashion, where it skips the intermediate data and then parses +% whatever output data it encounters in the order in which it +% exists in the file, rather than searching for data based on +% the particular method. (this makes it independent of method.) +% it also can read and parse the .dat tabular_output file. +% +% this data would typically be used for plotting and other +% post-processing within matlab or excel. +% +% "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 [method,dresp,scm,pcm,srcm,prcm]=dakota_out_parse(filei) + +if ~nargin + help dakota_out_parse + return +end + +if ~exist('filei' ,'var') || isempty(filei) + filei=input('Input file? ','s'); +end +fidi=fopen(sprintf('%s',filei),'r'); +if (fidi < 0) + error('''%s'' could not be opened.',filei); +end + +%% check the first line for the Dakota tabular output file + +method=[]; +fline=fgetl(fidi); +if ~ischar(fline) + error('File ''%s'' is empty.',filei); +end + +if strncmpi(fline,'%eval_id',8) + method='unknown'; + [dresp]=dak_tab_out(fidi,fline); + return +else + fseek(fidi,0,'bof'); +end + +%% loop through the file to find the Dakota method name + +[fline]=findline(fidi,'methodName = '); +if ~ischar(fline) + return +end +% display([' ' deblank(fline)]); + +[ntokens,tokens]=fltokens(fline); +method=tokens{1}{3}; +display(sprintf('Dakota methodName=''%s''.',method)); + +dresp=struct([]); +scm =struct([]); +pcm =struct([]); +srcm=struct([]); +prcm=struct([]); + +%% loop through the file to find the function evaluation summary + +fline=''; +[nfeval]=nfeval_read(fidi,fline); +fline=fgetl(fidi); + +%% process each results section based on content of the file + +while ischar(fline) +% ipos=ftell(fidi); + if isempty(fline) + elseif strncmp(fline,'<<<<< Function evaluation summary',33) + [nfeval]=nfeval_read(fidi,fline); + elseif strncmp(fline,'Statistics based on ',20) + [nsamp]=nsamp_read(fidi,fline); + elseif strncmp(fline,'Moments for each response function',34) + [dresp]=moments_read(fidi,dresp,fline); + elseif strncmp(fline,'95% confidence intervals for each response function',51) + [dresp]=cis_read(fidi,dresp,fline); + elseif strncmp(fline,'Probabilities for each response function',40) + [dresp]=cdfs_read(fidi,dresp,fline); + elseif strncmp(fline,'Simple Correlation Matrix',25) + [scm]=corrmat_read(fidi,'Simple Correlation Matrix',fline); + elseif strncmp(fline,'Partial Correlation Matrix',26) + [pcm]=corrmat_read(fidi,'Partial Correlation Matrix',fline); + elseif strncmp(fline,'Simple Rank Correlation Matrix',30) + [srcm]=corrmat_read(fidi,'Simple Rank Correlation Matrix',fline); + elseif strncmp(fline,'Partial Rank Correlation Matrix',31) + [prcm]=corrmat_read(fidi,'Partial Rank Correlation Matrix',fline); + elseif strncmp(fline,'MV Statistics for ',18) + [dresp]=mvstats_read(fidi,dresp,fline); + elseif strncmp(fline,'<<<<< Best ',11) + [dresp]=best_read(fidi,dresp,fline); + 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) + else + display(['Unexpected line: ' deblank(fline)]); + fline=fgetl(fidi); + end + fline=fgetl(fidi); +% fseek(fidi,ipos,'bof'); +end + +%% loop through the file to verify the end + +% [fline]=findline(fidi,'<<<<< Single Method Strategy completed'); +% if ~ischar(fline) +% return +% end +display('End of file successfully reached.'); +fclose(fidi); + +end + +%% function to parse the dakota tabular output file + +function [dresp]=dak_tab_out(fidi,fline) + +display('Reading Dakota tabular output file.'); + +% process column headings of matrix (skipping eval_id) + +[ntokens,tokens]=fltokens(fline); +desc=cell (1,ntokens-1); +data=zeros(1,ntokens-1); + +for i=1:ntokens-1 + desc(1,i)=cellstr(tokens{1}{i+1}); +end +display(sprintf('Number of columns (Dakota V+R)=%d.',ntokens-1)); + +% process rows of matrix + +nrow=0; +while 1 + fline=fgetl(fidi); + if ~ischar(fline) || isempty(fline) + break; + end + [ntokens,tokens]=fltokens(fline); + +% add row values to matrix (skipping eval_id) + + nrow=nrow+1; + for i=1:ntokens-1 + data(nrow,i)=tokens{1}{i+1}; + end +end +display(sprintf('Number of rows (Dakota func evals)=%d.',nrow)); + +% calculate statistics + +%dmean =mean (data); +%dstddev=std (data,0); +[dmean,dstddev,dmeanci,dstddevci]=... + normfit(data,0.05); + +dmin =min (data); +dquart1=prctile(data,25); +dmedian=median (data); +dquart3=prctile(data,75); +dmax =max (data); + +% same as Dakota scm, Excel correl +dcorrel=corrcoef(data); + +% divide the data into structures for consistency + +for i=1:length(desc) + dresp(i).descriptor=char(desc(i)); + dresp(i).sample =data(:,i); + dresp(i).mean =dmean(i); + dresp(i).stddev =dstddev(i); + dresp(i).meanci =dmeanci(:,i); + dresp(i).stddevci =dstddevci(:,i); + dresp(i).min =dmin(i); + dresp(i).quart1 =dquart1(i); + dresp(i).median =dmedian(i); + dresp(i).quart3 =dquart3(i); + dresp(i).max =dmax(i); +end + +% draw box plot + +% figure +% subplot(2,1,1) +% plot_boxplot(dresp); + +% draw normal probability plot + +% subplot(2,1,2) +% plot_normplot(dresp); + +end + +%% function to find and read the number of function evaluations + +function [nfeval]=nfeval_read(fidi,fline) + +if ~exist('fline','var') || isempty(fline) || ~ischar(fline) + [fline]=findline(fidi,'<<<<< Function evaluation summary'); + if ~ischar(fline) + return + end +end + +[ntokens,tokens]=fltokens(fline); +nfeval=tokens{1}{5}; +display(sprintf(' Dakota function evaluations=%d.',nfeval)); + +end + +%% function to find and read the number of samples + +function [nsamp]=nsamp_read(fidi,fline) + +if ~exist('fline','var') || isempty(fline) || ~ischar(fline) + [fline]=findline(fidi,'Statistics based on '); + if ~ischar(fline) + return + end +end + +[ntokens,tokens]=fltokens(fline); +nsamp=tokens{1}{4}; +display(sprintf(' Dakota samples=%d.',nsamp)); + +end + +%% function to find and read the moments + +function [dresp]=moments_read(fidi,dresp,fline) + +if ~exist('fline','var') || isempty(fline) || ~ischar(fline) + [fline]=findline(fidi,'Moments for each response function'); + if ~ischar(fline) + return + end +end + +display('Reading moments for response functions:'); + +while 1 + fline=fgetl(fidi); + if isempty(fline) + break; + end + [ntokens,tokens]=fltokens(fline); + +% add new response function and moments + + dresp(end+1).descriptor=tokens{1}{ 1}; + display(sprintf(' %s',dresp(end).descriptor)); + dresp(end ).mean =tokens{1}{ 4}; + dresp(end ).stddev =tokens{1}{ 8}; + dresp(end ).coefvar =tokens{1}{13}; +end + +display(sprintf(' Number of Dakota response functions=%d.',... + length(dresp))); + +end + +%% function to find and read the confidence intervals + +function [dresp]=cis_read(fidi,dresp,fline) + +if ~exist('fline','var') || isempty(fline) || ~ischar(fline) + [fline]=findline(fidi,... + '95% confidence intervals for each response function'); + if ~ischar(fline) + return + end +end + +display('Reading 95% confidence intervals for response functions:'); + +while 1 + fline=fgetl(fidi); + if isempty(fline) + break; + end + [ntokens,tokens]=fltokens(fline); + +% find response function associated with confidence intervals + + idresp=0; + for i=1:length(dresp) + if strcmpi(tokens{1}{ 1},dresp(i).descriptor) + idresp=i; + break; + end + end + if ~idresp + idresp=length(dresp)+1; + dresp(idresp).descriptor=tokens{1}{ 1}; + display(sprintf(' %s',dresp(idresp).descriptor)); + end + +% add confidence intervals to response functions + + dresp(i).meanci (1,1)=tokens{1}{ 5}; + dresp(i).meanci (2,1)=tokens{1}{ 6}; + dresp(i).stddevci(1,1)=tokens{1}{12}; + dresp(i).stddevci(2,1)=tokens{1}{13}; +end + +display(sprintf(' Number of Dakota response functions=%d.',... + length(dresp))); + +end + +%% function to find and read the cdf's + +function [dresp]=cdfs_read(fidi,dresp,fline) + +if ~exist('fline','var') || isempty(fline) || ~ischar(fline) + [fline]=findline(fidi,'Probabilities for each response function'); + if ~ischar(fline) + return + end +end + +display('Reading CDF''s for response functions:'); + +while ischar(fline) && ~isempty(fline) + fline=fgetl(fidi); + if ~ischar(fline) + break; + end + +% process header line of cdf + + while ischar(fline) && ~isempty(fline) + [ntokens,tokens]=fltokens(fline); + +% find response function associated with cdf + + idresp=0; + for i=1:length(dresp) + if strcmpi(tokens{1}{ 6},dresp(i).descriptor) + idresp=i; + break; + end + end + if ~idresp + idresp=length(dresp)+1; + dresp(idresp).descriptor=tokens{1}{ 6}; + display(sprintf(' %s',dresp(idresp).descriptor)); + end + +% skip column headings of cdf + + fline=fgetl(fidi); + fline=fgetl(fidi); + +% read and add cdf table to response function + + fline=fgetl(fidi); + icdf=0; + while ischar(fline) && ~isempty(fline) && ... + ~strncmpi(fline,'Cumulative Distribution Function',32) + [ntokens,tokens]=fltokens(fline); + icdf=icdf+1; + dresp(idresp).cdf(icdf,1:4)=NaN; +% in later versions of Dakota, uncalculated columns are now blank + itoken=0; + for i=1:length(fline)/19 + if ~isempty(deblank(fline((i-1)*19+1:i*19))) + itoken=itoken+1; + dresp(idresp).cdf(icdf,i)=tokens{1}{itoken}; + end + end; + fline=fgetl(fidi); + end + end +end + +display(sprintf(' Number of Dakota response functions=%d.',... + length(dresp))); + +end + +%% function to find and read a correlation matrix + +function [cmat]=corrmat_read(fidi,cmstr,fline) + +if ~exist('fline','var') || isempty(fline) || ~ischar(fline) + [fline]=findline(fidi,cmstr); + if ~ischar(fline) + cmat=struct([]); + return + end +end + +display(['Reading ''' fline '''.']); + +cmat.title=fline; + +while ~isempty(fline) + fline=fgetl(fidi); + if ~ischar(fline) + break; + end + +% process column headings of matrix + + [ntokens,tokens]=fltokens(fline); + cmat.column=cell(1,ntokens); + cmat.row =cell(1,1); + cmat.matrix=zeros(1,ntokens); + + for i=1:ntokens + cmat.column(1,i)=cellstr(tokens{1}{i}); + end + +% process rows of matrix, reading until blank line + + nrow=0; + while 1 + fline=fgetl(fidi); + if isempty(fline) + break; + end + [ntokens,tokens]=fltokens(fline); + +% add row heading to matrix + + nrow=nrow+1; + cmat.row (nrow,1)=cellstr(tokens{1}{1}); + +% add row values to matrix + + for i=2:ntokens + cmat.matrix(nrow,i-1)=tokens{1}{i}; + end + end +end + +end + +%% function to find and read the MV statistics + +function [dresp]=mvstats_read(fidi,dresp,fline) + +if ~exist('fline','var') || isempty(fline) || ~ischar(fline) + [fline]=findline(fidi,'MV Statistics for '); + if ~ischar(fline) + return + end +end + +display('Reading MV statistics for response functions:'); +ndresp=0; + +while ischar(fline) && ~isempty(fline) && ... + strncmpi(fline,'MV Statistics for ',18) + +% add new response function and moments + + [ntokens,tokens]=fltokens(fline); + dresp(end+1).descriptor=tokens{1}{4}; + display(sprintf(' %s',dresp(end).descriptor)); + fline=fgetl(fidi); + [ntokens,tokens]=fltokens(fline); + dresp(end ).mean =tokens{1}{5}; + fline=fgetl(fidi); + [ntokens,tokens]=fltokens(fline); + dresp(end ).stddev =tokens{1}{7}; + +% read and add importance factors to response function + + idvar=0; + fline=fgetl(fidi); + if ~ischar(fline) + break; + end + + while ischar(fline) && ~isempty(fline) && ... + strncmpi(fline,' Importance Factor for variable ',33) + [ntokens,tokens]=fltokens(fline); + idvar=idvar+1; + dresp(end).var (idvar,1)=cellstr(tokens{1}{5}); + dresp(end).impfac(idvar,1)= tokens{1}{7}; + + fline=fgetl(fidi); + end + +% if importance factors missing, skip to cdf + + if ~idvar + display(' Importance Factors not available.'); + dresp(end).var ={}; + dresp(end).impfac=[]; + while ischar(fline) && ... + ~strncmpi(fline,'Cumulative Distribution Function',32) && ... + ~strncmpi(fline,'MV Statistics for ',18) && ... + ~strncmp (fline,'-',1) + fline=fgetl(fidi); + end + end + +% process header line of cdf + + icdf=0; + + while ischar(fline) && ~isempty(fline) && ... + strncmpi(fline,'Cumulative Distribution Function',32) + [ntokens,tokens]=fltokens(fline); + +% find response function associated with cdf + + idresp=0; + for i=1:length(dresp) + if strcmpi(tokens{1}{ 6},dresp(i).descriptor) + idresp=i; + break; + end + end + if ~idresp + idresp=length(dresp)+1; + dresp(idresp).descriptor=tokens{1}{ 6}; + display(sprintf(' %s',dresp(idresp).descriptor)); + end + +% skip column headings of cdf + + fline=fgetl(fidi); + fline=fgetl(fidi); + +% read and add cdf table to response function + + fline=fgetl(fidi); + while ~isempty(fline) && ... + ~strncmpi(fline,'MV Statistics for ',18) && ... + ~strncmp (fline,'-',1) + [ntokens,tokens]=fltokens(fline); + icdf=icdf+1; + dresp(idresp).cdf(icdf,1)=tokens{1}{1}; + dresp(idresp).cdf(icdf,2)=tokens{1}{2}; + if (ntokens == 4) + dresp(idresp).cdf(icdf,3)=tokens{1}{3}; + dresp(idresp).cdf(icdf,4)=tokens{1}{4}; + else + dresp(idresp).cdf(icdf,3)=NaN; + dresp(idresp).cdf(icdf,4)=NaN; + end + fline=fgetl(fidi); + end + end + +% if cdf missing, skip to end of response function + + if ~icdf + display(' Cumulative Distribution Function not available.'); + dresp(ndresp).cdf=[]; + while ischar(fline) && ... + ~strncmpi(fline,'MV Statistics for ',18) && ... + ~strncmp (fline,'-',1) + fline=fgetl(fidi); + end + end + +end + +display(sprintf(' Number of Dakota response functions=%d.',... + length(dresp))); + +end + +%% function to find and read the best evaluation + +function [dresp]=best_read(fidi,dresp,fline) + +if ~exist('fline','var') || isempty(fline) || ~ischar(fline) + [fline]=findline(fidi,'<<<<< Best '); + if ~ischar(fline) + return + end +end + +if isempty(dresp) + dresp(end+1).best=[]; +end +display('Reading values for best function evaluation:'); + +while ischar(fline) && ~isempty(fline) && ... + strncmpi(fline,'<<<<< Best ',11) + [ntokens,tokens]=fltokens(fline); + +% read and add best parameter(s) + + if strncmpi(cellstr(tokens{1}{3}),'parameter', 9) + display([' ' deblank(fline)]); + + fline=fgetl(fidi); + dresp.best.param =[]; + dresp.best.descriptor={}; + + while ischar(fline) && ~isempty(fline) && ... + ~strncmpi(fline,'<<<<< Best ',11) + [ntokens,tokens]=fltokens(fline); + dresp.best.param (end+1,1)= tokens{1}{1}; + dresp.best.descriptor(end+1,1)=cellstr(tokens{1}{2}); + fline=fgetl(fidi); + end + +% read and add best objective function(s) + + elseif strncmpi(cellstr(tokens{1}{3}),'objective', 9) && ... + strncmpi(cellstr(tokens{1}{4}),'function' , 8) + display([' ' deblank(fline)]); + + fline=fgetl(fidi); + dresp.best.of=[]; + + while ischar(fline) && ~isempty(fline) && ... + ~strncmpi(fline,'<<<<< Best ',11) + [ntokens,tokens]=fltokens(fline); + dresp.best.of(end+1,1)= tokens{1}{1}; + fline=fgetl(fidi); + 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) + + elseif strncmpi(cellstr(tokens{1}{3}),'residual', 8) && ... + strncmpi(cellstr(tokens{1}{4}),'term' , 4) + display([' ' deblank(fline)]); + + fline=fgetl(fidi); + dresp.best.res=[]; + + while ischar(fline) && ~isempty(fline) && ... + ~strncmpi(fline,'<<<<< Best ',11) + [ntokens,tokens]=fltokens(fline); + dresp.best.res(end+1,1)= tokens{1}{1}; + fline=fgetl(fidi); + end + +% read and add best constraint value(s) + + elseif strncmpi(cellstr(tokens{1}{3}),'constraint',10) && ... + strncmpi(cellstr(tokens{1}{4}),'value' , 5) + display([' ' deblank(fline)]); + + fline=fgetl(fidi); + dresp.best.nc=[]; + + while ischar(fline) && ~isempty(fline) && ... + ~strncmpi(fline,'<<<<< Best ',11) + [ntokens,tokens]=fltokens(fline); + dresp.best.nc(end+1,1)= tokens{1}{1}; + fline=fgetl(fidi); + end + +% read and add best data captured + + elseif strncmpi(cellstr(tokens{1}{3}),'data' , 4) && ... + strncmpi(cellstr(tokens{1}{4}),'captured', 8) + display([' ' deblank(fline)]); + dresp.best.eval= tokens{1}{8}; + + fline=fgetl(fidi); + + 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) ' (ignored)']); + + fline=fgetl(fidi); + + while ischar(fline) && ~isempty(fline) && ... + ~strncmpi(fline,'<<<<< Best ',11) + fline=fgetl(fidi); + end + end +end + +end + +%% function to find and read the volumetric uniformity measures + +function [dresp]=vum_read(fidi,dresp,fline) + +if ~exist('fline','var') || isempty(fline) || ~ischar(fline) + [fline]=findline(fidi,'The following lists volumetric uniformity measures'); + if ~ischar(fline) + return + end +end + +if isempty(dresp) + dresp(end+1).vum=[]; +end +display('Reading measures for volumetric uniformity.'); + +fline=fgetl(fidi); +fline=fgetl(fidi); + +while ischar(fline) && ~isempty(fline) + [ntokens,tokens]=fltokens(fline); + switch lower(tokens{1}{1}) + case 'chi' + dresp.vum.chi=tokens{1}{4}; + case 'd' + dresp.vum.d =tokens{1}{4}; + case 'h' + dresp.vum.h =tokens{1}{4}; + case 'tau' + dresp.vum.tau=tokens{1}{4}; + end + fline=fgetl(fidi); +end + +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 + +function [fline]=findline(fidi,string) + +ipos=ftell(fidi); + +while 1 + fline=fgetl(fidi); + if ~ischar(fline) + break; + else + if (strncmpi(fline,string,length(string))) + return; + end + end +end + +% issue warning and reset file position + +warning('findline:str_not_found',... + 'String ''%s'' not found in file.',string); +fseek(fidi,ipos,'bof'); + +end + +%% function to parse a file line into tokens + +function [ntokens,tokens]=fltokens(fline) + +if ~ischar(fline) + ntokens=-1; + tokens={}; + return; +end +if isempty(fline) + ntokens=0; + tokens={}; + return; +end + +strings=textscan(fline,'%s','delimiter',' :'); +%for i=1:length(strings{1}) +% display(sprintf('i=%d; strings{1}{%d}=%s',i,i,strings{1}{i})) +%end +ntokens=0; +tokens{1}{length(strings)}=''; + +for i=1:length(strings{1}) + if isempty(strings{1}{i}) + continue + end + ntokens=ntokens+1; + inum=sscanf(strings{1}{i},'%f'); + if isempty(inum) + tokens{1}{ntokens}=strings{1}{i}; +% display(sprintf('i=%d; tokens{1}{%d}=%s',... +% i,ntokens,tokens{1}{ntokens})) + else + tokens{1}{ntokens}=inum; +% display(sprintf('i=%d; tokens{1}{%d}=%f',... +% i,ntokens,tokens{1}{ntokens})) + end +end + +end Index: /issm/trunk/src/m/dakota/dakota_resp_uconv.m =================================================================== --- /issm/trunk/src/m/dakota/dakota_resp_uconv.m (revision 4129) +++ /issm/trunk/src/m/dakota/dakota_resp_uconv.m (revision 4129) @@ -0,0 +1,112 @@ +% +% convert the units for dakota responses. +% +% [dresp]=dakota_resp_uconv(dresp) +% +% where the required input is: +% dresp (structure array, responses) +% +% the required output is: +% dresp (structure array, responses) +% +% this function reads through a dakota response structure, and +% for those quantities whose descriptors are recognized, converts +% the units of all the applicable fields. a "unit" field is added +% to the response structure. +% +% this data would typically be read by dakota_out_parse and be used +% for plotting and other post-processing within matlab or excel. +% +% "Copyright 2010, 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 [dresp]=dakota_resp_uconv(dresp) + +if ~nargin + help dakota_resp_uconv + return +end + +if ~isstruct(dresp) + error('''%s'' is not a structure array.',inputname(1)); +end +if ~isfield(dresp,'descriptor') + error('''%s'' does not have a descriptor field.',inputname(1)); +end + +%% define the conversion factors + +sec_per_yr=365.2425*24*60*60; % mean gregorian year +m_per_km=1000; +kg_per_gton=10^12; + +%% loop through the response array + +for i=1:numel(dresp) + dresp(i).unit=''; + if ~isempty(strfind(dresp(i).descriptor,'vel')) % in m/sec + dresp(i)=drespi_conv(dresp(i),sec_per_yr,'m/yr'); + elseif ~isempty(strfind(dresp(i).descriptor,'misfit')) % in m^2*(m/sec)^2 + dresp(i)=drespi_conv(dresp(i),1/m_per_km^2*sec_per_yr^2,'km^2*(m/yr)^2'); + elseif ~isempty(strfind(dresp(i).descriptor,'mass_flux')) % in kg/sec + dresp(i)=drespi_conv(dresp(i),1/kg_per_gton*sec_per_yr,'Gton/yr'); + else + disp(['Skipping response ''' dresp(i).descriptor '''.']); + end +end + +end + +%% function to convert the units of a dakota response + +function [dresp]=drespi_conv(dresp,unew_per_uold,ulab) + +disp(['Converting response ''' dresp.descriptor ''' to ' ulab '.']); + +% loop over the fields, converting only the appropriate ones + +fnames=fieldnames(dresp); +for i=1:length(fnames) + switch fnames{i} + case {'sample',... + 'mean',... + 'stddev',... + 'meanci',... + 'stddevci',... + 'min',... + 'quart1',... + 'median',... + 'quart3',... + 'max'} % appropriate to convert + dresp.(fnames{i})=dresp.(fnames{i})*unew_per_uold; + case {'cdf'} % only responses, not probs or reliabilities + dresp.cdf(:,1)=dresp.cdf(:,1)*unew_per_uold; + case {'descriptor',... + 'coefvar',... + 'var',... + 'impfac'} % unitless (or non-numeric) + continue; + case {'best',... + 'vum',... + 'unit'} % other + continue; + otherwise + disp(['Unrecognized field ''' fnames{i} '''.']); + end +end + +if isfield(dresp,'unit') + dresp.unit=ulab; +end + +end Index: /issm/trunk/src/m/dakota/examples/dakota_rosenbrock_2d.m =================================================================== --- /issm/trunk/src/m/dakota/examples/dakota_rosenbrock_2d.m (revision 4129) +++ /issm/trunk/src/m/dakota/examples/dakota_rosenbrock_2d.m (revision 4129) @@ -0,0 +1,60 @@ +% +% multi-dimensional parameter study for rosenbrock case +% (see Users4.2.pdf, Sec. 2.4.1) +% +% "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 [dout,ddat]=dakota_rosenbrock_2d() + +% define dakota variables as continuous design +% (may use set 1 or set 2 of variables, but not both) +dvar(1).cdv(1)=continuous_design('x1',0,-2,2); +dvar(1).cdv(2)=continuous_design('x2',0,-2,2); +dvar(2).x1=continuous_design('',0,-2,2); +dvar(2).x2=continuous_design('',0,-2,2); + +% define dakota response as objective function +% (may use set 1 or set 2 of responses, but not both) +dresp(1).of=objective_function('f'); +dresp(2).f=objective_function(''); + +% define dakota method and specify method-dependent parameters +dmeth=dakota_method('multidim'); +dmeth=dmeth_params_set(dmeth,'partitions',[8 8]); + +% specify method-independent parameters +% (dakota_in_params does not need to be called, but provides a template) +dparams=dakota_in_params([]); +dparams.direct=true; +dparams.analysis_driver='rosenbrock'; +dparams.tabular_graphics_data=true; +dparams.tabular_graphics_file='dakota_rosenbrock_2d.dat'; + +% write out dakota input file +dakota_in_write(dmeth,dvar(2),dresp(2),dparams,'dakota_rosenbrock_2d.in') + +% execute dakota +!dakota -i dakota_rosenbrock_2d.in -o dakota_rosenbrock_2d.out + +% read dakota output and tabular data files +% (output file for parameter studies has no interesting info) +[method,dout]=dakota_out_parse('dakota_rosenbrock_2d.out'); +[~ ,ddat]=dakota_out_parse('dakota_rosenbrock_2d.dat'); + +% perform any desired plotting +plot_rvsv_surf(ddat,{'x1','x2'},ddat,{'f'}) + +end + Index: /issm/trunk/src/m/dakota/examples/dakota_rosenbrock_ls.m =================================================================== --- /issm/trunk/src/m/dakota/examples/dakota_rosenbrock_ls.m (revision 4129) +++ /issm/trunk/src/m/dakota/examples/dakota_rosenbrock_ls.m (revision 4129) @@ -0,0 +1,58 @@ +% +% least squares study for rosenbrock case +% (see Users4.2.pdf, Sec. 2.4.1) +% +% "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 [dout,ddat]=dakota_rosenbrock_ls() + +% define dakota variables as continuous design +dvar.x1=continuous_design('',-1.2,-2,2); +dvar.x2=continuous_design('', 1.0,-2,2); + +% define dakota response as least-squares terms +dresp.f1sq=least_squares_term(''); +dresp.f2sq=least_squares_term(''); + +% define dakota method and specify method-dependent parameters +dmeth=dakota_method('nl2sol'); +dmeth=dmeth_params_set(dmeth,'max_iterations',100,... + 'convergence_tolerance',1.e-4); + +% specify method-independent parameters +% (dakota_in_params does not need to be called, but provides a template) +dparams=dakota_in_params([]); +dparams.direct=true; +dparams.analysis_driver='rosenbrock'; +dparams.tabular_graphics_data=true; +dparams.tabular_graphics_file='dakota_rosenbrock_ls.dat'; +dparams.analytic_gradients=true; + +% write out dakota input file +dakota_in_write(dmeth,dvar,dresp,dparams,'dakota_rosenbrock_ls.in') + +% execute dakota +!dakota -i dakota_rosenbrock_ls.in -o dakota_rosenbrock_ls.out + +% read dakota output and tabular data files +% (output file for parameter studies has no interesting info) +[method,dout]=dakota_out_parse('dakota_rosenbrock_ls.out'); +[~ ,ddat]=dakota_out_parse('dakota_rosenbrock_ls.dat'); + +% perform any desired plotting +% plot_rvsv_surf(ddat,{'x1','x2'},ddat,{'f'}) + +end + Index: /issm/trunk/src/m/dakota/examples/dakota_rosenbrock_nond.m =================================================================== --- /issm/trunk/src/m/dakota/examples/dakota_rosenbrock_nond.m (revision 4129) +++ /issm/trunk/src/m/dakota/examples/dakota_rosenbrock_nond.m (revision 4129) @@ -0,0 +1,70 @@ +% +% monte carlo sampling for rosenbrock case +% (see Users4.2.pdf, Sec. 2.4.9) +% +% "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 [dout,ddat,scm,pcm,srcm,prcm]=dakota_rosenbrock_nond() + +% define dakota variables as uniform uncertain +% (may use set 1 or set 2 of variables, but not both) +dvar(1).uuv(1)=uniform_uncertain('x1',-2,2); +dvar(1).uuv(2)=uniform_uncertain('x2',-2,2); +dvar(2).x1=uniform_uncertain('',-2,2); +dvar(2).x2=uniform_uncertain('',-2,2); + +% define dakota response as response function +% (may use set 1 or set 2 of responses, but not both) +dresp(1).rf=response_function('f',[100]); +dresp(2).f=response_function('',[100]); + +% define dakota method and specify method-dependent parameters +dmeth=dakota_method('nond_samp'); +dmeth=dmeth_params_set(dmeth,'samples',200,... + 'seed',17,... + 'sample_type','random'); + +% specify method-independent parameters +% (dakota_in_params does not need to be called, but provides a template) +dparams=dakota_in_params([]); +dparams.direct=true; +dparams.analysis_driver='rosenbrock'; +dparams.tabular_graphics_data=true; +dparams.tabular_graphics_file='dakota_rosenbrock_nond.dat'; + +% write out dakota input file +dakota_in_write(dmeth,dvar(2),dresp(2),dparams,'dakota_rosenbrock_nond.in') + +% execute dakota +!dakota -i dakota_rosenbrock_nond.in -o dakota_rosenbrock_nond.out + +% read dakota output and tabular data files +[method,dout,scm,pcm,srcm,prcm]=dakota_out_parse('dakota_rosenbrock_nond.out'); +[~ ,ddat ]=dakota_out_parse('dakota_rosenbrock_nond.dat'); + +% perform any desired plotting +plot_boxplot(ddat,{'f'}) +plot_normplot(ddat,{'f'}) +plot_sampdist_bars(ddat,{'f'}) +plot_normdist_bars(ddat,{'f'}) +plot_hist_norm(ddat,{'f'}) +plot_hist_norm(ddat,{'f'},'hmin',0,'hmax',200) +plot_hist_norm_ci(ddat,{'f'}) +plot_hist_norm_ci(ddat,{'f'},'ciplt','line','cdfplt','off','ymin1',0,'ymax1',0.15) +plot_prob_bars(dout,{'f'}) +plot_rlev_bars_ci(ddat,{'f'},'xtlrot',90) + +end + Index: /issm/trunk/src/m/dakota/examples/dakota_rosenbrock_vector.m =================================================================== --- /issm/trunk/src/m/dakota/examples/dakota_rosenbrock_vector.m (revision 4129) +++ /issm/trunk/src/m/dakota/examples/dakota_rosenbrock_vector.m (revision 4129) @@ -0,0 +1,64 @@ +% +% vector parameter study for rosenbrock case +% (see Users4.2.pdf, Sec. 2.4.2) +% +% "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 [dout,ddat]=dakota_rosenbrock_vector() + +% define dakota variables as continuous design +% (may use set 1 or set 2 of variables, but not both) +dvar(1).cdv(1)=continuous_design('x1',-0.3); +dvar(1).cdv(2)=continuous_design('x2', 0.2); +dvar(2).x1=continuous_design('',-0.3); +dvar(2).x2=continuous_design('', 0.2); + +% define dakota response as objective function +% (may use set 1 or set 2 of responses, but not both) +dresp(1).of=objective_function('f'); +dresp(2).f=objective_function(''); + +% define dakota method and specify method-dependent parameters +dmeth=dakota_method('vector'); +% (ref. manual says step_length is required, but user's doesn't) +dmeth=dmeth_params_set(dmeth,'final_point',[1.1 1.3],... + 'num_steps',10); + +% specify method-independent parameters +% (dakota_in_params does not need to be called, but provides a template) +dparams=dakota_in_params([]); +dparams.direct=true; +dparams.analysis_driver='rosenbrock'; +dparams.tabular_graphics_data=true; +dparams.tabular_graphics_file='dakota_rosenbrock_vector.dat'; + +% write out dakota input file +dakota_in_write(dmeth,dvar(2),dresp(2),dparams,'dakota_rosenbrock_vector.in') + +% execute dakota +!dakota -i dakota_rosenbrock_vector.in -o dakota_rosenbrock_vector.out + +% read dakota output and tabular data files +% (output file for parameter studies has no interesting info) +[method,dout]=dakota_out_parse('dakota_rosenbrock_vector.out'); +[~ ,ddat]=dakota_out_parse('dakota_rosenbrock_vector.dat'); +clear dum + +% perform any desired plotting +plot_rvsv_line(ddat,{'x1','x2'},ddat,{'f'}) +plot_rvsv_line(ddat,{'x1','x2'},ddat,{'f'},'nplotr',1,'nplotc',2) + +end + Index: /issm/trunk/src/m/dakota/plot_boxplot.m =================================================================== --- /issm/trunk/src/m/dakota/plot_boxplot.m (revision 4129) +++ /issm/trunk/src/m/dakota/plot_boxplot.m (revision 4129) @@ -0,0 +1,145 @@ +% +% plot a box plot of the responses. +% +% []=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: +% 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: +% 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 +% the same length. +% +% this data would typically be contained in the dakota tabular +% output file 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_boxplot(varargin) + +if ~nargin + help plot_boxplot + return +end + +%% process input data and assemble into matrices 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 + + 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 + +figure +boxplot(sampr,'labels',descr,'notch','on') +ax1=gca; + +% add the annotation + +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/dakota/plot_hist_norm.m =================================================================== --- /issm/trunk/src/m/dakota/plot_hist_norm.m (revision 4129) +++ /issm/trunk/src/m/dakota/plot_hist_norm.m (revision 4129) @@ -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/dakota/plot_hist_norm_ci.m =================================================================== --- /issm/trunk/src/m/dakota/plot_hist_norm_ci.m (revision 4129) +++ /issm/trunk/src/m/dakota/plot_hist_norm_ci.m (revision 4129) @@ -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/dakota/plot_if_bars.m =================================================================== --- /issm/trunk/src/m/dakota/plot_if_bars.m (revision 4129) +++ /issm/trunk/src/m/dakota/plot_if_bars.m (revision 4129) @@ -0,0 +1,223 @@ +% +% plot a stacked bar chart of the importance factors. +% +% []=plot_if_bars(dresp ,params) +% []=plot_if_bars(dresp,descr,params) +% +% where the required input is: +% dresp (structure array, responses) +% 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) +% xtlrot (numeric, rotation in degrees of x-tick 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 importance factors, and annotates it with the +% description. the legend labels are constructed from the +% variable list. +% +% this data would typically be contained in the dakota output +% file from a local reliability 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_if_bars(varargin) + +if ~nargin + help plot_if_bars + return +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; +end + +if ~exist('isort','var') || isempty(isort) + isort=0; +end + +%% sort the data, if necessary + +if (isort) + ifmean=mean(ifr,1); + if (isort > 0) + [ifmean,index]=sort(ifmean,'descend'); + else + [ifmean,index]=sort(ifmean,'ascend' ); + end + clear ifmean + + dvar=dvar(index); + ifr =ifr (:,index); +end + +%% filter the data, if necessary + +if (ifmin > 0) + nif=length(dvar); + dvar(nif+1,1)=cellstr(sprintf('others < %f',ifmin)); + ifr (:,nif+1)=0.; + + nif2=0; + dvar2=cell (size(dvar)); + ifr2 =zeros(size(ifr )); + +% sum filtered rows and copy unfiltered rows + + for i=1:nif + if (max(ifr(:,i)) < ifmin) + ifr(:,nif+1)=ifr(:,nif+1)+ifr(:,i); + else + nif2=nif2+1; + dvar2(nif2) =dvar(i); + ifr2 (:,nif2)=ifr (:,i); + end + end + +% copy sums + + dvar2(nif2+1) =dvar(nif+1); + ifr2 (:,nif2+1)=ifr (:,nif+1); + +% copy back and truncate filtered rows + + clear dvar ifr + if (isort >= 0) + dvar(1:nif2+1) =dvar2(1:nif2+1); + ifr (:,1:nif2+1)=ifr2 (:,1:nif2+1); + else + dvar(1 ) =dvar2( nif2+1); + dvar(2:nif2+1) =dvar2(1:nif2 ); + ifr (:,1 )=ifr2 (:, nif2+1); + ifr (:,2:nif2+1)=ifr2 (:,1:nif2 ); + end + clear nif nif2 dvar2 ifr2 +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(dresp) == 1 + ifr=[ifr; ifr]; +end + +figure +hl1=bar(ifr,'stacked'); + +ax1=gca; +if length(dresp) == 1 + set(ax1,'xlim',[0.5 1.5]) +end + +% 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 + +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') + +hleg1=legend(ax1,dvar,'Location','EastOutside',... + 'Orientation','vertical','Interpreter','none'); + +end Index: /issm/trunk/src/m/dakota/plot_normdist_bars.m =================================================================== --- /issm/trunk/src/m/dakota/plot_normdist_bars.m (revision 4129) +++ /issm/trunk/src/m/dakota/plot_normdist_bars.m (revision 4129) @@ -0,0 +1,207 @@ +% +% plot a stacked bar chart of the sample distributions. +% +% []=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) +% 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) +% +% and the optional fields of dresp 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: +% 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 +% a stacked bar plot of the list of samples, where the bars +% are stacked by the given or default probability levels +% calculated from a normal distribution, and annotates it with +% the description. the mean and standard deviation will be +% calculated from the samples if they do not already exist. +% the legend labels can be given or constructed from the +% probability levels. +% +% this data would typically be contained in the dakota tabular +% output file 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_normdist_bars(varargin) + +if ~nargin + help plot_normdist_bars + 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=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) + [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) + 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 + +%% 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(dresp) == 1 + data=[data; data]; +end + +figure +hl1=bar(data,'stacked'); +% set barseries properties for lowest value +whitebg('white') +set(hl1(1),'FaceColor','white') +set(hl1(1),'Visible','off') + +ax1=gca; +if length(dresp) == 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('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(prob)); + for i=1:length(prob) + lstr(i)=cellstr(sprintf('%g%%',100*prob(i))); + end +end + +hleg1=legend(ax1,lstr,'Location','EastOutside',... + 'Orientation','vertical','Interpreter','none'); + +end Index: /issm/trunk/src/m/dakota/plot_normplot.m =================================================================== --- /issm/trunk/src/m/dakota/plot_normplot.m (revision 4129) +++ /issm/trunk/src/m/dakota/plot_normplot.m (revision 4129) @@ -0,0 +1,148 @@ +% +% plot a normal probability plot of the responses. +% +% []=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 +% a matlab normal probability plot of the list of samples +% and annotates it with the description. the lists of samples +% need not all be the same length. +% +% this data would typically be contained in the dakota tabular +% output file 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_normplot(varargin) + +if ~nargin + help plot_normplot + return +end + +%% process input data and assemble into matrices 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 + + 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 + +figure +normplot(sampr) +ax1=gca; + +% add the annotation + +xlim('auto') +[xlims]=xlim; +if exist('xmin','var') + xlims(1)=xmin; +end +if exist('xmax','var') + xlims(2)=xmax; +end +xlim(xlims) + +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'); + +end Index: /issm/trunk/src/m/dakota/plot_prob_bars.m =================================================================== --- /issm/trunk/src/m/dakota/plot_prob_bars.m (revision 4129) +++ /issm/trunk/src/m/dakota/plot_prob_bars.m (revision 4129) @@ -0,0 +1,191 @@ +% +% plot a stacked bar chart of the probabilities in the CDF. +% +% []=plot_prob_bars(dresp ,params) +% []=plot_prob_bars(dresp,descr,params) +% +% where the required input is: +% dresp (structure array, responses) +% 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 +% a stacked bar plot of the responses, where the bars are +% stacked by the probabilities corresponding to the given +% response levels in the CDF, and annotates it with the +% description. the legend labels can be given or constructed +% from the response levels. +% +% this data would typically be contained in the dakota output +% file 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_prob_bars(varargin) + +if ~nargin + help plot_prob_bars + return +end + +%% assemble the data into a matrix and calculate the increments + + +%% 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 + +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 + +% if there's only one row, Matlab 7.5 interprets it as a column, +% so add an extra row, then reduce xlim + +if length(dresp) == 1 + cdfr=[cdfr; cdfr]; +end + +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; +if length(dresp) == 1 + set(ax1,'xlim',[0.5 1.5]) +end + +% 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 + +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(lcdfr)); + for i=1:max(lcdfr) + lstr(i)=cellstr(sprintf('%g',... + dresp(find(lcdfr == max(lcdfr),1,'first')).cdf(i,1))); + end + if ~isempty(find(lcdfr < max(lcdfr))) + warning('Variable number of levels for responses.'); + end +end + +hleg1=legend(ax1,lstr,'Location','EastOutside',... + 'Orientation','vertical','Interpreter','none'); + +end Index: /issm/trunk/src/m/dakota/plot_rlev_bars.m =================================================================== --- /issm/trunk/src/m/dakota/plot_rlev_bars.m (revision 4129) +++ /issm/trunk/src/m/dakota/plot_rlev_bars.m (revision 4129) @@ -0,0 +1,172 @@ +% +% plot a stacked bar chart of the response levels in the cdf. +% +% []=plot_rlev_bars(dresp ,params) +% []=plot_rlev_bars(dresp,descr,params) +% +% where the required input is: +% dresp (structure array, responses) +% 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 +% 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 legend labels can be given or constructed +% from the probabilities. +% +% this data would typically be contained in the dakota output +% file 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_rlev_bars(varargin) + +if ~nargin + help plot_rlev_bars + return +end + +%% 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,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 + +%% 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(dresp) == 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(dresp) == 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)) + 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/dakota/plot_rlev_bars_ci.m =================================================================== --- /issm/trunk/src/m/dakota/plot_rlev_bars_ci.m (revision 4129) +++ /issm/trunk/src/m/dakota/plot_rlev_bars_ci.m (revision 4129) @@ -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/dakota/plot_rvsv_line.m =================================================================== --- /issm/trunk/src/m/dakota/plot_rvsv_line.m (revision 4129) +++ /issm/trunk/src/m/dakota/plot_rvsv_line.m (revision 4129) @@ -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/dakota/plot_rvsv_scat3.m =================================================================== --- /issm/trunk/src/m/dakota/plot_rvsv_scat3.m (revision 4129) +++ /issm/trunk/src/m/dakota/plot_rvsv_scat3.m (revision 4129) @@ -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/dakota/plot_rvsv_surf.m =================================================================== --- /issm/trunk/src/m/dakota/plot_rvsv_surf.m (revision 4129) +++ /issm/trunk/src/m/dakota/plot_rvsv_surf.m (revision 4129) @@ -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/dakota/plot_sampdist_bars.m =================================================================== --- /issm/trunk/src/m/dakota/plot_sampdist_bars.m (revision 4129) +++ /issm/trunk/src/m/dakota/plot_sampdist_bars.m (revision 4129) @@ -0,0 +1,205 @@ +% +% plot a stacked bar chart of the sample distributions. +% +% []=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) +% 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) +% +% and the optional fields of dresp are: +% min (double, minimum of sample) +% quart1 (double, first quartile of sample) +% median (double, median of sample) +% quart3 (double, third quartile of sample) +% 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 +% are stacked by the four quartiles, and annotates it with +% the description. the quartiles will be calculated from the +% samples if they do not already exist. +% +% this data would typically be contained in the dakota tabular +% output file 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_sampdist_bars(varargin) + +if ~nargin + help plot_sampdist_bars + 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) + +if ~isfield(dresp,'min') || ~isfield(dresp,'quart1') || ... + ~isfield(dresp,'median') || ~isfield(dresp,'quart3') || ... + ~isfield(dresp,'max') + for i=1:length(dresp) + dresp(i).min =min (dresp(i).sample); + dresp(i).quart1=prctile(dresp(i).sample,25); + dresp(i).median=median (dresp(i).sample); + dresp(i).quart3=prctile(dresp(i).sample,75); + dresp(i).max =max (dresp(i).sample); + end +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) + descr(i)=cellstr(dresp(i).descriptor); + data(i,1)=dresp(i).min; + data(i,2)=dresp(i).quart1-dresp(i).min; + data(i,3)=dresp(i).median-dresp(i).quart1; + data(i,4)=dresp(i).quart3-dresp(i).median; + data(i,5)=dresp(i).max -dresp(i).quart3; +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(dresp) == 1 + data=[data; data]; +end + +figure +hl1=bar(data,'stacked'); +% set barseries properties for lowest value +whitebg('white') +set(hl1(1),'FaceColor','white') +set(hl1(1),'Visible','off') + +ax1=gca; +if length(dresp) == 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('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') + +if ~exist('lstr','var') || isempty(lstr) + lstr={'minimum' 'quartile 1' 'median' 'quartile 3' 'maximum'}; +end + +hleg1=legend(ax1,lstr,'Location','EastOutside',... + 'Orientation','vertical','Interpreter','none'); + +end Index: /issm/trunk/src/m/dakota/postqmu.m =================================================================== --- /issm/trunk/src/m/dakota/postqmu.m (revision 4129) +++ /issm/trunk/src/m/dakota/postqmu.m (revision 4129) @@ -0,0 +1,65 @@ +function md=postqmu(md) + +%INPUT function md=postqmu(md,qmufile,qmudir) +%Deal with dakota output results in files. + +global ISSM_DIR; + +% check to see if dakota returned errors in the err file +qmuerrfile=[md.name '.qmu.err']; + +if exist(qmuerrfile,'file') + fide=fopen(qmuerrfile,'r'); + fline=fgetl(fide); + if ischar(fline) + while ischar(fline) + disp(sprintf('%s',fline)); + fline=fgetl(fide); + end + status=fclose(fide); + cd ../ + error(['Dakota returned error in ''' qmuerrfile ' file. ''' qmudir ''' directory retained.']) + end + status=fclose(fide); +else + cd ../ + error(['Dakota did not generate ''' qmuerrfile ' file. ''' qmudir ''' directory retained.']) +end + +%parse inputs and results from dakota +qmuinfile=[md.name '.qmu.in']; +qmuoutfile=[md.name '.qmu.out']; + +%[method,dvar,dresp_in]=dakota_in_parse(qmuinfile); +%dakotaresults.method =method; +%dakotaresults.dvar =dvar; +%dakotaresults.dresp_in =dresp_in; + +[method,dresp_out,scm,pcm,srcm,prcm]=dakota_out_parse(qmuoutfile); +dakotaresults.dresp_out=dresp_out; +dakotaresults.scm =scm; +dakotaresults.pcm =pcm; +dakotaresults.srcm =srcm; +dakotaresults.prcm =prcm; + +if exist('dakota_tabular.dat','file') + [method,dresp_dat ]=dakota_out_parse('dakota_tabular.dat'); + dakotaresults.dresp_dat=dresp_dat; +end + +%put dakotaresults in their right location. +md.results.dakota=dakotaresults; + +%save input and output files into model +%md.dakotain =readfile([qmufile '.in']); +%md.dakotaout=readfile(qmuoutfile); +%if exist('dakota_tabular.dat','file') +% md.dakotadat=readfile('dakota_tabular.dat'); +%end + +% move all the individual function evalutations into zip files +if ~md.qmu_analysis, + system('zip -mq params.in.zip params.in.[1-9]*'); + system('zip -mq results.out.zip results.out.[1-9]*'); + system('zip -mq matlab.out.zip matlab*.out.[1-9]*'); +end Index: /issm/trunk/src/m/dakota/preqmu.m =================================================================== --- /issm/trunk/src/m/dakota/preqmu.m (revision 4129) +++ /issm/trunk/src/m/dakota/preqmu.m (revision 4129) @@ -0,0 +1,111 @@ +function md=preqmu(md,options) +%QMU - apply Quantification of Margins and Uncertainties techniques +% to a solution sequence (like diagnostic.m, progonstic.m, etc ...), +% using the Dakota software from Sandia. +% +% options come from the solve.m routine. They can include Dakota options: +% +% qmudir: any directory where to run the qmu analysis +% qmufile: input file for Dakota +% ivar: selection number for variables input (if several are specified in variables) +% iresp: same thing for response functions +% imethod: same thing for methods +% iparams: same thing for params +% overwrite: overwrite qmudir +% outfiles: (John?) +% rstfile: backup file name +% rundakota: (John?) +% runmpi: (John?) + +global ISSM_DIR; + +displaystring(md.verbose,'\n%s\n','preprocessing dakota inputs'); + +%first create temporary directory in which we will work +if exist(options.qmudir,'dir') + if ~isfield(options,'overwrite') + %if exist(options.qmudir)==7, + % options.overwrite=input(['Overwrite existing ''' options.qmudir ''' directory? Y/N [N]: '], 's'); + %else + options.overwrite='y'; + %end + + end + if strncmpi(options.overwrite,'y',1) + system(['rm -rf ' options.qmudir '/*']); +% else +% error('Existing ''%s'' directory not overwritten.',qmudir); + end +end +mkdir(options.qmudir) +cd(options.qmudir) + +%when running in library mode, the in file needs to be called md.name.qmu.in +options.qmufile=[md.name ]; + +%retrieve variables and resposnes for this particular analysis. +variables=md.variables(options.ivar); +responses=md.responses(options.iresp); + +%create m and in files for dakota +dakota_in_data(md.qmu_method(options.imethod),variables,md.responses,md.qmu_params(options.iparams),options.qmufile,md); + +%in library mode, we only need the dakota in file +system(['rm -rf ' md.name '.m']); + +%figure out number of variables and responses, it's not straightforwared +numvariables=0; +variable_fieldnames=fieldnames(variables); +for i=1:length(variable_fieldnames), + field_name=variable_fieldnames{i}; + fieldvariables=variables.(field_name); + numvariables=numvariables+numel(variables.(field_name)); +end + +numresponses=0; +response_fieldnames=fieldnames(responses); +for i=1:length(response_fieldnames), + field_name=response_fieldnames{i}; + fieldresponses=responses.(field_name); + numresponses=numresponses+numel(responses.(field_name)); +end + +%ok, now, for this particular qmu analysis, iresp and ivar specifiy the variables and responses. +%The Qmu module will need a list of variable descriptors and response descriptors. +%For ease of use, we gather this list here. + +count=0; +variable_fieldnames=fieldnames(variables); +variabledescriptors={}; +for i=1:length(variable_fieldnames), + field_name=variable_fieldnames{i}; + fieldvariables=variables.(field_name); + for j=1:length(fieldvariables), + descriptor=fieldvariables(j).descriptor; + variabledescriptors{end+1}=descriptor; + count=count+1; + end +end + +count=0; +response_fieldnames=fieldnames(responses); +responsedescriptors={}; +for i=1:length(response_fieldnames), + field_name=response_fieldnames{i}; + fieldresponses=responses.(field_name); + for j=1:length(fieldresponses), + descriptor=fieldresponses(j).descriptor; + responsedescriptors{end+1}=descriptor; + count=count+1; + end +end + +%register the fields that will be needed by the Qmu model. +md.numberofvariables=numvariables; +md.numberofresponses=numresponses; +md.variabledescriptors=variabledescriptors; +md.responsedescriptors=responsedescriptors; + +%now, we have to provide all the info necessary for the solutions to compute the responses. For ex, if mass_flux +%is a response, we need a profile of points. For max_vel, we don't need anything. +md=process_qmu_response_data(md); Index: /issm/trunk/src/m/dakota/process_qmu_options.m =================================================================== --- /issm/trunk/src/m/dakota/process_qmu_options.m (revision 4129) +++ /issm/trunk/src/m/dakota/process_qmu_options.m (revision 4129) @@ -0,0 +1,129 @@ +function outoptions=process_qmu_options(options) +%PROCESS_QMU_OPTIONS - set up default options for qmu phase +% +% Usage: +% options=process_qmu_options(options) +% +% See also: QMU,RECOVER_QMU_OPTIONS + +%analysis_type: check on this option, error out otherwise +found=0; +for i=1:size(options,1), + if strcmpi(options{i,1},'analysis_type'), + analysis_type=options{i,2}; + found=1; + end +end +if ~found, + error('recover_qmu_options error message: no ''analysis_type'' was provided'); +end + +%package: is there one? default to ''JPL'' +found=0; +for i=1:size(options,1), + if strcmpi(options{i,1},'package'), + package=options{i,2}; + found=1; + end +end +if ~found, + disp('recover_qmu_options info message: no ''package'' was provided, defaulting to ''JPL'''); + options(end+1,:)={'package' 'JPL'}; + package='JPL'; +end + +if ~ischar(package), + error(['process_qmu_options error message: package ' package ' not supported yet']); +end + +%sub_analysis_type: check on it, not mandatory +found=0; +for i=1:size(options,1), + if strcmpi(options{i,1},'sub_analysis_type'), + sub_analysis_type=options{i,2}; + found=1; + end +end +if ~found + if ~strcmpi(analysis_type,'thermal'), + sub_analysis_type='none'; + else + disp('recover_qmu_options info message: no ''sub_analysis_type'' was provided, defaulting to ''steady'''); + sub_analysis_type='steady'; + end +end + +%check solution type is supported +if ~(strcmpi(analysis_type,'control') | ... + strcmpi(analysis_type,'diagnostic') | ... + strcmpi(analysis_type,'prognostic') | ... + strcmpi(analysis_type,'thermal') | ... + strcmpi(analysis_type,'parameters') | ... + strcmpi(analysis_type,'mesh') | ... + strcmpi(analysis_type,'mesh2grid') | ... + strcmpi(analysis_type,'transient') ), + error(['process_qmu_options error message: analysis_type ' analysis_type ' not supported yet!']); +end +if ~(strcmpi(sub_analysis_type,'none') | ... + strcmpi(sub_analysis_type,'steady') | ... + strcmpi(sub_analysis_type,'horiz') | ... + strcmpi(sub_analysis_type,'vert') | ... + strcmpi(sub_analysis_type,'') | ... + strcmpi(sub_analysis_type,'transient') ), + error(['process_qmu_options error message: sub_analysis_type ' sub_analysis_type ' not supported yet!']); +end + + + +% process qmu arguments + +%first, the defaults +qmudir ='qmu';% qmudir =['qmu_' datestr(now,'yyyymmdd_HHMMSS')]; +qmufile='qmu'; +ivar =1; +iresp =1; +imethod=1; +iparams=1; +runmpi =false; + +for i=1:size(options,1), + switch options{i,1}, + case 'qmudir' + qmudir=options{i,2}; + case 'qmufile' + qmufile=options{i,2}; + case 'ivar' + ivar=options{i,2}; + case 'iresp' + iresp=options{i,2}; + case 'imethod' + imethod=options{i,2}; + case 'iparams' + iparams=options{i,2}; + case 'overwrite' + outoptions.overwrite=options{i,2}; + case 'outfiles' + outoptions.outfiles=options{i,2}; + case 'rstfile' + outoptions.rstfile=options{i,2}; + case 'rundakota' + outoptions.rundakota=options{i,2}; + case 'runmpi' + runmpi=options{i,2}; + otherwise + %nothing + end +end + +%setup final options structure +outoptions.analysis_type=analysis_type; +outoptions.package=package; +outoptions.sub_analysis_type=sub_analysis_type; +outoptions.qmudir=qmudir; +outoptions.qmufile=qmufile; +outoptions.ivar=ivar; +outoptions.iresp=iresp; +outoptions.imethod=imethod; +outoptions.iparams=iparams; +outoptions.runmpi=runmpi; + Index: /issm/trunk/src/m/dakota/process_qmu_response_data.m =================================================================== --- /issm/trunk/src/m/dakota/process_qmu_response_data.m (revision 4129) +++ /issm/trunk/src/m/dakota/process_qmu_response_data.m (revision 4129) @@ -0,0 +1,22 @@ +function md=process_qmu_response_data(md) +%PROCESS_QMU_RESPONSE_DATA - process any data necessary for the solutions to process the data. +% +% Usage: md=process_qmu_response_data(md) +% +% See also PREQMU, PRESOLVE + + + +for i=1:numel(md.responsedescriptors), + if strcmpi(md.responsedescriptors{i},'mass_flux'), + %we need a profile of points on which to compute the mass_flux, is it here? + if isnans(md.qmu_mass_flux_profile), + error('process_qmu_response_data error message: could not find a mass_flux exp profile!'); + end + if ~ischar(md.qmu_mass_flux_profile), + error('process_qmu_response_data error message: mass_flux exp profile should be a domain outline name'); + end + %ok, process the qmu_mass_flux_profile to build a list of segments: + md.qmu_mass_flux_segments=MassFluxProcessProfile(md); + end +end Index: /issm/trunk/src/m/dakota/qmuisdistributed.m =================================================================== --- /issm/trunk/src/m/dakota/qmuisdistributed.m (revision 4129) +++ /issm/trunk/src/m/dakota/qmuisdistributed.m (revision 4129) @@ -0,0 +1,12 @@ +function found=qmuisdistribted(string) +%QMUISDISTRIBTED - figure out if a string is a decriptor with a numerical postfix. Like thickness1, or drag10 + + +%just take last string element, and see if it is numeric. +last=string(end); + +if ((double(last)<=57) & (double(last)>=49)), + found=1; +else + found=0; +end Index: /issm/trunk/src/m/dakota/qmumarshall.m =================================================================== --- /issm/trunk/src/m/dakota/qmumarshall.m (revision 4129) +++ /issm/trunk/src/m/dakota/qmumarshall.m (revision 4129) @@ -0,0 +1,99 @@ +function qmumarshall(md,variables,responses) +%QMUMARSHALL - output ISSM compatible binary file with qmu fields. This is +% in addition to the marshall routine for regular solve routines. +% Usage: +% qmumarshall(md,variables,responses) +% +% where variables and responses are the Dakota variables and responses found in the model @md. + +%some checks on list of arguments +if ((nargin~=3) & (nargout~=0)) + qmumarshallusage; + error('marshall error message'); +end + +disp(['qmu marshalling file ' md.name '.bin']); + +%open file for binary adding +fid=fopen([ md.name '.bin'],'ab'); +if fid==-1, + error(['qmumarshall error message: could not open ' [md.name '.bin'],' file for binary adding']); +end + +%count how many variables we have +numvariables=0; +variable_fieldnames=fieldnames(variables); +for i=1:length(variable_fieldnames), + field_name=variable_fieldnames{i}; + fieldvariables=variables.(field_name); + numvariables=numvariables+numel(variables.(field_name)); +end +%write number of variables to disk +WriteData(fid,numvariables,'Integer','numberofvariables'); + +%now, for each variable, write descriptor +count=0; +for i=1:length(variable_fieldnames), + field_name=variable_fieldnames{i}; + fieldvariables=variables.(field_name); + for j=1:length(fieldvariables), + descriptor=fieldvariables(j).descriptor; + WriteData(fid,descriptor,'String',['variabledescriptor' num2str(count)]); + count=count+1; + end +end + +%deal with responses + +%count how many responses we have +numresponses=0; +response_fieldnames=fieldnames(responses); +for i=1:length(response_fieldnames), + field_name=response_fieldnames{i}; + fieldresponses=responses.(field_name); + numresponses=numresponses+numel(responses.(field_name)); +end +%write number of responses to disk +WriteData(fid,numresponses,'Integer','numberofresponses'); + +%now, for each response, write descriptor +count=0; +for i=1:length(response_fieldnames), + field_name=response_fieldnames{i}; + fieldresponses=responses.(field_name); + for j=1:length(fieldresponses), + descriptor=fieldresponses(j).descriptor; + WriteData(fid,descriptor,'String',['responsedescriptor' num2str(count)]); + count=count+1; + end +end + +%write response specific data +count=0; +for i=1:length(response_fieldnames), + field_name=response_fieldnames{i}; + fieldresponses=responses.(field_name); + for j=1:length(fieldresponses), + descriptor=fieldresponses(j).descriptor; + if strcmpi(descriptor,'mass_flux'), + WriteData(fid,md.qmu_mass_flux_segments,'Mat','qmu_mass_flux_segments'); + end + end +end + +%write part and npart to disk +WriteData(fid,md.npart,'Integer','npart'); +WriteData(fid,md.part,'Mat','part'); + +%close file +st=fclose(fid); +if st==-1, + error(['qmumarshall error message: could not close file ' [md.name '.bin']]); +end + +end + +function qmumarshallusage(); +disp(' '); +disp('function qmumarshall(md,variables,responses)'); +end Index: /issm/trunk/src/m/dakota/qmuname.m =================================================================== --- /issm/trunk/src/m/dakota/qmuname.m (revision 4129) +++ /issm/trunk/src/m/dakota/qmuname.m (revision 4129) @@ -0,0 +1,15 @@ +function md=qmuname(md,varargin) +%INPUT function md=qmuname(md) +%Pick up the number from a file, or get it directly from the Dakota structure. Then modify the name of this +%model to reflect this new number. + +if nargin==1, + fid=fopen('number','r'); + number=fscanf(fid,'%i',1) + fclose(fid); +else + number=varargin{1}; +end + +%modify model name by appending number to the name +md.name=[md.name num2str(number)]; Index: /issm/trunk/src/m/dakota/qmuresponse.m =================================================================== --- /issm/trunk/src/m/dakota/qmuresponse.m (revision 4129) +++ /issm/trunk/src/m/dakota/qmuresponse.m (revision 4129) @@ -0,0 +1,35 @@ +function response=qmuresponse(models,results,processedresults,descriptor) +%QMURESPONSE - compute response function from model results. + +if strcmpi(descriptor,'max_vel'), + response=max(processedresults.vel); +elseif strcmpi(descriptor,'min_vel'), + response=min(processedresults.vel); +elseif strcmpi(descriptor,'max_vx'), + response=max(processedresults.vx); +elseif strcmpi(descriptor,'min_vx'), + response=min(processedresults.vx); +elseif strcmpi(descriptor,'max_vy'), + response=max(processedresults.vy); +elseif strcmpi(descriptor,'min_vy'), + response=min(processedresults.vy); +elseif strcmpi(descriptor,'mass_flux'), + %call mass flux module. + m_dh=models.dh; + m_dhu=models.dhu; + m_ds=models.ds; + ishutter=m_dh.parameters.ishutter; + ismacayealpattyn=m_dhu.parameters.ismacayealpattyn; + isstokes=m_ds.parameters.isstokes; + if ishutter, + response=MassFlux(m_dhu.elements,m_dhu.nodes,m_dhu.vertices,m_dhu.loads,m_dhu.materials,m_dhu.parameters,results.u_g); + elseif ismacayealpattyn, + response=MassFlux(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,results.u_g); + elseif isstokes, + response=MassFlux(m_ds.elements,m_ds.nodes,m_ds.vertices,m_ds.loads,m_ds.materials,m_ds.parameters,results.u_g); + else + error('qmuresponse error message: unsupported analysis type for mass_flux computation!'); + end +else + error(['qmuresponse error message: unknow descriptor ' descriptor]); +end Index: /issm/trunk/src/m/dakota/qmuroot.m =================================================================== --- /issm/trunk/src/m/dakota/qmuroot.m (revision 4129) +++ /issm/trunk/src/m/dakota/qmuroot.m (revision 4129) @@ -0,0 +1,13 @@ +function root=qmuroot(string) +%QMUROOT - return root of a distributed descriptor + +root=''; +found=0; +for i=1:length(string), + if ((49<=double(string(i))) && (double(string(i)<=57))) + break; + else + root=[root string(i)]; + end +end + Index: /issm/trunk/src/m/dakota/recover_qmu_options.m =================================================================== --- /issm/trunk/src/m/dakota/recover_qmu_options.m (revision 4129) +++ /issm/trunk/src/m/dakota/recover_qmu_options.m (revision 4129) @@ -0,0 +1,22 @@ +function options=recover_qmu_options(md,varargin) +%RECOVER_SOLVE_OPTIONS - recover solution options for qmu runs. +% +% Usage: +% options=recover_qmu_options(md,varargin); +% +% See also: SOLVE + +%initialize options. +options=cell(0,2); + +%make sure length(varargin) is even, ie options come in pairs. +if mod(length(varargin),2), + error('recover_qmu_options error message: an even number of options is necessary'); +end + +%go through varargin, extract options +for i=1:length(varargin)/2, + + options(end+1,:)={varargin{2*i-1} varargin{2*i}}; + +end Index: /issm/trunk/src/m/dakota/responsedataprocessing/MassFluxProcessProfile.m =================================================================== --- /issm/trunk/src/m/dakota/responsedataprocessing/MassFluxProcessProfile.m (revision 4129) +++ /issm/trunk/src/m/dakota/responsedataprocessing/MassFluxProcessProfile.m (revision 4129) @@ -0,0 +1,14 @@ +function segments=MassFluxProcessProfile(md); +%MASSFLUXPROCESSPROFILE: process an argus domain outlien profile into a list of segments. +% +% Usage: segments=MassFluxProcessProfile(md); +% +% +% See also: PROCESS_QMU_RESPONSE_DATA, PREQMU + + +%first read the profile points. +profile=expread(['../' md.qmu_mass_flux_profile],1); %go to ../ because the exp file is in the root directory. + +%project this profile onto mesh. +segments=ProfileProjectOntoMesh(md,profile); Index: /issm/trunk/src/m/dakota/setupdesign/QmuSetupDesign.m =================================================================== --- /issm/trunk/src/m/dakota/setupdesign/QmuSetupDesign.m (revision 4129) +++ /issm/trunk/src/m/dakota/setupdesign/QmuSetupDesign.m (revision 4129) @@ -0,0 +1,51 @@ +function dvar=QmuSetupDesign(dvar,variables,params,varargin) + +%get descriptor +descriptor=variables.descriptor; + +%loop on descriptor +if strcmpi(descriptor,'rho_ice') + + dvar=setuprhoice(dvar,variables,params,varargin{:}); + +elseif strcmpi(descriptor,'rho_water') + + dvar=setuprhowater(dvar,variables,params,varargin{:}); + +elseif strcmpi(descriptor,'heatcapacity') + + dvar=setupheatcapacity(dvar,variables,params,varargin{:}); + +elseif strcmpi(descriptor,'thermalconductivity') + + dvar=setupthermalconductivity(dvar,variables,params,varargin{:}); + +elseif strcmpi(descriptor,'gravity') + + dvar=setupgravity(dvar,variables,params,varargin{:}); + +elseif strcmpi(descriptor,'thickness') + + dvar=setupthickness(dvar,variables,params,varargin{:}); + +elseif strcmpi(descriptor,'drag') + + dvar=setupdrag(dvar,variables,params,varargin{:}); + +elseif strncmpi(descriptor,'drag_node',9) + + dvar=setupdrag_node(dvar,variables,params,varargin{:}); + +elseif strncmpi(descriptor,'thickness_node',14) + + dvar=setupthickness_node(dvar,variables,params,varargin{:}); + +elseif strcmpi(descriptor,'riftsfriction') + + dvar=setupriftsfriction(dvar,variables,params,varargin{:}); + +else + error(['QmuSetupDesign warning message: could not find ' descriptor ' setup design function']); +end + +end Index: /issm/trunk/src/m/dakota/setupdesign/setupdrag.m =================================================================== --- /issm/trunk/src/m/dakota/setupdesign/setupdrag.m (revision 4129) +++ /issm/trunk/src/m/dakota/setupdesign/setupdrag.m (revision 4129) @@ -0,0 +1,28 @@ +function dvar=setupdrag(dvar,variables,params,varargin) + +for i=1:length(varargin) + if strcmp(class(varargin{i}),'model') + md=varargin{i}; + break; + end +end + +%capture stddev: +stddev=variables.stddev; + +if length(stddev)>md.npart, + error('setupdrag error message: stddev should be either a scalar or a ''npart'' length vector'); +end + +%ok, dealing with semi-discrete distributed variable. Distribute according to how many +%partitions we want + +for j=1:md.npart + dvar(end+1) =variables; + dvar(end ).descriptor=sprintf('%s%d',variables.descriptor,j); + if length(stddev)>1, + dvar(end ).stddev=stddev(j); + end +end + +end Index: /issm/trunk/src/m/dakota/setupdesign/setupdrag_node.m =================================================================== --- /issm/trunk/src/m/dakota/setupdesign/setupdrag_node.m (revision 4129) +++ /issm/trunk/src/m/dakota/setupdesign/setupdrag_node.m (revision 4129) @@ -0,0 +1,35 @@ +function dvar=setupdrag_node(dvar,variables,params,varargin) + +for i=1:length(varargin) + if strcmp(class(varargin{i}),'model') + md=varargin{i}; + break; + end +end + +%first, recover id of node: +descriptor=variables.descriptor(); +node=descriptor(10:end); + +%this variable we are trying to setup for qmu analysis is shadowing +%a variable in another pool of distributed variables. Ex: drag_node2 should +%owershadow the drag2 parameter. So go through the dvars, and figure out +%if the drag was already distributed. If not, error out. It so, then +%plug drag_nodei variable into the i'th variable of the distributed drag. + + +found=0; +for i=1:numel(dvar), + if strcmpi(dvar(i).descriptor,['drag' num2str(node)]), + found=i; + break; + end +end +if found==0, + error('setupdrag_node error message: could not find distributed drag parameters!'); +end + +%overshadow dvar(found) variable. +variables.descriptor=['drag' num2str(node)]; +dvar(found)=variables; + Index: /issm/trunk/src/m/dakota/setupdesign/setupgravity.m =================================================================== --- /issm/trunk/src/m/dakota/setupdesign/setupgravity.m (revision 4129) +++ /issm/trunk/src/m/dakota/setupdesign/setupgravity.m (revision 4129) @@ -0,0 +1,5 @@ +function dvar=setupgravity(dvar,variables,params,varargin) + +dvar(end+1)=variables; + +end Index: /issm/trunk/src/m/dakota/setupdesign/setupheatcapacity.m =================================================================== --- /issm/trunk/src/m/dakota/setupdesign/setupheatcapacity.m (revision 4129) +++ /issm/trunk/src/m/dakota/setupdesign/setupheatcapacity.m (revision 4129) @@ -0,0 +1,5 @@ +function dvar=setupheatcapacity(dvar,variables,params,varargin) + +dvar(end+1)=variables; + +end Index: /issm/trunk/src/m/dakota/setupdesign/setuprhoice.m =================================================================== --- /issm/trunk/src/m/dakota/setupdesign/setuprhoice.m (revision 4129) +++ /issm/trunk/src/m/dakota/setupdesign/setuprhoice.m (revision 4129) @@ -0,0 +1,5 @@ +function dvar=setuprhoice(dvar,variables,params,varargin) + +dvar(end+1)=variables; + +end Index: /issm/trunk/src/m/dakota/setupdesign/setuprhowater.m =================================================================== --- /issm/trunk/src/m/dakota/setupdesign/setuprhowater.m (revision 4129) +++ /issm/trunk/src/m/dakota/setupdesign/setuprhowater.m (revision 4129) @@ -0,0 +1,5 @@ +function dvar=setuprhowater(dvar,variables,params,varargin) + +dvar(end+1)=variables; + +end Index: /issm/trunk/src/m/dakota/setupdesign/setupriftsfriction.m =================================================================== --- /issm/trunk/src/m/dakota/setupdesign/setupriftsfriction.m (revision 4129) +++ /issm/trunk/src/m/dakota/setupdesign/setupriftsfriction.m (revision 4129) @@ -0,0 +1,17 @@ +function dvar=setupriftsfriction(dvar,variables,params,varargin) + +for i=1:length(varargin) + if strcmp(class(varargin{i}),'model') + md=varargin{i}; + break; + end +end + +%we have several rifts. + +for j=1:md.numrifts + dvar(end+1) =variables; + dvar(end ).descriptor=sprintf('%s%d',variables.descriptor,j); +end + +end Index: /issm/trunk/src/m/dakota/setupdesign/setupthermalconductivity.m =================================================================== --- /issm/trunk/src/m/dakota/setupdesign/setupthermalconductivity.m (revision 4129) +++ /issm/trunk/src/m/dakota/setupdesign/setupthermalconductivity.m (revision 4129) @@ -0,0 +1,5 @@ +function dvar=setupthermalconductivity(dvar,variables,params,varargin) + +dvar(end+1)=variables; + +end Index: /issm/trunk/src/m/dakota/setupdesign/setupthickness.m =================================================================== --- /issm/trunk/src/m/dakota/setupdesign/setupthickness.m (revision 4129) +++ /issm/trunk/src/m/dakota/setupdesign/setupthickness.m (revision 4129) @@ -0,0 +1,27 @@ +function dvar=setupthickness(dvar,variables,params,varargin) + +for i=1:length(varargin) + if strcmp(class(varargin{i}),'model') + md=varargin{i}; + break; + end +end +%capture stddev: +stddev=variables.stddev; + +if length(stddev)>md.npart, + error('setupthickness error message: stddev should be either a scalar or a ''npart'' length vector'); +end + +%ok, dealing with semi-discrete distributed variable. Distribute according to how many +%partitions we want + +for j=1:md.npart + dvar(end+1) =variables; + dvar(end ).descriptor=sprintf('%s%d',variables.descriptor,j); + if length(stddev)>1, + dvar(end ).stddev=stddev(j); + end +end + +end Index: /issm/trunk/src/m/dakota/setupdesign/setupthickness_node.m =================================================================== --- /issm/trunk/src/m/dakota/setupdesign/setupthickness_node.m (revision 4129) +++ /issm/trunk/src/m/dakota/setupdesign/setupthickness_node.m (revision 4129) @@ -0,0 +1,34 @@ +function dvar=setupthickness_node(dvar,variables,params,varargin) + +for i=1:length(varargin) + if strcmp(class(varargin{i}),'model') + md=varargin{i}; + break; + end +end + +%first, recover id of node: +descriptor=variables.descriptor(); +node=descriptor(15:end); + +%this variable we are trying to setup for qmu analysis is shadowing +%a variable in another pool of distributed variables. Ex: thickness_node2 should +%owershadow the thickness2 parameter. So go through the dvars, and figure out +%if the thickness was already distributed. If not, error out. It so, then +%plug thickness_nodei variable into the i'th variable of the distributed thickness. + + +found=0; +for i=1:numel(dvar), + if strcmpi(dvar(i).descriptor,['thickness' num2str(node)]), + found=i; + break; + end +end +if found==0, + error('setupthickness_node error message: could not find distributed thickness parameters!'); +end + +%overshadow dvar(found) variable. +variables.descriptor=['thickness' num2str(node)]; +dvar(found)=variables;