Changeset 32

Timestamp:

04/24/09 10:38:58 (16 years ago)

Author:

Mathieu Morlighem

Message:

initial input from ice1

Location:

issm/trunk/src/m/solutions

Files:

• dakota/dakota_in_data.m (modified) (3 diffs)
• dakota/dakota_in_params.m (modified) (5 diffs)
• dakota/dakota_in_parse.m (modified) (10 diffs)
• dakota/dakota_in_write.m (modified) (20 diffs)
• dakota/dakota_m_write.m (modified) (12 diffs)
• dakota/qmu.m (modified) (3 diffs)
• dakota/setupdesign/QmuSetupDesign.m (modified) (2 diffs)
• dakota/setupdesign/setupdrag.m (modified) (1 diff)
• dakota/setupdesign/setupgravity.m (added)
• dakota/setupdesign/setupheatcapacity.m (added)
• dakota/setupdesign/setuprhoice.m (modified) (1 diff)
• dakota/setupdesign/setuprhowater.m (added)
• dakota/setupdesign/setupriftsfriction.m (modified) (1 diff)
• dakota/setupdesign/setupthermalconductivity.m (added)
• dakota/setupdesign/setupthickness.m (modified) (1 diff)
• dakota/updatefunctions/QmuUpdateFunctions.m (modified) (1 diff)
• dakota/updatefunctions/updatedrag.m (modified) (2 diffs)
• dakota/updatefunctions/updategravity.m (modified) (1 diff)
• dakota/updatefunctions/updateheatcapacity.m (modified) (1 diff)
• dakota/updatefunctions/updaterho_ice.m (modified) (1 diff)
• dakota/updatefunctions/updaterho_water.m (modified) (1 diff)
• dakota/updatefunctions/updateriftsfriction.m (modified) (2 diffs)
• dakota/updatefunctions/updatethermalconductivity.m (modified) (1 diff)
• dakota/updatefunctions/updatethickness.m (modified) (2 diffs)
• ice/ModelProcessorDiagnosticHoriz.m (modified) (2 diffs)
• ice/ModelProcessorPrognostic.m (modified) (1 diff)
• ice/StrainRateCompute.m (modified) (3 diffs)
• ice/icediagnostic.m (added)
• ice/icediagnostic_wrapper.m (added)
• macayeal/macayealcontrol.m (modified) (6 diffs)

issm/trunk/src/m/solutions/dakota/dakota_in_data.m

@@ -2 +2 @@
 %  define the data to write the dakota .in file
 %
-%  []=dakota_in_data(md,filei,package)
+%  []=dakota_in_data(variables,responses,dmeth,dparams,filei,package,varargin)
 %
-function []=dakota_in_data(md,filei,package)
+function []=dakota_in_data(dmeth,variables,responses,dparams,filei,package,varargin)
 
 if ~nargin
@@ -11 +11 @@
 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 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(md.variables);
+fnames=fieldnames(variables);
+
 for i=1:length(fnames)
-    fhandle=str2func([class(md.variables.(fnames{i})) '.empty']);
-    dvar.(fnames{i})=fhandle();
-    for j=1:length(md.variables.(fnames{i}))
-        %call setupdesign
-        dvar.(fnames{i})=QmuSetupDesign(md,dvar.(fnames{i}),md.variables.(fnames{i})(j)); 
+    
+%  for linear constraints, just copy
+
+    if strcmp(class(variables.(fnames{i})),'linear_inequality_constraint') || ...
+       strcmp(class(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
@@ -25 +68 @@
 %%  responses
 
-fnames=fieldnames(md.responses);
+fnames=fieldnames(responses);
+
 for i=1:length(fnames)
-    fhandle=str2func([class(md.responses.(fnames{i})) '.empty']);
-    dresp.(fnames{i})=fhandle();
-    for j=1:length(md.responses.(fnames{i}))
-        dresp.(fnames{i})(j)=md.responses.(fnames{i})(j);
-    end
+%     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
 
-%%  run parameters
-
-params=dakota_in_params(struct());
-
-params.evaluation_concurrency=md.evaluation_concurrency;
-params.npart=md.npart;
-disp(md.method);
-
-if     strncmpi(md.method,'nond_l',6),
-        params.method               ='nond_local_reliability';
-        params.interval_type        =md.interval_type;
-        params.fd_gradient_step_size=md.fd_gradient_step_size;
-elseif strncmpi(md.method,'nond_s',6),
-        params.method               ='nond_sampling';
-        params.seed                 =md.seed;
-        params.samples              =md.samples;
-        params.sample_type          =md.sample_type;
-elseif strncmpi(md.method,'conmin_f',8),
-        params.method               ='conmin_frcg';
-elseif strncmpi(md.method,'conmin_m',8),
-        params.method               ='conmin_mfd';
-else
-        error('dakota_in_data: error message, method not supported yet!');
-end
-
-if     strcmpi(md.analysis_driver,'matlab')
-    params.analysis_driver='matlab';
-    if ~isempty(md.analysis_components)
-        params.analysis_components=md.analysis_components;
-    else
-        params.analysis_components=filei;
-    end
-elseif ~isempty(md.analysis_driver)
-    params.analysis_driver=md.analysis_driver;
-else
-    md.analysis_driver=params.analysis_driver;
-end
+%%  write files
 
 %Write in file
-dakota_in_write(params.method,dvar,dresp,filei,params);
+dakota_in_write(dmeth.method,dmeth,dvar,dresp,params,filei,varargin{:});
 
 %Write m file
-dakota_m_write(md,params.method,dvar,dresp,filei,params,package);
+dakota_m_write(dmeth.method,dmeth,dvar,dresp,params,filei,package,varargin{:});
 
 end

issm/trunk/src/m/solutions/dakota/dakota_in_params.m

@@ -5 +5 @@
 %
 function [params]=dakota_in_params(params)
-
-global ISSM_DIR;
 
 if ~nargin
@@ -34 +32 @@
 %%  method section
 
-% if isfield(params,'method')
-%     method=params.method;
-% end
+%  all method parameters are in the dakota_method class
 
-%  method-independent
-
-if ~isfield(params,'output')
-    params.output=false;
-end
-if ~isfield(params,'max_iterations')
-    params.max_iterations=100;
-end
-if ~isfield(params,'max_function_evaluations')
-    params.max_function_evaluations=100;
-end
-if ~isfield(params,'convergence_tolerance')
-    params.convergence_tolerance=1.e-4;
-end
-
-%  nondeterministic methods
-
-if ~isfield(params,'distribution')
-    params.distribution='cumulative';
-end
-if ~isfield(params,'compute')
-    params.compute='probabilities';
-end
-
-%  nond_sampling
-
-if ~isfield(params,'seed')
-    params.seed=1234;
-end
-if ~isfield(params,'samples')
-    params.samples=100;
-end
-if ~isfield(params,'sample_type')
-    params.sample_type='lhs';
-end
-if ~isfield(params,'all_variables')
-    params.all_variables=false;
-end
-if ~isfield(params,'variance_based_decomp')
-    params.variance_based_decomp=false;
-end
-
-%  nond_local_reliability
+%%  model section
 
 %%  interface section
@@ -91 +45 @@
 end
 if ~isfield(params,'analysis_driver')
-    params.analysis_driver=[ISSM_DIR '/src/m/solutions/dakota/cielo_ice_script.sh'];
+    params.analysis_driver='';
+end
+if ~isfield(params,'analysis_components')
+    params.analysis_components='';
 end
 if ~isfield(params,'parameters_file')
@@ -105 +62 @@
     params.file_save=true;
 end
-if ~isfield(params,'analysis_components')
-    params.analysis_components='';
-end
 
 %%  responses section
 
-%  gradient specifications dependent on method
-% if ~isfield(params,'no_gradients')
-%     params.no_gradients=false;
-% end
 if ~isfield(params,'numerical_gradients')
     params.numerical_gradients=true;
@@ -140 +90 @@
     params.id_numerical_gradients=false;
 end
-% if ~isfield(params,'no_hessians')
-%     params.no_hessians=false;
-% end
+%  hessians not fully implemented
+if ~isfield(params,'numerical_hessians')
+    params.numerical_gradients=true;
+end
+if ~isfield(params,'hessian_gradient_step_size')
+    params.fd_gradient_step_size=0.001;
+end
 
 end

issm/trunk/src/m/solutions/dakota/dakota_in_parse.m

@@ -108 +108 @@
             ncdv=tlist;
             dvar.cdv=[];
-            display(sprintf('  Number of Dakota %s variables=%d.',...
+            display(sprintf('    Number of Dakota %s variables=%d.',...
                     'continuous_design',ncdv));
         case 'cdv_initial_point'
@@ -135 +135 @@
             nnuv=tlist;
             dvar.nuv=[];
-            display(sprintf('  Number of Dakota %s variables=%d.',...
+            display(sprintf('    Number of Dakota %s variables=%d.',...
                     'normal_uncertain',nnuv));
         case 'nuv_means'
@@ -167 +167 @@
             ncsv=tlist;
             dvar.csv=[];
-            display(sprintf('  Number of Dakota %s variables=%d.',...
+            display(sprintf('    Number of Dakota %s variables=%d.',...
                     'continuous_state',ncsv));
         case 'csv_initial_state'
@@ -236 +236 @@
 dresp=[];
 nof =0;
+nlst=0;
 nnic=0;
 nnec=0;
-nlst=0;
 nrf =0;
 
@@ -262 +262 @@
             nof =tlist;
             dresp.of =[];
-            display(sprintf('  Number of Dakota %s=%d.',...
+            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',nof));
+            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;
-            display(sprintf('  Number of Dakota %s=%d.',...
-                    'nonlinear_equality_constraints',nof));
-        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',nof));
+            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.',...
+            display(sprintf('    Number of Dakota %s=%d.',...
                     'response_functions',nrf));
+
         case 'response_descriptors'
             [tlist,tokens,itoken]=readtlist(fidi,tokens,itoken);
-            for i=1:length(tlist)
-                dresp.rf(i).descriptor=char(tlist(i));
-            end
+            desc=tlist;
         otherwise
             warning('responses_parse:unrec_key',...
@@ -308 +377 @@
        strncmpi(tokens{1}{itoken},'responses',9)
 
-%  supply default descriptors if necessary
-
-        if isfield(dresp,'of' ) && ~isfield(dresp.of,'descriptor')
-            for i=1:nof
-                dresp.of(i).descriptor=sprintf('obj_fn_%d',i);
-            end
-        end
-        if isfield(dresp,'nic') && ~isfield(dresp.nic,'descriptor')
-            for i=1:nnic
-                dresp.nic(i).descriptor=sprintf('nln_ineq_con_%d',i);
-            end
-        end
-        if isfield(dresp,'nec') && ~isfield(dresp.nec,'descriptor')
-            for i=1:nnec
-                dresp.nec(i).descriptor=sprintf('nln_eq_con_%d',i);
-            end
-        end
-        if isfield(dresp,'lst') && ~isfield(dresp.lst,'descriptor')
-            for i=1:nlst
-                dresp.lst(i).descriptor=sprintf('least_sq_term_%d',i);
-            end
-        end
-        if isfield(dresp,'rf' ) && ~isfield(dresp.rf,'descriptor')
-            for i=1:nrf
-                dresp.rf(i).descriptor=sprintf('response_fn_%d',i);
+%  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
@@ -413 +517 @@
 %  assemble the lists by response function
 
-        if exist('nrespl','var')
+        if exist('nrespl','var') && isfield(dresp,'rf')
             if ~exist('nresp','var')
                 nresp(1:length(dresp.rf))=floor(length(nrespl)/length(dresp.rf));
@@ -424 +528 @@
         end
 
-        if exist('nprobl','var')
+        if exist('nprobl','var') && isfield(dresp,'rf')
             if ~exist('nprob','var')
                 nprob(1:length(dresp.rf))=floor(length(nprobl)/length(dresp.rf));
@@ -435 +539 @@
         end
 
-        if exist('nrell' ,'var')
+        if exist('nrell' ,'var') && isfield(dresp,'rf')
             if ~exist('nrel' ,'var')
                 nrel (1:length(dresp.rf))=floor(length(nrell )/length(dresp.rf));
@@ -446 +550 @@
         end
 
-        if exist('ngrell','var')
+        if exist('ngrell','var') && isfield(dresp,'rf')
             if ~exist('ngrel','var')
                 ngrel(1:length(dresp.rf))=floor(length(ngrell)/length(dresp.rf));

issm/trunk/src/m/solutions/dakota/dakota_in_write.m

@@ -2 +2 @@
 %  write a Dakota .in file.
 %
-%  []=dakota_in_write(method,dvar,dresp,filei,params)
+%  []=dakota_in_write(method,dmeth,dvar,dresp,params,filei,varargin)
 %
-function []=dakota_in_write(method,dvar,dresp,filei,params)
+function []=dakota_in_write(method,dmeth,dvar,dresp,params,filei,varargin)
 
 if ~nargin
@@ -13 +13 @@
 %  process the input parameters
 
-if ~exist('method' ,'var') || isempty(method)
-    method=input('Method: sampling (s) or local reliability (l)?  ','s');
-    if     strncmpi(method,'s',1)
-        method='nond_sampling';
-    elseif strncmpi(method,'l',1)
-        method='nond_local_reliability';
-    end
-end
-
-switch method
-%     case {'dot_bfgs','dot_frcg','dot_mmfd','dot_slp','dot_sqp'}
-%     case {'npsol_sqp'}
-    case {'conmin_frcg','conmin_mfd'}
-%     case {'optpp_cg','optpp_q_newton','optpp_fd_newton',...
-%             'optpp_newton','optpp_pds'}
-%     case {'asynch_pattern_search'}
-%     case {'coliny_cobyla','coliny_direct','coliny_ea',...
-%             'coliny_pattern_search','coliny_solis_wets'}
-%     case {'ncsu_direct'}
-%     case {'moga','soga'}
-%     case {'nl2sol','nlssol_sqp','optpp_g_newton'}
-    case {'nond_sampling'}
-    case {'nond_local_reliability'}
-%     case {'dace','fsu_quasi_mc','fsu_cvt'}
-%     case {'vector_parameter_study','list_parameter_study',...
-%             'centered parameter_study','multidim_parameter_study'}
-    otherwise
-        error('Unrecognized method: ''%s''.',method);
+if ~exist('method','var') || isempty(method)
+    method=input('Method?  ','s');
+end
+
+if ~exist('dmeth' ,'var') || isempty(dmeth)
+    dmeth=dakota_method(method);
 end
 
@@ -69 +47 @@
 %  write the method section
 
-method_write(fidi,method,dresp,params);
+method_write(fidi,dmeth,dresp,params);
 
 %  write the model section
@@ -77 +55 @@
 %  write the variables section
 
-variables_write(fidi,method,dvar);
+variables_write(fidi,dmeth,dvar);
 
 %  write the interface section
@@ -85 +63 @@
 %  write the responses section
 
-responses_write(fidi,method,dresp,params);
+responses_write(fidi,dmeth,dresp,params);
 
 fclose(fidi);
@@ -109 +87 @@
 %%  function to write the method section of the file
 
-function []=method_write(fidi,method,dresp,params)
+function []=method_write(fidi,dmeth,dresp,params)
 
 display('Writing method section of Dakota input file.');
 
 fprintf(fidi,'method,\n');
-fprintf(fidi,'\t%s\n',method);
-param_write(fidi,'\t  ','output',' ','\n',params);
-param_write(fidi,'\t  ','max_iterations','           = ','\n',params);
-param_write(fidi,'\t  ','max_function_evaluations',' = ','\n',params);
-param_write(fidi,'\t  ','convergence_tolerance','    = ','\n',params);
-
-switch method
-%     case {'dot_bfgs','dot_frcg','dot_mmfd','dot_slp','dot_sqp'}
-%     case {'npsol_sqp'}
-    case {'conmin_frcg','conmin_mfd'}
-%     case {'optpp_cg','optpp_q_newton','optpp_fd_newton',...
-%             'optpp_newton','optpp_pds'}
-%     case {'asynch_pattern_search'}
-%     case {'coliny_cobyla','coliny_direct','coliny_ea',...
-%             'coliny_pattern_search','coliny_solis_wets'}
-%     case {'ncsu_direct'}
-%     case {'moga','soga'}
-%     case {'nl2sol','nlssol_sqp','optpp_g_newton'}
-    case {'nond_sampling'}
-        param_write(fidi,'\t  ','seed','      = ','\n',params);
-        param_write(fidi,'\t  ','samples','   = ','\n',params);
-        param_write(fidi,'\t  ','sample_type',' ','\n',params);
-    case {'nond_local_reliability'}
-%     case {'dace','fsu_quasi_mc','fsu_cvt'}
-%     case {'vector_parameter_study','list_parameter_study',...
-%             'centered parameter_study','multidim_parameter_study'}
-end
+fprintf(fidi,'\t%s\n',dmeth.method);
+
+dmeth_params_write(dmeth,fidi);
 
 %  write response levels
 
-if isfield(dresp,'rf')
+if strcmp(dmeth.type,'nond') && isfield(dresp,'rf')
     param_write(fidi,'\t  ','distribution',' ','\n',params);
-    [respl,probl,rell,grell]=dresp_levels(dresp.rf);
+    [respl,probl,rell,grell]=prop_levels(dresp.rf);
     if ~isempty(respl)
         rlev_write(fidi,'response_levels',respl);
@@ -198 +152 @@
 %%  function to write the variables section of the file
 
-function []=variables_write(fidi,method,dvar)
+function []=variables_write(fidi,dmeth,dvar)
 
 display('Writing variables section of Dakota input file.');
@@ -206 +160 @@
 %  variables vary by method
 
-switch method
-%     case {'dot_bfgs','dot_frcg','dot_mmfd','dot_slp','dot_sqp'}
-%     case {'npsol_sqp'}
-    case {'conmin_frcg','conmin_mfd'}
-        vsets_write(fidi,dvar,'cdv','csv');
-%     case {'optpp_cg','optpp_q_newton','optpp_fd_newton',...
-%             'optpp_newton','optpp_pds'}
-%     case {'asynch_pattern_search'}
-%     case {'coliny_cobyla','coliny_direct','coliny_ea',...
-%             'coliny_pattern_search','coliny_solis_wets'}
-%     case {'ncsu_direct'}
-%     case {'moga','soga'}
-%     case {'nl2sol','nlssol_sqp','optpp_g_newton'}
-    case {'nond_sampling'}
-        vsets_write(fidi,dvar,'nuv','csv');
-    case {'nond_local_reliability'}
-        vsets_write(fidi,dvar,'nuv','csv');
-%     case {'dace','fsu_quasi_mc','fsu_cvt'}
-%     case {'vector_parameter_study','list_parameter_study',...
-%             'centered parameter_study','multidim_parameter_study'}
-end
+vsets_write(fidi,dvar,dmeth.variables);
+
+%  linear constraints vary by method
+
+lcsets_write(fidi,dvar,dmeth.lcspec);
 
 fprintf(fidi,'\n');
@@ -234 +172 @@
 %%  function to write variable sets
 
-function []=vsets_write(fidi,dvar,varargin)
-
-for i=1:length(varargin)
-    switch(varargin{i})
+function []=vsets_write(fidi,dvar,variables)
+
+for i=1:length(variables)
+    switch(variables{i})
         case 'cdv'
-            vstring='continuous_design';
+            cstring='continuous_design';
         case 'nuv'
-            vstring='normal_uncertain';
+            cstring='normal_uncertain';
         case 'csv'
-            vstring='continuous_state';
+            cstring='continuous_state';
         otherwise
             warning('vsets_write:unrec_var',...
-                'Unrecognized variable ''%s''.',varargin{i});
-    end
-
-    if isfield(dvar,varargin{i})
-        vlist_write(fidi,vstring,varargin{i},dvar.(varargin{i}));
+                'Unrecognized variable ''%s''.',variables{i});
+    end
+
+    if isfield(dvar,variables{i})
+        vlist_write(fidi,cstring,variables{i},dvar.(variables{i}));
     end
 end
@@ -258 +196 @@
 %%  function to write variable list
 
-function []=vlist_write(fidi,vstring,vstring2,dvar)
+function []=vlist_write(fidi,cstring,cstring2,dvar)
 
 %  put variables into lists for writing
 
-vinitpt=dvar_initpt(dvar);
-vlower =dvar_lower (dvar);
-vupper =dvar_upper (dvar);
-vmean  =dvar_mean  (dvar);
-vstddev=dvar_stddev(dvar);
-vinitst=dvar_initst(dvar);
-vdesc  =dvar_desc  (dvar);
+pinitpt=prop_initpt(dvar);
+plower =prop_lower (dvar);
+pupper =prop_upper (dvar);
+pmean  =prop_mean  (dvar);
+pstddev=prop_stddev(dvar);
+pinitst=prop_initst(dvar);
+pstype =prop_stype (dvar);
+pscale =prop_scale (dvar);
+pdesc  =prop_desc  (dvar);
 
 %  write variables
@@ -275 +215 @@
 disp(sprintf('  Writing %d %s variables.',length(dvar),class(dvar)));
 
-fprintf(fidi,'\t%s = %d\n',vstring,length(dvar));
-if ~isempty(vinitpt)
-    fprintf(fidi,'\t  %s_initial_point =\n',vstring2);
-    vector_write(fidi,sprintf('\t    '),vinitpt,6,76);
-end
-if ~isempty(vlower)
-    fprintf(fidi,'\t  %s_lower_bounds =\n',vstring2);
-    vector_write(fidi,sprintf('\t    '),vlower ,6,76);
-end
-if ~isempty(vupper)
-    fprintf(fidi,'\t  %s_upper_bounds =\n',vstring2);
-    vector_write(fidi,sprintf('\t    '),vupper ,6,76);
-end
-if ~isempty(vmean)
-    fprintf(fidi,'\t  %s_means =\n',vstring2);
-    vector_write(fidi,sprintf('\t    '),vmean  ,6,76);
-end
-if ~isempty(vstddev)
-    fprintf(fidi,'\t  %s_std_deviations =\n',vstring2);
-    vector_write(fidi,sprintf('\t    '),vstddev,6,76);
-end
-if ~isempty(vinitst)
-    fprintf(fidi,'\t  %s_initial_state =\n',vstring2);
-    vector_write(fidi,sprintf('\t    '),vinitst,6,76);
-end
-if ~isempty(vdesc)
-    fprintf(fidi,'\t  %s_descriptors =\n',vstring2);
-    vector_write(fidi,sprintf('\t    '),vdesc  ,6,76);
+fprintf(fidi,'\t%s = %d\n',cstring,length(dvar));
+if ~isempty(pinitpt)
+    fprintf(fidi,'\t  %s_initial_point =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pinitpt,6,76);
+end
+if ~isempty(plower)
+    fprintf(fidi,'\t  %s_lower_bounds =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),plower ,6,76);
+end
+if ~isempty(pupper)
+    fprintf(fidi,'\t  %s_upper_bounds =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pupper ,6,76);
+end
+if ~isempty(pmean)
+    fprintf(fidi,'\t  %s_means =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pmean  ,6,76);
+end
+if ~isempty(pstddev)
+    fprintf(fidi,'\t  %s_std_deviations =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pstddev,6,76);
+end
+if ~isempty(pinitst)
+    fprintf(fidi,'\t  %s_initial_state =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pinitst,6,76);
+end
+if ~isempty(pstype)
+    fprintf(fidi,'\t  %s_scale_types =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pstype ,6,76);
+end
+if ~isempty(pscale)
+    fprintf(fidi,'\t  %s_scales =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pscale ,6,76);
+end
+if ~isempty(pdesc)
+    fprintf(fidi,'\t  %s_descriptors =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pdesc  ,6,76);
+end
+
+end
+
+%%  function to write linear constraint sets
+
+function []=lcsets_write(fidi,dvar,lcspec)
+
+for i=1:length(lcspec)
+    switch(lcspec{i})
+        case 'lic'
+            cstring ='linear_inequality_constraints';
+            cstring2='linear_inequality';
+        case 'lec'
+            cstring ='linear_equality_constraints';
+            cstring2='linear_equality';
+        otherwise
+            warning('lcspec_write:unrec_lcspec',...
+                'Unrecognized linear constraint ''%s''.',lcspec{i});
+    end
+    
+    if isfield(dvar,lcspec{i})
+        lclist_write(fidi,cstring,cstring2,dvar.(lcspec{i}));
+    end
+end
+
+end
+
+%%  function to write linear constraint list
+
+function []=lclist_write(fidi,cstring,cstring2,dvar)
+
+%  put linear constraints into lists for writing
+
+pmatrix=prop_matrix(dvar);
+plower =prop_lower (dvar);
+pupper =prop_upper (dvar);
+ptarget=prop_target(dvar);
+pstype =prop_stype (dvar);
+pscale =prop_scale (dvar);
+
+%  write linear constraints
+
+disp(sprintf('  Writing %d %s linear constraints.',...
+    length(dvar),class(dvar)));
+
+if ~isempty(pmatrix)
+    fprintf(fidi,'\t  %s_matrix =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pmatrix,6,76);
+end
+if ~isempty(plower)
+    fprintf(fidi,'\t  %s_lower_bounds =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),plower ,6,76);
+end
+if ~isempty(pupper)
+    fprintf(fidi,'\t  %s_upper_bounds =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pupper ,6,76);
+end
+if ~isempty(ptarget)
+    fprintf(fidi,'\t  %s_targets =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),ptarget,6,76);
+end
+if ~isempty(pstype)
+    fprintf(fidi,'\t  %s_scale_types =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pstype ,6,76);
+end
+if ~isempty(pscale)
+    fprintf(fidi,'\t  %s_scales =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pscale ,6,76);
 end
 
@@ -341 +358 @@
 %%  function to write the responses section of the file
 
-function []=responses_write(fidi,method,dresp,params)
+function []=responses_write(fidi,dmeth,dresp,params)
 
 display('Writing responses section of Dakota input file.');
@@ -349 +366 @@
 %  functions, gradients, and hessians vary by method
 
-switch method
-%     case {'dot_bfgs','dot_frcg','dot_mmfd','dot_slp','dot_sqp'}
-%     case {'npsol_sqp'}
-    case {'conmin_frcg','conmin_mfd'}
-        rsets_write(fidi,dresp,'of','nic','nec');
-        ghspec_write(fidi,params,'grad');
-%     case {'optpp_cg','optpp_q_newton','optpp_fd_newton',...
-%             'optpp_newton','optpp_pds'}
-%     case {'asynch_pattern_search'}
-%     case {'coliny_cobyla','coliny_direct','coliny_ea',...
-%             'coliny_pattern_search','coliny_solis_wets'}
-%     case {'ncsu_direct'}
-%     case {'moga','soga'}
-%     case {'nl2sol','nlssol_sqp','optpp_g_newton'}
-    case {'nond_sampling'}
-        rsets_write(fidi,dresp,'rf');
-        ghspec_write(fidi,params);
-    case {'nond_local_reliability'}
-        rsets_write(fidi,dresp,'rf');
-        ghspec_write(fidi,params,'grad');
-%     case {'dace','fsu_quasi_mc','fsu_cvt'}
-%     case {'vector_parameter_study','list_parameter_study',...
-%             'centered parameter_study','multidim_parameter_study'}
-end
+rsets_write(fidi,dresp,dmeth.responses);
+ghspec_write(fidi,params,dmeth.ghspec);
 
 fprintf(fidi,'\n');
@@ -380 +375 @@
 %%  function to write response sets
 
-function []=rsets_write(fidi,dresp,varargin)
+function []=rsets_write(fidi,dresp,responses)
 
 rdesc={};
 
-for i=1:length(varargin)
-    switch(varargin{i})
+for i=1:length(responses)
+    switch(responses{i})
         case 'of'
-            rstring ='objective_functions';
-            rstring2='objective_function';
+            cstring ='objective_functions';
+            cstring2='objective_function';
         case 'lst'
-            rstring ='least_squares_terms';
-            rstring2='least_squares_term';
+            cstring ='least_squares_terms';
+            cstring2='least_squares_term';
         case 'nic'
-            rstring ='nonlinear_inequality_constraints';
-            rstring2='nonlinear_inequality';
+            cstring ='nonlinear_inequality_constraints';
+            cstring2='nonlinear_inequality';
         case 'nec'
-            rstring ='nonlinear_equality_constraints';
-            rstring2='nonlinear_equality';
+            cstring ='nonlinear_equality_constraints';
+            cstring2='nonlinear_equality';
         case 'rf'
-            rstring ='response_functions';
-            rstring2='response_function';
+            cstring ='response_functions';
+            cstring2='response_function';
         otherwise
             warning('rsets_write:unrec_resp',...
-                'Unrecognized response ''%s''.',varargin{i});
+                'Unrecognized response ''%s''.',responses{i});
     end
     
-    if isfield(dresp,varargin{i})
-        [rdesc]=rlist_write(fidi,rstring,rstring2,dresp.(varargin{i}),rdesc);
+    if isfield(dresp,responses{i})
+        [rdesc]=rlist_write(fidi,cstring,cstring2,dresp.(responses{i}),rdesc);
     end
 end
@@ -422 +417 @@
 %%  function to write response list
 
-function [rdesc]=rlist_write(fidi,rstring,rstring2,dresp,rdesc)
+function [rdesc]=rlist_write(fidi,cstring,cstring2,dresp,rdesc)
 
 %  put responses into lists for writing
 
-rstype =dresp_stype (dresp);
-rscale =dresp_scale (dresp);
-rweight=dresp_weight(dresp);
-rlower =dresp_lower (dresp);
-rupper =dresp_upper (dresp);
-rtarget=dresp_target(dresp);
+pstype =prop_stype (dresp);
+pscale =prop_scale (dresp);
+pweight=prop_weight(dresp);
+plower =prop_lower (dresp);
+pupper =prop_upper (dresp);
+ptarget=prop_target(dresp);
 
 %  accumulate descriptors into list for subsequent writing
 
-rdesc(end+1:end+numel(dresp))=dresp_desc  (dresp);
+rdesc(end+1:end+numel(dresp))=prop_desc  (dresp);
 
 %  write responses
@@ -441 +436 @@
 disp(sprintf('  Writing %d %s responses.',length(dresp),class(dresp)));
 
-fprintf(fidi,'\tnum_%s = %d\n',rstring,length(dresp));
-if ~isempty(rstype)
-    fprintf(fidi,'\t  %s_scale_types =\n',rstring2);
-    vector_write(fidi,sprintf('\t    '),rstype ,6,76);
-end
-if ~isempty(rscale)
-    fprintf(fidi,'\t  %s_scales =\n',rstring2);
-    vector_write(fidi,sprintf('\t    '),rscale ,6,76);
-end
-if ~isempty(rweight)
-    switch rstring2
+fprintf(fidi,'\tnum_%s = %d\n',cstring,length(dresp));
+if ~isempty(pstype)
+    fprintf(fidi,'\t  %s_scale_types =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pstype ,6,76);
+end
+if ~isempty(pscale)
+    fprintf(fidi,'\t  %s_scales =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pscale ,6,76);
+end
+if ~isempty(pweight)
+    switch cstring2
         case 'objective_function'
-            wtype='multi_objective';
+            fprintf(fidi,'\t  %s_weights =\n','multi_objective');
+            vector_write(fidi,sprintf('\t    '),pweight,6,76);
         case 'least_squares_term'
-            wtype='least_squares';
-    end
-    fprintf(fidi,'\t  %s_weights =\n',wtype);
-    vector_write(fidi,sprintf('\t    '),rweight,6,76);
-end
-if ~isempty(rlower)
-    fprintf(fidi,'\t  %s_lower_bounds =\n',rstring2);
-    vector_write(fidi,sprintf('\t    '),rlower ,6,76);
-end
-if ~isempty(rupper)
-    fprintf(fidi,'\t  %s_upper_bounds =\n',rstring2);
-    vector_write(fidi,sprintf('\t    '),rupper ,6,76);
-end
-if ~isempty(rtarget)
-    fprintf(fidi,'\t  %s_targets =\n',rstring2);
-    vector_write(fidi,sprintf('\t    '),rtarget,6,76);
+            fprintf(fidi,'\t  %s_weights =\n','least_squares');
+            vector_write(fidi,sprintf('\t    '),pweight,6,76);
+    end
+end
+if ~isempty(plower)
+    fprintf(fidi,'\t  %s_lower_bounds =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),plower ,6,76);
+end
+if ~isempty(pupper)
+    fprintf(fidi,'\t  %s_upper_bounds =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),pupper ,6,76);
+end
+if ~isempty(ptarget)
+    fprintf(fidi,'\t  %s_targets =\n',cstring2);
+    vector_write(fidi,sprintf('\t    '),ptarget,6,76);
 end
 
@@ -477 +472 @@
 %%  function to write gradient and hessian specifications
 
-function []=ghspec_write(fidi,params,varargin)
+function []=ghspec_write(fidi,params,ghspec)
 
 %  gradients
 
-if findstring(varargin,'grad')
+if find_string(ghspec,'grad')
     if     params.numerical_gradients
         param_write(fidi,'\t','numerical_gradients','','\n',params);
@@ -497 +492 @@
 %  hessians (no implemented methods use them yet)
 
-if findstring(varargin,'hess')
+if find_string(ghspec,'hess')
     error('Hessians needed by method but not provided.');
 else
     fprintf(fidi,'\tno_hessians\n');
-end
-
-end
-
-%%  function to find a string in a cell array
-
-function [ifound]=findstring(cells,str)
-
-ifound=false;
-
-for i=1:length(cells)
-    if ischar(cells{i}) && strcmp(cells{i},str)
-        ifound=i;
-        return
-    end
 end
 
@@ -543 +523 @@
 
 end
+
 %%  function to write a vector on multiple lines
 
@@ -560 +541 @@
 lsvec=cmax;
 
+%  transpose vector from column-wise to row-wise
+
+vec=vec';
+
 %  assemble each line, flushing when necessary
 
-for i=1:length(vec)
+for i=1:numel(vec)
     if isnumeric(vec(i))
         sitem=sprintf('%g'    ,vec(i));

issm/trunk/src/m/solutions/dakota/dakota_m_write.m

@@ -3 +3 @@
 %  driver.
 %
-%  []=dakota_m_write(md,method,dvar,dresp,filem,params,package)
+%  []=dakota_m_write(method,dmeth,dvar,dresp,params,filem,package)
 %
-function []=dakota_m_write(md,method,dvar,dresp,filem,params,package)
+function []=dakota_m_write(method,dmeth,dvar,dresp,params,filem,package,varargin)
 
 if ~nargin
@@ -14 +14 @@
 %  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 ~isfield(params,'npart')
-    params.npart=10;
 end
 
@@ -42 +47 @@
 %  write variables into the Matlab m-file
 
-variables_write(md,fidm,method,dvar);
+variables_write(fidm,dmeth,dvar,params,varargin{:});
 
 %  write solution into the Matlab m-file
@@ -50 +55 @@
 %  write responses into the Matlab m-file
 
-responses_write(fidm,method,params,dresp);
+responses_write(fidm,dmeth,params,dresp);
 
 %  write end of the Matlab m-file
@@ -84 +89 @@
 end
 fprintf(fidm,'\tloadmodel(infile);\n\n');
-fprintf(fidm,'\tmd=qmuname(md);\n\n');
+% fprintf(fidm,'\tmd=qmuname(md);\n\n');
 
 end
@@ -90 +95 @@
 %%  function to write design variables into the Matlab m-file
 
-function []=variables_write(md,fidm,method,dvar)
+function []=variables_write(fidm,dmeth,dvar,params,varargin)
 
 display('Writing variables for Matlab m-file.');
 
-fprintf(fidm,'%%  Apply the design variables.\n\n');
+fprintf(fidm,'%%  Apply the variables.\n\n');
 ixc=0;
 
 %  variables vary by method
 
-switch method
-%     case {'dot_bfgs','dot_frcg','dot_mmfd','dot_slp','dot_sqp'}
-%     case {'npsol_sqp'}
-    case {'conmin_frcg','conmin_mfd'}
-        ixc=vsets_write(md,fidm,ixc,dvar,'cdv','csv');
-%     case {'optpp_cg','optpp_q_newton','optpp_fd_newton',...
-%             'optpp_newton','optpp_pds'}
-%     case {'asynch_pattern_search'}
-%     case {'coliny_cobyla','coliny_direct','coliny_ea',...
-%             'coliny_pattern_search','coliny_solis_wets'}
-%     case {'ncsu_direct'}
-%     case {'moga','soga'}
-%     case {'nl2sol','nlssol_sqp','optpp_g_newton'}
-    case {'nond_sampling'}
-        ixc=vsets_write(md,fidm,ixc,dvar,'nuv','csv');
-    case {'nond_local_reliability'}
-        ixc=vsets_write(md,fidm,ixc,dvar,'nuv','csv');
-%     case {'dace','fsu_quasi_mc','fsu_cvt'}
-%     case {'vector_parameter_study','list_parameter_study',...
-%             'centered parameter_study','multidim_parameter_study'}
-end
+ixc=vsets_write(fidm,ixc,dvar,dmeth.variables,params,varargin{:});
 
 end
@@ -125 +110 @@
 %%  function to write variable sets into the Matlab m-file
 
-function [ixc]=vsets_write(md,fidm,ixc,dvar,varargin)
-
-for i=1:length(varargin)
-    if isfield(dvar,varargin{i})
-        ixc=vlist_write(md,fidm,ixc,varargin{i},dvar.(varargin{i}));
+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
@@ -137 +122 @@
 %%  function to write variable list into the Matlab m-file
 
-function [ixc]=vlist_write(md,fidm,ixc,vtype,dvar)
+function [ixc]=vlist_write(fidm,ixc,vtype,dvar,params,varargin)
 
 disp(sprintf('  Writing %d %s variables.',length(dvar),class(dvar)));
@@ -156 +141 @@
                 %now, we need a string to put in the matlab file, which will update all the design variables 
                 %for  this descriptor.
-                [string,ixc]=QmuUpdateFunctions(md,ixc,descriptor,dvar,i);
+                [string,ixc]=QmuUpdateFunctions(ixc,descriptor,dvar,params,i,varargin{:});
 
                 %dump this string in the matlab file.
@@ -178 +163 @@
 %%  function to write responses into the Matlab m-file
 
-function []=responses_write(fidm,method,params,dresp)
+function []=responses_write(fidm,dmeth,params,dresp)
 
 display('Writing responses for Matlab m-file.');
 
-fprintf(fidm,'%%  Calculate the response functions.\n\n');
+fprintf(fidm,'%%  Calculate the responses.\n\n');
 ifnvals=0;
 
@@ -191 +176 @@
 %  responses vary by method
 
-switch method
-%     case {'dot_bfgs','dot_frcg','dot_mmfd','dot_slp','dot_sqp'}
-%     case {'npsol_sqp'}
-    case {'conmin_frcg','conmin_mfd'}
-        ifnvals=rsets_write(fidm,ifnvals,params,dresp,'of','nic','nec');
-%     case {'optpp_cg','optpp_q_newton','optpp_fd_newton',...
-%             'optpp_newton','optpp_pds'}
-%     case {'asynch_pattern_search'}
-%     case {'coliny_cobyla','coliny_direct','coliny_ea',...
-%             'coliny_pattern_search','coliny_solis_wets'}
-%     case {'ncsu_direct'}
-%     case {'moga','soga'}
-%     case {'nl2sol','nlssol_sqp','optpp_g_newton'}
-    case {'nond_sampling'}
-        ifnvals=rsets_write(fidm,ifnvals,params,dresp,'rf');
-    case {'nond_local_reliability'}
-        ifnvals=rsets_write(fidm,ifnvals,params,dresp,'rf');
-%     case {'dace','fsu_quasi_mc','fsu_cvt'}
-%     case {'vector_parameter_study','list_parameter_study',...
-%             'centered parameter_study','multidim_parameter_study'}
-end
+ifnvals=rsets_write(fidm,ifnvals,params,dresp,dmeth.responses);
 
 fprintf(fidm,'\n');
@@ -222 +187 @@
 %%  function to write response sets into the Matlab m-file
 
-function [ifnvals]=rsets_write(fidm,ifnvals,params,dresp,varargin)
-
-for i=1:length(varargin)
-    if isfield(dresp,varargin{i})
-        ifnvals=rlist_write(fidm,ifnvals,params,varargin{i},dresp.(varargin{i}));
+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

issm/trunk/src/m/solutions/dakota/qmu.m

@@ -2 +2 @@
 %INPUT function md=qmu(md,package)
 %Deal with coupled ISSM or Cielo/ Dakota runs, to do sensitivity analyses.
+
+global ISSM_DIR;
 
 %first create temporary directory in which we will work
@@ -7 +9 @@
 qmudir='qmu';
 if exist(qmudir,'dir')
-    overwrite=input('Overwrite existing ''qmu'' directory? Y/N [N]: ', 's');
+    overwrite=input(['Overwrite existing ''' qmudir ''' directory? Y/N [N]: '], 's');
     if strncmpi(overwrite,'y',1)
         system(['rm -rf ' qmudir]);
     else
-        error('Existing ''qmu'' directory not overwritten');
+        error('Existing ''%s'' directory not overwritten.',qmudir);
     end
 end
@@ -22 +24 @@
 
 %create m and in files for dakota
-dakota_in_data(md,'qmu',package);
+if ~isfield(md.qmu_params,'analysis_driver') || ...
+    isempty(md.qmu_params.analysis_driver)
+    md.qmu_params.analysis_driver=[ISSM_DIR '/src/m/solutions/dakota/cielo_ice_script.sh'];
+end
+
+if (numel(md.qmu_method) > 1)
+    imeth=input('Which method? ');
+else
+    imeth=1;
+end
+dakota_in_data(md.qmu_method(imeth),md.variables,md.responses,md.qmu_params,'qmu',package,md);
+cd ..
+return
 
 %call dakota
 system('dakota -i qmu.in -o qmu.out -e qmu.err');
+% system('export MPIRUN_NPROCS=8;mpirun -np 4 dakota -i qmu.in -o qmu.out -e qmu.err');
 
 %parse inputs and results from dakota

issm/trunk/src/m/solutions/dakota/setupdesign/QmuSetupDesign.m

@@ -1 @@
-function dvar=QmuSetupDesign(md,dvar,variables)
+function dvar=QmuSetupDesign(dvar,variables,params,varargin)
 
 %get descriptor
@@ -7 +7 @@
 if strcmpi(descriptor,'rho_ice')
 
-        dvar=setuprhoice(md,dvar,variables);
+        dvar=setuprhoice(dvar,variables,params,varargin{:});
+
+elseif strcmpi(descriptor,'rho_water')
+
+        dvar=setuprhowater(dvar,variables,params,varargin{:});
+
+elseif strcmpi(descriptor,'heat_capacity')
+
+        dvar=setupheatcapacity(dvar,variables,params,varargin{:});
+
+elseif strcmpi(descriptor,'thermal_conductivity')
+
+        dvar=setupthermalconductivity(dvar,variables,params,varargin{:});
+
+elseif strcmpi(descriptor,'gravity')
+
+        dvar=setupgravity(dvar,variables,params,varargin{:});
 
 elseif strcmpi(descriptor,'thickness')
 
-        dvar=setupthickness(md,dvar,variables);
+        dvar=setupthickness(dvar,variables,params,varargin{:});
 
 elseif strcmpi(descriptor,'drag')
 
-        dvar=setupdrag(md,dvar,variables);
+        dvar=setupdrag(dvar,variables,params,varargin{:});
 
 elseif strcmpi(descriptor,'riftsfriction')
 
-        dvar=setupriftsfriction(md,dvar,variables);
+        dvar=setupriftsfriction(dvar,variables,params,varargin{:});
 
 else

issm/trunk/src/m/solutions/dakota/setupdesign/setupdrag.m

@@ -1 @@
-function dvar=setupdrag(md,dvar,variables)
+function dvar=setupdrag(dvar,variables,params,varargin)
+
+for i=1:length(varargin)
+    if strcmp(class(varargin{i}),'model')
+        md=varargin{i};
+        break;
+    end
+end
 
 %ok, dealing with semi-discrete distributed variable. Distribute according to how many 

issm/trunk/src/m/solutions/dakota/setupdesign/setuprhoice.m

@@ -1 @@
-function dvar=setuprhoice(md,dvar,variables)
+function dvar=setuprhoice(dvar,variables,params,varargin)
 
 dvar(end+1)=variables;

issm/trunk/src/m/solutions/dakota/setupdesign/setupriftsfriction.m

@@ -1 @@
-function dvar=setupriftsfriction(md,dvar,variables)
+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.

issm/trunk/src/m/solutions/dakota/setupdesign/setupthickness.m

@@ -1 @@
-function dvar=setupthickness(md,dvar,variables)
+function dvar=setupthickness(dvar,variables,params,varargin)
+
+for i=1:length(varargin)
+    if strcmp(class(varargin{i}),'model')
+        md=varargin{i};
+        break;
+    end
+end
 
 %ok, dealing with semi-discrete distributed variable. Distribute according to how many 

issm/trunk/src/m/solutions/dakota/updatefunctions/QmuUpdateFunctions.m

@@ -1 @@
-function [string,ixc]=QmuUpdateFunctions(md,ixc,descriptor,dvar,i)
+function [string,ixc]=QmuUpdateFunctions(ixc,descriptor,dvar,params,i,varargin)
 
 if strcmpi(descriptor,'rho_ice')
-        [string,ixc]=updaterho_ice(md,ixc,dvar,i);
+        [string,ixc]=updaterho_ice(ixc,dvar,params,i,varargin{:});
 
 elseif strcmpi(descriptor,'rho_water')
-        [string,ixc]=updaterho_water(md,ixc,dvar,i);
+        [string,ixc]=updaterho_water(ixc,dvar,params,i,varargin{:});
 
 elseif strcmpi(descriptor,'heatcapacity')
-        [string,ixc]=updateheatcapacity(md,ixc,dvar,i);
+        [string,ixc]=updateheatcapacity(ixc,dvar,params,i,varargin{:});
 
 elseif strcmpi(descriptor,'thermalconductivity')
-        [string,ixc]=updatethermalconductivity(md,ixc,dvar,i);
+        [string,ixc]=updatethermalconductivity(ixc,dvar,params,i,varargin{:});
 
 elseif strcmpi(descriptor,'gravity')
-        [string,ixc]=updategravity(md,ixc,dvar,i);
+        [string,ixc]=updategravity(ixc,dvar,params,i,varargin{:});
 
 elseif strcmpi(descriptor,'thickness')
-        [string,ixc]=updatethickness(md,ixc,dvar);
+        [string,ixc]=updatethickness(ixc,dvar,params,varargin{:});
 
 elseif strcmpi(descriptor,'drag')
-        [string,ixc]=updatedrag(md,ixc,dvar);
+        [string,ixc]=updatedrag(ixc,dvar,params,varargin{:});
 
 elseif strcmpi(descriptor,'riftsfriction')
-        [string,ixc]=updateriftsfriction(md,ixc,dvar);
+        [string,ixc]=updateriftsfriction(ixc,dvar,params,varargin{:});
 
 else
-        warning('QmuUpdateFunctions warning message: ','Unrecognized design variable: %s',descriptor)
+        warning('QmuUpdateFunctions:warning','Unrecognized design variable: %s',descriptor)
 end

issm/trunk/src/m/solutions/dakota/updatefunctions/updatedrag.m

@@ -1 @@
-function [string,ixc]=updatedrag(md,ixc,dvar)
+function [string,ixc]=updatedrag(ixc,dvar,params,varargin)
+
+for i=1:length(varargin)
+    if strcmp(class(varargin{i}),'model')
+        md=varargin{i};
+        break;
+    end
+end
 
 string=sprintf('\n\t[epart,npart]=MeshPartition(md,%d);\n\n',md.npart);
@@ -7 +14 @@
         ixc=ixc+1;
                 idrag =sscanf(dvar(j).descriptor(5:end),'%d');
-        if strcmpi(md.analysis_driver,'matlab')
+        if strcmpi(params.analysis_driver,'matlab')
             string=[string sprintf('\tdragi(%d,1)=Dakota.xC(%d);\n',idrag,ixc)];
         else

issm/trunk/src/m/solutions/dakota/updatefunctions/updategravity.m

@@ -1 @@
-function [string,ixc]=updategravity(md,ixc,dvar,i)
+function [string,ixc]=updategravity(ixc,dvar,params,i,varargin)
         
 ixc=ixc+1;
-if strcmpi(md.analysis_driver,'matlab')
+if strcmpi(params.analysis_driver,'matlab')
     string=sprintf('\tmd.g=Dakota.xC(%d);\n',ixc);
 else

issm/trunk/src/m/solutions/dakota/updatefunctions/updateheatcapacity.m

@@ -1 @@
-function [string,ixc]=updateheatcapacity(md,ixc,dvar,i)
+function [string,ixc]=updateheatcapacity(ixc,dvar,params,i,varargin)
 
 ixc=ixc+1;
-if strcmpi(md.analysis_driver,'matlab')
+if strcmpi(params.analysis_driver,'matlab')
     string=sprintf('\tmd.heatcapacity=Dakota.xC(%d);\n',ixc);
 else

issm/trunk/src/m/solutions/dakota/updatefunctions/updaterho_ice.m

@@ -1 @@
-function [string,ixc]=updaterho_ice(md,ixc,dvar,i)
+function [string,ixc]=updaterho_ice(ixc,dvar,params,i,varargin)
 
 ixc=ixc+1;
-if strcmpi(md.analysis_driver,'matlab')
+if strcmpi(params.analysis_driver,'matlab')
     string=sprintf('\tmd.rho_ice=Dakota.xC(%d);\n',ixc);
 else

issm/trunk/src/m/solutions/dakota/updatefunctions/updaterho_water.m

@@ -1 @@
-function [string,ixc]=updaterho_water(md,ixc,dvar,i)
+function [string,ixc]=updaterho_water(ixc,dvar,params,i,varargin)
 
 ixc=ixc+1;
-if strcmpi(md.analysis_driver,'matlab')
+if strcmpi(params.analysis_driver,'matlab')
     string=sprintf('\tmd.rho_water=Dakota.xC(%d);\n',ixc);
 else

issm/trunk/src/m/solutions/dakota/updatefunctions/updateriftsfriction.m

@@ -1 @@
-function [string,ixc]=updateriftsfriction(md,ixc,dvar)
+function [string,ixc]=updateriftsfriction(ixc,dvar,params,varargin)
+
+for i=1:length(varargin)
+    if strcmp(class(varargin{i}),'model')
+        md=varargin{i};
+        break;
+    end
+end
 
 string=sprintf('\n');
@@ -6 +13 @@
         ixc=ixc+1;
                 irift=sscanf(dvar(j).descriptor(14:end),'%d');
-        if strcmpi(md.analysis_driver,'matlab')
+        if strcmpi(params.analysis_driver,'matlab')
             string=[string sprintf('\tmd.rifts(%d).friction=Dakota.xC(%d);\n',irift,ixc)];
         else

issm/trunk/src/m/solutions/dakota/updatefunctions/updatethermalconductivity.m

@@ -1 @@
-function [string,ixc]=updatethermalconductivity(md,ixc,dvar,i)
+function [string,ixc]=updatethermalconductivity(ixc,dvar,params,i,varargin)
 
 ixc=ixc+1;
-if strcmpi(md.analysis_driver,'matlab')
+if strcmpi(params.analysis_driver,'matlab')
     string=sprintf('\tmd.thermalconductivity=Dakota.xC(%d);\n',ixc);
 else

issm/trunk/src/m/solutions/dakota/updatefunctions/updatethickness.m

@@ -1 @@
-function [string,ixc]=updatethickness(md,ixc,dvar)
+function [string,ixc]=updatethickness(ixc,dvar,params,varargin)
+
+for i=1:length(varargin)
+    if strcmp(class(varargin{i}),'model')
+        md=varargin{i};
+        break;
+    end
+end
 
 string=sprintf('\n\t[epart,npart]=MeshPartition(md,%d);\n',md.npart);
@@ -7 +14 @@
         ixc=ixc+1;
                 ithick=sscanf(dvar(j).descriptor(10:end),'%d');
-        if strcmpi(md.analysis_driver,'matlab')
+        if strcmpi(params.analysis_driver,'matlab')
             string=[string sprintf('\tthicki(%d,1)=Dakota.xC(%d);\n',ithick,ixc)];
         else

issm/trunk/src/m/solutions/ice/ModelProcessorDiagnosticHoriz.m

@@ -312 +312 @@
 if ~isempty(md.penalties) & ~isnan(md.penalties),
         for i=1:size(md.penalties,1),
-                for j=1:(md.numlayers-1),
+                for j=2:size(md.penalties,2),
 
                         %constrain first dof
                         constraints(count).constraint=rgb;
                         constraints(count).constraint.grid1=md.penalties(i,1);
-                        constraints(count).constraint.grid2=md.penalties(i,j+1);
+                        constraints(count).constraint.grid2=md.penalties(i,j);
                         constraints(count).constraint.dof=1;
                         count=count+1;
@@ -324 +324 @@
                         constraints(count).constraint=rgb;
                         constraints(count).constraint.grid1=md.penalties(i,1);
-                        constraints(count).constraint.grid2=md.penalties(i,j+1);
+                        constraints(count).constraint.grid2=md.penalties(i,j);
                         constraints(count).constraint.dof=2;
                         count=count+1;

issm/trunk/src/m/solutions/ice/ModelProcessorPrognostic.m

@@ -116 +116 @@
 loads=struct('load',cell([0,1]));
 
-%Single point constraints: none;
-constraints=struct('constraint',cell(0,0));
+
+%Initialize constraints structure
+constraints=struct('constraint',cell(length(find(md.gridondirichlet_prog)),1));
+
+pos=find(md.gridondirichlet_prog);
+count=1;
+for i=1:length(find(md.gridondirichlet_prog)),
+
+                %constrain the thickness
+                constraints(count).constraint=spc;
+                constraints(count).constraint.grid=pos(i);
+                constraints(count).constraint.dof=1;
+                constraints(count).constraint.value=md.dirichletvalues_prog(pos(i)); %in m
+                count=count+1;
+end
 
 %Last thing, return a partitioning vector, and its transpose.

issm/trunk/src/m/solutions/ice/StrainRateCompute.m

@@ -34 +34 @@
 
                         strainratematrix=[strainratevector(1) strainratevector(3) 
-                                      strainratevector(3)  strainratevector(2)];
+                                                                        strainratevector(3)  strainratevector(2)];
 
                         %eigen values and vectors
@@ -41 +41 @@
                         %Plug into global vectors
                         strainrate1(n,:)=strainratevector;
-                        A1(n,1)=value(1,1); A2(n,1)=value(2,2);
-                        Vx1(n,1)=directions(1,1); Vx2(n,1)=directions(1,2);
-                        Vy1(n,1)=directions(2,1); Vy2(n,1)=directions(2,2);
+                        A1(n,1)=value(1,1); A2(n,1)=value(2,2);
+                        Vx1(n,1)=directions(1,1); Vx2(n,1)=directions(1,2);
+                        Vy1(n,1)=directions(2,1); Vy2(n,1)=directions(2,2);
                 end
         end
@@ -54 +54 @@
         strainrate.xy=strainrate1(:,3)*(365.25*24*3600);
         %principal axis and values
-        strainrate.principalvalue2=A1;
+        strainrate.principalvalue2=A1*(365.25*24*3600);  %strain rate in 1/a instead of 1/s
         strainrate.principalaxis2=[Vx1 Vy1];
-        strainrate.principalvalue1=A2;
+        strainrate.principalvalue1=A2*(365.25*24*3600);
         strainrate.principalaxis1=[Vx2 Vy2];
         %Norm or effective value

issm/trunk/src/m/solutions/macayeal/macayealcontrol.m

@@ -142 +142 @@
           index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,B_bar);
 
-      if niteration==1
-          scrsz = get(0,'ScreenSize');
-          figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3])
-      end
-        plotmodel(md,'data',sqrt(u.^2+v.^2)*yts,'title','Modeled velocity',...
-                'data',sqrt(adjointu.^2+adjointv.^2)*yts,'title','Adjoint vectors','colorbar#all', 'on',...
-                'data',sqrt(vx_obs.^2+vy_obs.^2)*yts,'title','Observed velocity',...
-                'data',100*abs(sqrt(vx_obs.^2+vy_obs.^2)-sqrt(u.^2+v.^2))./sqrt(vx_obs.^2+vy_obs.^2),'title','Relative misfit','caxis#4',[0 100],'figure',1);pause(1);
+                 if md.plot==1,
+                         if niteration==1
+                                 scrsz = get(0,'ScreenSize');
+                                 figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3])
+                         end
+                         plotmodel(md,'data',sqrt(u.^2+v.^2)*yts,'title','Modeled velocity',...
+                                 'data',sqrt(adjointu.^2+adjointv.^2)*yts,'title','Adjoint vectors','colorbar#all', 'on',...
+                                 'data',sqrt(vx_obs.^2+vy_obs.^2)*yts,'title','Observed velocity',...
+                                 'data',100*abs(sqrt(vx_obs.^2+vy_obs.^2)-sqrt(u.^2+v.^2))./sqrt(vx_obs.^2+vy_obs.^2),'title','Relative misfit','caxis#4',[0 100],'figure',1);drawnow;
+                 end
 
        %Keep track of u and v to use in objectivefunction_C: 
@@ -164 +166 @@
        old_direction=direction;
 
-       %visualize direction 
-       if niteration==1
-           scrsz = get(0,'ScreenSize');
-           figure('Position',[10 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3])
-       end              
-       plotmodel(md,'data',direction,'title','Orthogonal direction','figure',2);pause(1);
+                 %visualize direction 
+                 if md.plot==1,
+                         if niteration==1
+                                 scrsz = get(0,'ScreenSize');
+                                 figure('Position',[10 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3])
+                         end              
+                         plotmodel(md,'data',direction,'title','Orthogonal direction','figure',2);drawnow;
+                 end
 
 
@@ -208 +212 @@
 
        %visualize new distribution. 
-       if niteration==1
-           scrsz = get(0,'ScreenSize');
-           figure('Position',[scrsz(3)*1.97/3 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3])
-       end
-       plotmodel(md,'data',B*B_to_p,'title',['B at iteration ' num2str(niteration)],'caxis',[200 900],'colorbar','on','figure',3);pause(1);
+                 if md.plot==1,
+                         if niteration==1
+                                 scrsz = get(0,'ScreenSize');
+                                 figure('Position',[scrsz(3)*1.97/3 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3])
+                         end
+                         plotmodel(md,'data',B*B_to_p,'title',['B at iteration ' num2str(niteration)],'caxis',[200 900],'colorbar','on','figure',3);drawnow;
+                 end
 
        %evaluate new misfit. 
@@ -243 +249 @@
             index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,drag_coeff,drag_coeff_bar);
 
-      if niteration==1
-          scrsz = get(0,'ScreenSize');
-          figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3])
-      end
-        plotmodel(md,'data',sqrt(u.^2+v.^2),'title','Modeled velocity',...
-                'data',sqrt(adjointu.^2+adjointv.^2),'title','Adjoint vectors','colorbar', 'all',...
-                'data',sqrt(vx_obs.^2+vy_obs.^2),'title','Observed velocity',...
-                'data',100*abs(sqrt(vx_obs.^2+vy_obs.^2)-sqrt(u.^2+v.^2))./sqrt(vx_obs.^2+vy_obs.^2),'title','Relative misfit','caxis#3',[0 100],'figure',1);pause(1);
+                  if md.plot==1,
+                          if niteration==1
+                                  scrsz = get(0,'ScreenSize');
+                                  figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3])
+                          end
+                          plotmodel(md,'data',sqrt(u.^2+v.^2),'title','Modeled velocity',...
+                                  'data',sqrt(adjointu.^2+adjointv.^2),'title','Adjoint vectors','colorbar', 'all',...
+                                  'data',sqrt(vx_obs.^2+vy_obs.^2),'title','Observed velocity',...
+                                  'data',100*abs(sqrt(vx_obs.^2+vy_obs.^2)-sqrt(u.^2+v.^2))./sqrt(vx_obs.^2+vy_obs.^2),'title','Relative misfit','caxis#3',[0 100],'figure',1);drawnow;
+                  end
 
         %Keep track of u and v to use in objectivefunction_C:
@@ -268 +276 @@
 
        %visualize direction 
-       if niteration==1
-           scrsz = get(0,'ScreenSize');
-           figure('Position',[10 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3])
-       end 
-       plotmodel(md,'data',direction,'title','Orthogonal direction','figure',2);pause(1);
+                 if md.plot==1,
+                         if niteration==1
+                                 scrsz = get(0,'ScreenSize');
+                                 figure('Position',[10 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3])
+                         end 
+                         plotmodel(md,'data',direction,'title','Orthogonal direction','figure',2);drawnow;
+                 end
 
 
@@ -319 +329 @@
 
        %visualize new distribution. 
-       plotmodel(md,'data',drag_coeff,'title',['Drag at iteration ' num2str(niteration)],'caxis',[0 1000],'colorbar','on','figure',3);pause(1);
+                 if md.plot==1,
+                         plotmodel(md,'data',drag_coeff,'title',['Drag at iteration ' num2str(niteration)],'caxis',[0 1000],'colorbar','on','figure',3);drawnow;
+                 end
 
        %evaluate new misfit.