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 a misfit, we need the observed velocity, etc ...
md=process_qmu_response_data(md);