Index: /issm/trunk-jpl/src/m/classes/results.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/results.py	(revision 25725)
+++ /issm/trunk-jpl/src/m/classes/results.py	(revision 25726)
@@ -1,2 +1,6 @@
+from copy import deepcopy
+
+import numpy as np
+
 from fielddisplay import fielddisplay
 
@@ -7,68 +11,125 @@
     Usage:
         md.results = results()
+
+    TODO:
+    - Rework so that a solution of arbitrary length (normal or transient) can 
+    initialized from one call to the results class constructor.
+    """
+
+    def __init__(self): #{{{
+        pass
+    #}}}
+
+    def __repr__(self): #{{{
+        s = ''
+        for key, value in self.__dict__.items():
+            if isinstance(value, resultsdakota):
+                lengthvalue = 1
+            else:
+                lengthvalue = len(value)
+            s += '    {}: [1x{} struct]\n'.format(key, lengthvalue)
+
+        return s
+    #}}}
+
+    def setdefaultparameters(self): #{{{
+        #do nothing
+        return self
+    #}}}
+
+    def checkconsistency(self, md, solution, analyses): #{{{
+        return md
+    #}}}
+
+    def marshall(self, prefix, md, fid): #{{{
+        pass
+    #}}}
+#}}}
+
+class resultsdakota(object): #{{{
+    """RESULTSDAKOTA class definition - Used to store results from a run of 
+    Dakota.
+
+    Usage:
+        md.results.dakota = resultsdakota()
+
+    NOTE: Values of attributes can themselves be instances of solution class.
+    """
+
+    def __init__(self): #{{{
+        pass
+    #}}}
+
+    def __repr__(self): #{{{
+        s = ''
+        for key, value in self.__dict__.items():
+            s += '    {}: '.format(key)
+            if isinstance(value, list):
+                s += '[{} element list]'.format(len(value))
+            else:
+                s += '{}'.format(value)
+            s += '\n'
+        return s
+    #}}}
+
+    def __len__(self): #{{{
+        return len(self.__dict__.keys())
+    #}}}
+
+    def setdefaultparameters(self): #{{{
+        #do nothing
+        return self
+    #}}}
+
+    def checkconsistency(self, md, solution, analyses): #{{{
+        return md
+    #}}}
+
+    def marshall(self, prefix, md, fid): #{{{
+        pass
+    #}}}
+#}}}
+
+class solution(object): #{{{
+    """SOLUTION class definition - Value of an attribute (which should be 
+    a string representing a solution type) of an instance of results class
+
+    Elements of self.steps should be instances of solutionstep.
+
+    Usage:
+        setattr(md.results, 'SolutionType', solution())
+
+    NOTE:
+    - Under MATLAB, this is implemented as a structure array
+    - Only when this instance of solution represents a transient solution 
+    should self.steps have more than one element
     """
 
     def __init__(self, *args): #{{{
-        pass
-    #}}}
-
-    def __repr__(self): #{{{
-        s = ''
-        if 'step' in self.__dict__:
-            s += '{}\n'.format(fielddisplay(self, 'step', "step number"))
-        if 'time' in self.__dict__:
-            s += '{}\n'.format(fielddisplay(self, 'time', "time value"))
-        if 'SolutionType' in self.__dict__:
-            s += '{}\n'.format(fielddisplay(self, 'SolutionType', "solution type"))
-
-        for name in list(self.__dict__.keys()):
-            if name not in ['step', 'time', 'SolutionType', 'errlog', 'outlog']:
-                if isinstance(getattr(self, name), list):
-                    s += '{}\n'.format(fielddisplay(self, name, "model results list"))
-                elif isinstance(getattr(self, name), results):
-                    s += '{}\n'.format(fielddisplay(self, name, "model results case"))
-                else:
-                    s += '{}\n'.format(fielddisplay(self, name, ""))
-
-        if 'errlog' in self.__dict__:
-            s += '{}\n'.format(fielddisplay(self, 'errlog', "error log file"))
-        if 'outlog' in self.__dict__:
-            s += '{}\n'.format(fielddisplay(self, 'outlog', "output log file"))
-
-        return s
-    #}}}
-
-    def setdefaultparameters(self): #{{{
-        #do nothing
-        return self
-    #}}}
-
-    def checkconsistency(self, md, solution, analyses): #{{{
-        return md
-    #}}}
-
-    def marshall(self, prefix, md, fid): #{{{
-        pass
-    #}}}
-#}}}
-
-class resultselement(object): #{{{
-    """RESULTSELEMENT class definition - Value of attribute of results
-
-    Usage:
-        setattr(md.results, '<solution_type>', resultselement())
-    """
-
-    def __init__(self, *args): #{{{
-        self.elements = [result()]
-    #}}}
-
-    def __repr__(self): #{{{
-        s = '  1x{} struct array with fields:\n'.format(len(self.elements))
-        s += '\n'
-
-        if len(self.elements) >= 1:
-            for name in list(self.elements[0].__dict__.keys()):
-                s += '    {}\n'.format(name)
+        if len(args) == 1:
+            arg = args[0]
+            if isinstance(arg, list):
+                self.steps = arg
+            else:
+                raise Exception('solution class error: if initializing with an argument, that argument should be an empty list or a list of instances of solutionstep')
+        else:
+            self.steps = [solutionstep()]
+    #}}}
+
+    def __deepcopy__(self, memo): #{{{
+        return solution(deepcopy(self.steps, memo))
+    #}}}
+
+    def __repr__(self): #{{{
+        s = ''
+        numsteps = len(self.steps)
+        if numsteps == 1:
+            for key, value in self.steps[0].__dict__.items():
+                s += '    {}: {}\n'.format(key, value)
+        else:
+            s = '  1x{} struct array with fields:\n'.format(numsteps)
+            s += '\n'
+            for fieldname in self.steps[0].getfieldnames():
+                s += '    {}\n'.format(fieldname)
 
         return s
@@ -76,5 +137,12 @@
 
     def __len__(self): #{{{
-        return len(self.elements)
+        return len(self.steps)
+    #}}}
+
+    def __getattr__(self, key): #{{{
+        if len(self.steps) == 1:
+            return getattr(self.steps[0], key)
+        else:
+            raise Exception('<results>.<solution> error: Currently, can only get a field if we are not working with a transient solution.')
     #}}}
 
@@ -82,36 +150,29 @@
         while True:
             try:
-                return self.elements[index]
+                return self.steps[index]
             except:
-                self.elements.append(result())
-    #}}}
-
-    def __setitem__(self, index, value): #{{{
-        while True:
-            try:
-                self.elements[index] = value
-                break
-            except:
-                self.elements.append(result())
-    #}}}
-
-    def setdefaultparameters(self): #{{{
-        return self
-    #}}}
-
-    def checkconsistency(self, md, solution, analyses): #{{{
-        return md
-    #}}}
-
-    def marshall(self, prefix, md, fid): #{{{
-        pass
-    #}}}
-#}}}
-
-class result(object): #{{{
-    """RESULT class definition - Element of resultselement::elements
-
-    Usage:
-        resultselement.elements.append(result())
+                self.steps.append(solutionstep())
+        else:
+            raise Exception('<results>.<solution>: either request a specific result by index or make sure that there is only a single result for this solution (cannot be a transient solution)')
+    #}}}
+
+    def setdefaultparameters(self): #{{{
+        return self
+    #}}}
+
+    def checkconsistency(self, md, solution, analyses): #{{{
+        return md
+    #}}}
+
+    def marshall(self, prefix, md, fid): #{{{
+        pass
+    #}}}
+#}}}
+
+class solutionstep(object): #{{{
+    """SOLUTIONSTEP class definition - Single element of <solution>.steps
+
+    Usage:
+        <solution>.steps.append(solutionstep())
     """
 
@@ -120,14 +181,37 @@
     #}}}
 
-    def setdefaultparameters(self): #{{{
-        return self
-    #}}}
-
-    def checkconsistency(self, md, solution, analyses): #{{{
-        return md
-    #}}}
-
-    def marshall(self, prefix, md, fid): #{{{
-        pass
-    #}}}
-#}}}
+    def __repr__(self): #{{{
+        s = ''
+        width = self.getlongestfieldname()
+        for key, value in self.__dict__.items():
+            s += '    {:{width}s}: {}\n'.format(key, value, width=width)
+
+        return s
+    #}}}
+
+    def getfieldnames(self): #{{{
+        return self.__dict__.keys()
+    #}}}
+
+    def getlongestfieldname(self): #{{{
+        maxlength = 0
+        for key in self.__dict__.keys():
+            length = len(key)
+            if length > maxlength:
+                maxlength = length
+
+        return maxlength
+    #}}}
+
+    def setdefaultparameters(self): #{{{
+        return self
+    #}}}
+
+    def checkconsistency(self, md, solution, analyses): #{{{
+        return md
+    #}}}
+
+    def marshall(self, prefix, md, fid): #{{{
+        pass
+    #}}}
+#}}}
Index: /issm/trunk-jpl/src/m/miscellaneous/fielddisplay.py
===================================================================
--- /issm/trunk-jpl/src/m/miscellaneous/fielddisplay.py	(revision 25725)
+++ /issm/trunk-jpl/src/m/miscellaneous/fielddisplay.py	(revision 25726)
@@ -22,5 +22,4 @@
     #string
     if isinstance(field, str):
-
         if len(field) > 30:
             string = displayunit(offset, name, "not displayed", comment)
Index: /issm/trunk-jpl/src/m/miscellaneous/prctile_issm.py
===================================================================
--- /issm/trunk-jpl/src/m/miscellaneous/prctile_issm.py	(revision 25725)
+++ /issm/trunk-jpl/src/m/miscellaneous/prctile_issm.py	(revision 25726)
@@ -3,61 +3,56 @@
 
 
-#  wrapper for prctile to avoid using the matlab statistics toolbox.
 def prctile_issm(x, p, dim):
+    # NumPy has no interpolation method that matches MATLAB's percentile function
+    #y = np.percentile(x, p, dim, interpolation = 'higher')
 
-    try:
-        raise RuntimeError('hello world')
+    if len(np.shape(x)) > 2:
+        raise RuntimeError('Number of dimensions {} not implemented.'.format(len(np.shape(x))))
 
-    #numpy has no interpolation method that matches matlab's percentile function
-    #y = np.percentile(x, p, dim, interpolation = 'higher')
-    except:
-        if len(np.shape(x)) > 2:
-            raise RuntimeError('Number of dimensions  #d not implemented.' + str(len(np.shape(x))))
+    # presumably at least 1 input value has been given
+    #    np.shape(integer) -> (), must be at least (1, )
+    psize = np.shape(p) or (1, )
+    if len(psize) > 1 and np.size(p, 1) > 1:
+        p = p.T
 
-        # presumably at least 1 input value has been given
-        #    np.shape(integer) -> (), must be at least (1, )
-        psize = np.shape(p) or (1, )
-        if len(psize) > 1 and np.size(p, 1) > 1:
-            p = p.T
+    xsize = np.shape(x) or (1, )
+    if dim == 2:
+        x = x.T
 
-        xsize = np.shape(x) or (1, )
-        if dim == 2:
-            x = x.T
+    #  check for any NaN in any columns
+    if not np.isnan(x).any():
+        x = np.sort(x, axis=0)
+        n = np.size(x, 0)
 
-        #  check for any NaN in any columns
-        if not np.isnan(x).any():
-            x = np.sort(x, axis=0)
-            n = np.size(x, 0)
+        #  branch based on number of elements
+        if n > 1:
+            #  set up percent values and interpolate
+            xi = [((i + 0.5) * 100 / n) for i in range(n)]
+            # scipy's interp1d returns a function
+            y = interp1d(xi, x, axis=dim, bounds_error=False)
+            y = y(p)
 
-            #  branch based on number of elements
-            if n > 1:
-                #  set up percent values and interpolate
-                xi = [((i + 0.5) * 100 / n) for i in range(n)]
-                # scipy's interp1d returns a function
-                y = interp1d(xi, x, axis=dim, bounds_error=False)
-                y = y(p)
+            #  fill in high and low values outside of interp range
+            if p > xi[n - 1]:
+                y = np.tile(x[n - 1, :], 1)
+            if p < xi[0]:
+                y = np.tile(x[0, :], 1)
 
-                #  fill in high and low values outside of interp range
-                if p > xi[n - 1]:
-                    y = np.tile(x[n - 1, :], 1)
-                if p < xi[0]:
-                    y = np.tile(x[0, :], 1)
+        #  if one value, just copy it
+        elif n == 1:
+            y = np.tile(x[0, :], (len(p), 1))
 
-            #  if one value, just copy it
-            elif n == 1:
-                y = np.tile(x[0, :], (len(p), 1))
+        #  if no values, use NaN
+        else:
+            y = np.tile(float('NaN'), (np.size(p, 0), np.size(x, 0)))
+    else:
+        #  must loop over columns, since number of elements could be different
+        y = np.zeros((np.size(p, 0), np.size(x, 1)))
+        for j in range(np.size(x, 1)):
+            #  remove any NaN and recursively call column
+            y[:, j] = prctile_issm(x[np.where(not np.isnan(x[:, j]), j)], p)
 
-            #  if no values, use NaN
-            else:
-                y = np.tile(float('NaN'), (np.size(p, 0), np.size(x, 0)))
-        else:
-            #  must loop over columns, since number of elements could be different
-            y = np.zeros((np.size(p, 0), np.size(x, 1)))
-            for j in range(np.size(x, 1)):
-                #  remove any NaN and recursively call column
-                y[:, j] = prctile_issm(x[np.where(not np.isnan(x[:, j]), j)], p)
-
-        if (np.min(xsize) == 1 and len(xsize) > 1 and xsize[dim] > 1 and len(p) > 1 and psize[1] > 1) or (np.min(xsize) > 1 and dim == 2):
-            y = y.T
+    if (np.min(xsize) == 1 and len(xsize) > 1 and xsize[dim] > 1 and len(p) > 1 and psize[1] > 1) or (np.min(xsize) > 1 and dim == 2):
+        y = y.T
 
     return y
Index: /issm/trunk-jpl/src/m/qmu/dakota_in_write.m
===================================================================
--- /issm/trunk-jpl/src/m/qmu/dakota_in_write.m	(revision 25725)
+++ /issm/trunk-jpl/src/m/qmu/dakota_in_write.m	(revision 25726)
@@ -1,45 +1,40 @@
-%
-%  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,corrin,varargin)
-
+%DAKOTA_IN_WRITE - Write a Dakota .in input file.
+%
+%   Usage:
+%      []=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 for 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. (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."
+%
 if ~nargin
     help dakota_in_write
@@ -57,5 +52,5 @@
     dmeth=method;
 else
-    error(['Method ''' inputname(1) ''' is unrecognized class ''' class(method) '''.']);
+    error(['Method ''' inputname(1) ''' is unrecognized class ''' class(method) '''']);
 end
 
@@ -70,8 +65,8 @@
 filei2=fullfile(pathstr,[name ext]);
 
-display(sprintf('Opening Dakota input file ''%s''.',filei2));
+display(sprintf('Opening Dakota input file ''%s''',filei2));
 fidi=fopen(sprintf('%s',filei2),'w');
 if (fidi < 0)
-    error('''%s'' could not be opened.',filei2);
+    error('''%s'' could not be opened',filei2);
 end
 
@@ -110,5 +105,5 @@
 
 fclose(fidi);
-display('End of file successfully written.');
+display('End of file successfully written');
 
 % Uncomment to print contents of Dakota input file (for debugging)
@@ -121,5 +116,5 @@
 function []=strategy_write(fidi,params)
 
-display('Writing strategy section of Dakota input file.');
+display('Writing strategy section of Dakota input file');
 
 fprintf(fidi,'strategy,\n');
@@ -136,5 +131,5 @@
 function []=environment_write(fidi,params)
 
-	display('Writing environment section of Dakota input file.');
+	display('Writing environment section of Dakota input file');
 
 	fprintf(fidi,'environment,\n');
@@ -150,5 +145,5 @@
 function []=method_write(fidi,dmeth,dresp,params)
 
-display('Writing method section of Dakota input file.');
+display('Writing method section of Dakota input file');
 
 fprintf(fidi,'method,\n');
@@ -172,5 +167,5 @@
 function []=model_write(fidi)
 
-display('Writing model section of Dakota input file.');
+display('Writing model section of Dakota input file');
 
 fprintf(fidi,'model,\n');
@@ -183,5 +178,5 @@
 function []=variables_write(fidi,dmeth,dvar,corrin)
 
-display('Writing variables section of Dakota input file.');
+display('Writing variables section of Dakota input file');
 
 fprintf(fidi,'variables,\n');
@@ -219,5 +214,5 @@
 function []=interface_write(fidi,params)
 
-display('Writing interface section of Dakota input file.');
+display('Writing interface section of Dakota input file');
 
 fprintf(fidi,'interface,\n');
@@ -226,5 +221,5 @@
     params.fork=true;
 elseif (params.system+params.fork+params.direct > 1)
-    error('Too many interfaces selected.')
+    error('Too many interfaces selected')
 end
 
@@ -296,5 +291,5 @@
 function []=responses_write(fidi,dmeth,dresp,params)
 
-display('Writing responses section of Dakota input file.');
+display('Writing responses section of Dakota input file');
 
 fprintf(fidi,'responses,\n');
@@ -333,5 +328,5 @@
         params.numerical_gradients=true;
     elseif (params.numerical_gradients+params.analytic_gradients > 1)
-        error('Too many gradients selected.')
+        error('Too many gradients selected')
     end
 
@@ -352,5 +347,5 @@
 
 if ismember('hess',ghspec)
-    error('Hessians needed by method but not provided.');
+    error('Hessians needed by method but not provided');
 else
     fprintf(fidi,'\tno_hessians\n');
Index: /issm/trunk-jpl/src/m/qmu/dakota_in_write.py
===================================================================
--- /issm/trunk-jpl/src/m/qmu/dakota_in_write.py	(revision 25725)
+++ /issm/trunk-jpl/src/m/qmu/dakota_in_write.py	(revision 25726)
@@ -13,32 +13,41 @@
 
 def dakota_in_write(method, dvar, dresp, params, filei, *args):
-    '''
-  write a Dakota .in input file.
-
-  [] = dakota_in_write(method, dvar, dresp, params, filei, args)
-  [] = dakota_in_write(dmeth , dvar, dresp, params, filei, args)
-
-  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 - indepent 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 args are not yet used.
-
-  this function writes a dakota .in input file to be used
-  by dakota.  this file is indepent 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.
-'''
+    """DAKOTA_IN_WRITE - Write a Dakota .in input file.
+
+    Usage:
+        [] = dakota_in_write(method, dvar, dresp, params, filei, args)
+        [] = dakota_in_write(dmeth , dvar, dresp, params, filei, args)
+
+    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 for if empty. params may be empty, in 
+    which case defaults will be used.
+
+    The optional args are not yet used.
+
+    This function writes a Dakota .in input file to be used by Dakota. This 
+    file is indepent 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. (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."
+    """
 
     #  process the input parameters
@@ -63,5 +72,5 @@
     filei2 = fullfile(pathstr, name + ext)
 
-    print('Opening Dakota input file \'' + filei2 + '\'.')
+    print('Opening Dakota input file \'' + filei2 + '\'')
     try:
         with open(filei2, 'w+') as fidi:
@@ -90,7 +99,7 @@
 
     except IOError:
-        print(filei2 + ' could not be opened.')
-
-    print('End of file successfully written.')
+        print(filei2 + ' could not be opened')
+
+    print('End of file successfully written')
 
     # Uncomment to print contents of Dakota input file (for debugging)
@@ -101,5 +110,5 @@
 def strategy_write(fidi, params):
 
-    print('Writing strategy section of Dakota input file.')
+    print('Writing strategy section of Dakota input file')
 
     fidi.write('strategy, \n')
@@ -114,5 +123,5 @@
 def environment_write(fidi, params):
 
-    print('Writing environment section of Dakota input file.')
+    print('Writing environment section of Dakota input file')
 
     fidi.write('environment, \n')
@@ -126,5 +135,5 @@
 def method_write(fidi, dmeth, dresp, params):
 
-    print('Writing method section of Dakota input file.')
+    print('Writing method section of Dakota input file')
 
     fidi.write('method, \n')
@@ -147,5 +156,5 @@
 def model_write(fidi):
 
-    print('Writing model section of Dakota input file.')
+    print('Writing model section of Dakota input file')
 
     fidi.write('model, \n')
@@ -156,5 +165,5 @@
 def variables_write(fidi, dmeth, dvar):
 
-    # print('Writing variables section of Dakota input file.')
+    # print('Writing variables section of Dakota input file')
 
     # fidi.write('variables, \n')
@@ -195,5 +204,5 @@
     # fidi.write('\n')
 
-    print('Writing variables section of Dakota input file.')
+    print('Writing variables section of Dakota input file')
 
     fidi.write('variables, \n')
@@ -218,5 +227,5 @@
 def interface_write(fidi, params):
 
-    print('Writing interface section of Dakota input file.')
+    print('Writing interface section of Dakota input file')
 
     fidi.write('interface, \n')
@@ -225,5 +234,5 @@
         params.fork = True
     elif params.system + params.fork + params.direct > 1:
-        raise RuntimeError('Too many interfaces selected.')
+        raise RuntimeError('Too many interfaces selected')
     if params.system or params.fork:
         param_write(fidi, '\t', 'asynchronous', '', '\n', params)
@@ -289,5 +298,5 @@
 def responses_write(fidi, dmeth, dresp, params):
 
-    print('Writing responses section of Dakota input file.')
+    print('Writing responses section of Dakota input file')
 
     fidi.write('responses, \n')
@@ -320,5 +329,5 @@
             params.numerical_gradients = True
         elif (params.numerical_gradients + params.analytic_gradients > 1):
-            raise RuntimeError('Too many gradients selected.')
+            raise RuntimeError('Too many gradients selected')
 
         if params.numerical_gradients:
@@ -335,5 +344,5 @@
     #  hessians (no implemented methods use them yet)
     if 'hess' in ghspec:
-        raise RuntimeError('Hessians needed by method but not provided.')
+        raise RuntimeError('Hessians needed by method but not provided')
     else:
         fidi.write('\tno_hessians\n')
Index: /issm/trunk-jpl/src/m/qmu/dakota_m_write.m
===================================================================
--- /issm/trunk-jpl/src/m/qmu/dakota_m_write.m	(revision 25725)
+++ /issm/trunk-jpl/src/m/qmu/dakota_m_write.m	(revision 25726)
@@ -1,54 +1,48 @@
-%
-%  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)
-
+%DAKOTA_M_WRITE - Write a Matlab .m function file to be called by Dakota for 
+% the MATLAB direct or external driver.
+%
+%   Usage:
+%      []=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 for 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."
+%
 if ~nargin
     help dakota_m_write
@@ -79,8 +73,8 @@
 filem2=fullfile(pathstr,[name ext versn]);
 
-display(sprintf('Opening Matlab m-file ''%s''.',filem2));
+display(sprintf('Opening Matlab m-file ''%s''',filem2));
 fidm=fopen(sprintf('%s',filem2),'w');
 if (fidm < 0)
-    error('''%s'' could not be opened.',filem2);
+    error('''%s'' could not be opened',filem2);
 end
 
@@ -106,5 +100,5 @@
 
 fclose(fidm);
-display('End of file successfully written.');
+display('End of file successfully written');
 
 end
@@ -114,5 +108,5 @@
 function []=begin_write(fidm,name,params)
 
-display('Writing beginning of Matlab m-file.');
+display('Writing beginning of Matlab m-file');
 
 fprintf(fidm,'%%\n');
@@ -147,5 +141,5 @@
 function []=variables_write(fidm,dmeth,dvar,params,varargin)
 
-display('Writing variables for Matlab m-file.');
+display('Writing variables for Matlab m-file');
 
 fprintf(fidm,'%%  Apply the variables.\n\n');
@@ -174,5 +168,5 @@
 function [ixc]=vlist_write(fidm,ixc,vtype,dvar,params,varargin)
 
-disp(sprintf('  Writing %d %s variables.',length(dvar),class(dvar)));
+disp(sprintf('  Writing %d %s variables',length(dvar),class(dvar)));
 
 for i=1:length(dvar)
@@ -204,5 +198,5 @@
 function []=solution_write(fidm,package)
 
-display('Writing solution for Matlab m-file.');
+display('Writing solution for Matlab m-file');
 fprintf(fidm,'%%  Run the solution.\n\n');
 
@@ -215,5 +209,5 @@
 function []=responses_write(fidm,dmeth,params,dresp)
 
-display('Writing responses for Matlab m-file.');
+display('Writing responses for Matlab m-file');
 
 fprintf(fidm,'%%  Calculate the responses.\n\n');
@@ -251,5 +245,5 @@
 function [ifnvals]=rlist_write(fidm,ifnvals,params,rtype,dresp)
 
-disp(sprintf('  Writing %d %s responses.',length(dresp),class(dresp)));
+disp(sprintf('  Writing %d %s responses',length(dresp),class(dresp)));
 for i=1:length(dresp)
     ifnvals=ifnvals+1;
@@ -267,5 +261,5 @@
 function []=end_write(fidm,name,params)
 
-display('Writing end of Matlab m-file.');
+display('Writing end of Matlab m-file');
 
 fprintf(fidm,'%%  Error condition.\n\n');
Index: /issm/trunk-jpl/src/m/qmu/dakota_out_parse.m
===================================================================
--- /issm/trunk-jpl/src/m/qmu/dakota_out_parse.m	(revision 25725)
+++ /issm/trunk-jpl/src/m/qmu/dakota_out_parse.m	(revision 25726)
@@ -1,51 +1,46 @@
 function [method,dresp,scm,pcm,srcm,prcm]=dakota_out_parse(filei) % {{{
+%DAKOTA_OUT_PARSE - Read a Dakota .out or .dat output file and parse it.
 %
-%  read a Dakota .out or .dat output file and parse it.
+%   Usage:
+%      [method,dresp,scm,pcm,srcm,prcm]=dakota_out_parse(filei)
 %
-%  [method,dresp,scm,pcm,srcm,prcm]=dakota_out_parse(filei)
+%   where the required input is,
+%      filei         (character, name of .out file)
 %
-%  where the required input is:
-%    filei         (character, name of .out file)
+%   the required output is,
+%      method        (character, dakota method name)
+%      dresp         (structure array, responses)
 %
-%  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)
 %
-%  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 for 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.
 %
-%  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: 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 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 postprocessing 
+%   within MATLAB or Excel.
 %
-%  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)
 %
-%  "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."
 %
-%  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."
-%
-
 	if ~nargin
 		help dakota_out_parse
@@ -58,5 +53,5 @@
 	fidi=fopen(sprintf('%s',filei),'r');
 	if (fidi < 0)
-		error('''%s'' could not be opened.',filei);
+		error('''%s'' could not be opened',filei);
 	end
 
@@ -66,5 +61,5 @@
 	fline=fgetl(fidi);
 	if ~ischar(fline)
-		error('File ''%s'' is empty.',filei);
+		error('File ''%s'' is empty',filei);
 	end
 
@@ -90,10 +85,10 @@
 			[ntokens,tokens]=fltokens(fline);
 			method=tokens{1}{1};
-			display(sprintf('Dakota method = ''%s''.',method));
+			display(sprintf('Dakota method = ''%s''',method));
 		elseif strcmp(tokens{1}{1}(7),'N')
 			[fline]=findline(fidi,'methodName = ');
 			[ntokens,tokens]=fltokens(fline);
 			method=tokens{1}{3};
-			display(sprintf('Dakota methodName = ''%s''.',method));
+			display(sprintf('Dakota methodName = ''%s''',method));
 		end
 	end
@@ -165,5 +160,5 @@
 	%     return
 	% end
-	display('End of file successfully reached.');
+	display('End of file successfully reached');
 	fclose(fidi);
 
@@ -172,5 +167,5 @@
 %%  function to parse the dakota tabular output file
 
-	display('Reading Dakota tabular output file.');
+	display('Reading Dakota tabular output file');
 
 	%  process column headings of matrix (skipping eval_id)
@@ -178,11 +173,10 @@
 	[ntokens,tokens]=fltokens(fline);
 
-	% New file DAKOTA versions>6
-	if strncmpi(fline,'%eval_id interface',18)
+	if strncmpi(fline,'%eval_id interface',18) % Dakota versions >= 6
 		offset=2;
-	else %DAKOTA versions<6
+	else % Dakota versions < 6
 		offset=1;
 	end
-	desc=cell (1,ntokens-offset);
+	desc=cell(1,ntokens-offset);
 	data=zeros(1,ntokens-offset);
 
@@ -190,5 +184,5 @@
 		desc(1,i)=cellstr(tokens{1}{i+offset});
 	end
-	display(sprintf('Number of columns (Dakota V+R)=%d.',ntokens-2));
+	display(sprintf('Number of columns (Dakota V + R) = %d',ntokens-2));
 
 	%  process rows of matrix
@@ -197,4 +191,5 @@
 	while 1
 		fline=fgetl(fidi);
+
 		if ~ischar(fline) || isempty(fline)
 			break;
@@ -209,5 +204,5 @@
 		end
 	end
-	display(sprintf('Number of rows (Dakota func evals)=%d.',nrow));
+	display(sprintf('Number of rows (Dakota func evals) = %d',nrow));
 
 	%  calculate statistics
@@ -231,11 +226,11 @@
 	end
 
-	dmin   =min         (data,[],1);
+	dmin   =min(data,[],1);
 	dquart1=prctile_issm(data,25,1);
-	dmedian=median      (data,1);
+	dmedian=median(data,1);
 	dquart3=prctile_issm(data,75,1);
-	dmax   =max         (data,[],1);
-	dmin95=prctile_issm(data,5,1);
-	dmax95=prctile_issm(data,95,1);
+	dmax   =max(data,[],1);
+	dmin95 =prctile_issm(data,5,1);
+	dmax95 =prctile_issm(data,95,1);
 
 	%  same as Dakota scm, Excel correl
@@ -285,5 +280,5 @@
 	[ntokens,tokens]=fltokens(fline);
 	nfeval=tokens{1}{5};
-	display(sprintf('  Dakota function evaluations=%d.',nfeval));
+	display(sprintf('  Dakota function evaluations = %d',nfeval));
 
 end % }}}
@@ -300,5 +295,5 @@
 	[ntokens,tokens]=fltokens(fline);
 	nsamp=tokens{1}{4};
-	display(sprintf('  Dakota samples=%d.',nsamp));
+	display(sprintf('  Dakota samples = %d',nsamp));
 
 end % }}}
@@ -331,5 +326,5 @@
 	end
 
-	display(sprintf('  Number of Dakota response functions=%d.',...
+	display(sprintf('  Number of Dakota response functions = %d',...
 		length(dresp)));
 
@@ -368,5 +363,5 @@
 	end
 
-	display(sprintf('  Number of Dakota response functions=%d.',...
+	display(sprintf('  Number of Dakota response functions = %d',...
 		length(dresp)));
 
@@ -432,5 +427,5 @@
 	end
 
-	display(sprintf('  Number of Dakota response functions=%d.',...
+	display(sprintf('  Number of Dakota response functions = %d',...
 		length(dresp)));
 
@@ -504,5 +499,5 @@
 	end
 
-	display(sprintf('  Number of Dakota response functions=%d.',...
+	display(sprintf('  Number of Dakota response functions = %d',...
 		length(dresp)));
 
@@ -568,5 +563,5 @@
 	end
 
-	display(sprintf('  Number of Dakota response functions=%d.',...
+	display(sprintf('  Number of Dakota response functions = %d',...
 		length(dresp)));
 
@@ -583,5 +578,5 @@
 	end
 
-	display(['Reading ''' fline '''.']);
+	display(['Reading ''' fline '''']);
 
 	cmat.title=fline;
@@ -682,5 +677,5 @@
 
 		if ~idvar
-			display('    Importance Factors not available.');
+			display('    Importance Factors not available');
 			dresp(end).var   ={};
 			dresp(end).impfac=[];
@@ -746,5 +741,5 @@
 
 		if ~icdf
-			display('    Cumulative Distribution Function not available.');
+			display('    Cumulative Distribution Function not available');
 			dresp(ndresp).cdf=[];
 			while ischar(fline) && ...
@@ -757,5 +752,5 @@
 	end
 
-	display(sprintf('  Number of Dakota response functions=%d.',...
+	display(sprintf('  Number of Dakota response functions = %d',...
 		length(dresp)));
 
@@ -902,5 +897,5 @@
 		dresp(end+1).vum=[];
 	end
-	display('Reading measures for volumetric uniformity.');
+	display('Reading measures for volumetric uniformity');
 
 	fline=fgetl(fidi);
@@ -941,5 +936,5 @@
 	[ntokens,tokens]=fltokens(fline);
 	method=tokens{1}{3};
-	display(sprintf('Dakota iterator ''%s'' completed.',method));
+	display(sprintf('Dakota iterator ''%s'' completed',method));
 
 end % }}}
@@ -963,5 +958,5 @@
 
 	warning('findline:str_not_found',...
-		'String ''%s'' not found in file.',string);
+		'String ''%s'' not found in file',string);
 	fseek(fidi,ipos,'bof');
 
Index: /issm/trunk-jpl/src/m/qmu/dakota_out_parse.py
===================================================================
--- /issm/trunk-jpl/src/m/qmu/dakota_out_parse.py	(revision 25725)
+++ /issm/trunk-jpl/src/m/qmu/dakota_out_parse.py	(revision 25726)
@@ -8,44 +8,58 @@
 from helpers import *
 
-#Note: this may be re-written later to take advantage of Python's file i / o mechanics
-#    as it is written now it is often difficult to work with, but is analagous to
-#    the Matlab version of dakota_out_parse
+# NOTE: May be rewritten later to take advantage of Python's file I/O 
+# mechanics. As it is written now, it is often difficult to work with, but is 
+# analagous to the MATLAB version of dakota_out_parse.
 
 
 def dakota_out_parse(filei):  # {{{
-    '''
-  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.
-'''
+    """DAKOTA_OUT_PARSE - read a Dakota .out or .dat output file and parse it.
+
+    Usage:
+        [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 for 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 Python 
+    runtime. 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 postprocessing 
+    within MATLAB or Excel.
+
+    TODO:
+    - Figure out why output from Dakota is different under MATLAB and Python 
+    (is it the input file that we write?)
+
+    "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."
+    """
+
     if filei is None:
         help(dakota_out_parse)
@@ -53,5 +67,5 @@
 
     if not isfile(filei) or getsize(filei) == 0:
-        filei = str(eval(input('Input file?  ')))
+        filei = str(eval(input('Input file? ')))
 
     #fidi = fopen(sprintf('%s', filei), 'r')
@@ -62,7 +76,7 @@
         fline = fidi.readline()
         if getsize(filei) == 0 or fline == '':
-            raise RuntimeError('File ' + filei + ' is empty.')
-
-        dresp = []  # of struct()
+            raise RuntimeError('File ' + filei + ' is empty')
+
+        dresp = [] # of struct()
         scm = struct()
         pcm = struct()
@@ -77,5 +91,5 @@
             fidi.seek(0, 0)
 
-    #  loop through the file to find the Dakota method name
+        # loop through the file to find the Dakota method name
         fline = findline(fidi, 'method', True)
         if fline is None:
@@ -87,19 +101,19 @@
                 [ntokens, tokens] = fltokens(fline)
                 method = tokens[0].strip()
-                print('Dakota method = \'' + method + '\'.')
+                print('Dakota method = \'' + method + '\'')
             elif fline[6] in ['N', 'n']:
                 fline = findline(fidi, 'methodName = ')
                 [ntokens, tokens] = fltokens(fline)
                 method = tokens[2].strip()
-                print('Dakota methodName = \'' + method + '\'.')
-
-    #  loop through the file to find the function evaluation summary
+                print('Dakota methodName = \'' + method + '\'')
+
+        # loop through the file to find the function evaluation summary
         counter = 0
         fline = ''
         nfeval = nfeval_read(fidi, fline)
 
-        #  process each results section based on content of the file
+        # process each results section based on content of the file
         while counter < 10:
-            # because python makes file i / o difficult
+            # because python makes file I/O difficult
             # if we see 10 + blank lines in a row then we have reached EOF
             # (tests show actual maximum number of blank lines is around 5)
@@ -108,5 +122,5 @@
             else:
                 counter = 0
-    #     ipos = ftell(fidi)
+            # ipos = ftell(fidi)
             fline = fidi.readline()
             if fline == '' or fline.isspace():
@@ -147,7 +161,7 @@
                 'Unexpected line: ' + str(fline)
 
-    #     fidi.seek(ipos, 0)
-
-    #  loop through the file to verify the end
+            # fidi.seek(ipos, 0)
+
+    # loop through the file to verify the end
 
     # fline = findline(fidi, '<<<<< Single Method Strategy completed')
@@ -155,10 +169,10 @@
     #     return
     #
-        print('End of file successfully reached.')
+        print('End of file successfully reached')
     #close(fidi)
     #except Exception as err:
     #print "ERROR in dakota_out_parse: " + err
     #raise err
-    #raise RuntimeError(filei + ' could not be opened.')
+    #raise RuntimeError(filei + ' could not be opened')
 
     return [method, dresp, scm, pcm, srcm, prcm]
@@ -167,29 +181,28 @@
 
 def dak_tab_out(fidi, fline):  # {{{
-    #  function to parse the dakota tabular output file
-
-    print('Reading Dakota tabular output file.')
-
-    #  process column headings of matrix (skipping eval_id)
+    """DAK_TAB_OUT - function to parse the Dakota tabular output file
+    """
+
+    print('Reading Dakota tabular output file')
+
+    # Process column headings of matrix (skipping eval_id)
     [ntokens, tokens] = fltokens(fline)
 
-    # New file DAKOTA versions > 6
-    if strncmpi(fline, '%eval_id interface', 18):
+    if strncmpi(fline, '%eval_id interface', 18): # Dakota versions >= 6
         offset = 2
-    else:  #DAKOTA versions < 6
+    else:  # Dakota versions < 6
         offset = 1
 
-    desc = [['' for i in range(ntokens - offset)]]
+    desc = ['' for i in range(ntokens - offset)]
     data = np.zeros((1, ntokens - offset))
 
     for i in range(ntokens - offset):
-        desc[0][i] = str(tokens[i + offset])
-
-    print("Number of columns (Dakota V + R)=" + str(ntokens - 2) + '.')
-
-    #  process rows of matrix
+        desc[i] = str(tokens[i + offset])
+
+    print('Number of columns (Dakota V + R) = {}'.format(ntokens - 2))
+
+    # Process rows of matrix
     nrow = 0
     while True:
-
         fline = fidi.readline()
 
@@ -202,6 +215,5 @@
         [ntokens, tokens] = fltokens(fline)
 
-    #  add row values to matrix (skipping eval_id)
-
+        # Add row values to matrix (skipping eval_id)
         for i in range(ntokens - offset):
             data[nrow, i] = tokens[i + offset]
@@ -209,15 +221,10 @@
         nrow = nrow + 1
 
-    print('Number of rows (Dakota func evals) = ' + str(nrow) + '.')
-
-    #  calculate statistics
-
-    #  since normfit doesn't have a dim argument, and matlab isvector is True
-    #  for a 1xn matrix, handle the case of one row explicitly
-
-    #  Update: normfit_issm.py does handle this case
+    print('Number of rows (Dakota func evals) = ' + str(nrow))
+
+    # Calculate statistics
     if (np.size(data, 0) > 1):
-        #dmean  =mean   (data)
-        #dstddev = std    (data, 0)
+        #dmean = mean(data)
+        #dstddev = std(data, 0)
         [dmean, dstddev, dmeanci, dstddevci] = normfit_issm(data, 0.05)
     else:
@@ -229,31 +236,30 @@
             [dmean[0, i], dstddev[0, i], dmeanci[:, i], dstddevci[:, i]] = normfit_issm(data[:, i], 0.05)
 
-    dmin = data.min(0)
+    dmin = data.min(axis=0)
     dquart1 = prctile_issm(data, 25, 0)
-    dmedian = np.median(data, 0)
+    dmedian = np.median(data, axis=0)
     dquart3 = prctile_issm(data, 75, 0)
-    dmax = data.max(0)
+    dmax = data.max(axis=0)
     dmin95 = prctile_issm(data, 5, 0)
     dmax95 = prctile_issm(data, 95, 0)
 
-    # Note: the following line may cause the following warning
-    #    (should not crash or invalidate results) when one of
-    #    the inputs does not change with respect to the
-    #    other / s causing an internal divide-by - zero error
-
-    # / usr / local / lib / python2.7 / dist - packages / numpy / lib / function_base.py:3163:
-    #    RuntimeWarning: invalid value encountered in true_divide
-    #    c / = stddev[:, None]
-
-    #    (and / or the same but with "c / = stddev[None, :]")
-
+    # NOTE: The following line may cause the following warning (should not 
+    # crash or invalidate results) when one of the inputs does not change with 
+    # respect to the other(s), causing an internal divide-by-zero error,
+    #
+    #       /usr/local/lib/python2.7/dist-packages/numpy/lib/function_base.py:3163:
+    #       RuntimeWarning: invalid value encountered in true_divide
+    #       c /= stddev[:, None]
+    #
+    # (and/or the same but with "c /= stddev[None, :]")
+
+    # Equivalent to Dakota scm, MATLAB corrcoef, and Excel correl
     dcorrel = np.corrcoef(data.T)
 
-    #  divide the data into structures for consistency
+    # Divide the data into structures for consistency
     dresp = []
-
     for i in range(len(desc)):
         dresp.append(struct())
-        dresp[i].descriptor = str(desc[i])
+        dresp[i].descriptor = desc[i]
         dresp[i].sample = data[:, i]
         dresp[i].mean = dmean[i]
@@ -294,5 +300,5 @@
     [ntokens, tokens] = fltokens(fline)
     nfeval = tokens[4]
-    print('  Dakota function evaluations = ' + str(int(nfeval)) + '.')
+    print('  Dakota function evaluations = ' + str(int(nfeval)))
 
     return nfeval
@@ -309,5 +315,5 @@
     [ntokens, tokens] = fltokens(fline)
     nsamp = tokens[3]
-    print('  Dakota samples = ' + str(int(nsamp)) + '.')
+    print('  Dakota samples = ' + str(int(nsamp)))
 
     return nsamp
@@ -340,5 +346,5 @@
         dresp[-1].coefvar = tokens[12]
 
-    print('  Number of Dakota response functions = ' + str(len(dresp)) + '.')
+    print('  Number of Dakota response functions = ' + str(len(dresp)))
 
     return dresp
@@ -350,8 +356,8 @@
 
     if fline is None or fline == '' or fline.isspace():
-        fline = findline(fidi, 'Moment - based statistics for each response function')
+        fline = findline(fidi, 'Moment-based statistics for each response function')
         return
 
-    print('Reading moment - based statistics for response functions:')
+    print('Reading moment-based statistics for response functions:')
 
     #  skip column headings of moment - based statistics
@@ -376,5 +382,5 @@
         dresp[-1].kurtosis = tokens[4]
 
-    print('  Number of Dakota response functions = ' + str(len(dresp)) + '.')
+    print('  Number of Dakota response functions = ' + str(len(dresp)))
 
     return dresp
@@ -433,5 +439,5 @@
             dresp[i].stddevci[1, 0] = tokens[4]
 
-    print('  Number of Dakota response functions = ' + str(len(dresp)) + '.')
+    print('  Number of Dakota response functions = ' + str(len(dresp)))
 
     return dresp
@@ -495,5 +501,5 @@
                 fline = fidi.readline()
 
-    print('  Number of Dakota response functions = ' + str(len(dresp)) + '.')
+    print('  Number of Dakota response functions = ' + str(len(dresp)))
 
     return dresp
@@ -550,5 +556,5 @@
                 fline = fidi.readline()
 
-    print('  Number of Dakota response functions = ' + str(len(dresp)) + '.')
+    print('  Number of Dakota response functions = ' + str(len(dresp)))
 
     return dresp
@@ -565,5 +571,5 @@
             return
 
-    print('Reading ' + fline + '.')
+    print('Reading ' + fline)
 
     cmat.title = fline
@@ -657,5 +663,5 @@
         #  if importance factors missing, skip to cdf
         if not idvar:
-            print('    Importance Factors not available.')
+            print('    Importance Factors not available')
             dresp[-1].var = []
             dresp[-1].impfac = []
@@ -716,10 +722,10 @@
         #  if cdf missing, skip to end of response function
         if not icdf:
-            print('    Cumulative Distribution Function not available.')
+            print('    Cumulative Distribution Function not available')
             dresp[ndresp].cdf = []
             while (fline != '' and not fline.isspace()) and not strncmpi(fline, 'MV Statistics for ', 18) and not strncmp(fline, ' - ', 1):
                 fline = fidi.readline()
 
-    print('  Number of Dakota response functions = ' + str(len(dresp)) + '.')
+    print('  Number of Dakota response functions = ' + str(len(dresp)))
 
     return dresp
@@ -845,5 +851,5 @@
         dresp[-1].vum = []
 
-    print('Reading measures for volumetric uniformity.')
+    print('Reading measures for volumetric uniformity')
     fline = fidi.readline()
     fline = fidi.readline()
@@ -881,5 +887,5 @@
     [ntokens, tokens] = fltokens(fline)
     method = tokens[2]
-    print('Dakota iterator \'' + str(method) + '\' completed.')
+    print('Dakota iterator \'' + str(method) + '\' completed')
 
     return method
@@ -905,5 +911,5 @@
 
     #  issue warning and reset file position
-    print('Warning: findline:str_not_found: String ' + str(string) + ' not found in file.')
+    print('Warning: findline:str_not_found: String ' + str(string) + ' not found in file')
     fidi.seek(ipos, 0)
     return None
Index: /issm/trunk-jpl/src/m/qmu/helpers.py
===================================================================
--- /issm/trunk-jpl/src/m/qmu/helpers.py	(revision 25725)
+++ /issm/trunk-jpl/src/m/qmu/helpers.py	(revision 25726)
@@ -1,15 +1,40 @@
 from collections import OrderedDict
+from copy import deepcopy
+
 import numpy as np
-from copy import deepcopy
 
 
 class struct(object):
-    """An empty struct that can be assigned arbitrary attributes"""
-    pass
+    """STRUCT class definition - An empty struct that can be assigned arbitrary 
+    attributes
+    """
+    def __init__(self): #{{{
+        pass
+    #}}}
+
+    def __repr__(self): #{{{
+        s = ''
+        for key, value in self.__dict__.items():
+            s += '    {}: '.format(key)
+            if isinstance(value, list):
+                s += '[{} element list]'.format(len(value))
+            elif isinstance(value, np.ndarray):
+                if len(value.shape) == 1:
+                    s += '[{} element numpy.ndarray]'.format(value.shape[0])
+                else:
+                    s += '[{}x{} numpy.ndarray]'.format(value.shape[0], value.shape[1])
+            else:
+                s += '{}'.format(value)
+            s += '\n'
+        return s
+    #}}}
+
+    def __len__(self): #{{{
+        return len(self.__dict__.keys())
+    #}}}
 
 
 class Lstruct(list):
-    """
-    An empty struct that can be assigned arbitrary attributes but can also be
+    """An empty struct that can be assigned arbitrary attributes but can also be
     accesed as a list. Eg. x.y = 'hello', x[:] = ['w', 'o', 'r', 'l', 'd']
 
Index: /issm/trunk-jpl/src/m/qmu/postqmu.m
===================================================================
--- /issm/trunk-jpl/src/m/qmu/postqmu.m	(revision 25725)
+++ /issm/trunk-jpl/src/m/qmu/postqmu.m	(revision 25726)
@@ -39,7 +39,8 @@
 	if strcmpi(md.qmu.method.method,'nond_sampling'),
 		dakotaresults.modelresults={};
-		md2=md; md2.qmu.isdakota=0;
+		md2=md;
+		md2.qmu.isdakota=0;
 		for i=1:md2.qmu.method.params.samples,
-			disp(['reading qmu file ' md2.miscellaneous.name '.outbin.' num2str(i)]);
+			disp(['Reading qmu file ' md2.miscellaneous.name '.outbin.' num2str(i)]);
 			md2=loadresultsfromdisk(md2,[md2.miscellaneous.name '.outbin.' num2str(i)]);
 			dakotaresults.modelresults{end+1}=md2.results;
Index: /issm/trunk-jpl/src/m/qmu/postqmu.py
===================================================================
--- /issm/trunk-jpl/src/m/qmu/postqmu.py	(revision 25725)
+++ /issm/trunk-jpl/src/m/qmu/postqmu.py	(revision 25726)
@@ -6,5 +6,6 @@
 from dakota_out_parse import *
 from helpers import *
-import loadresultsfromdisk as loadresults
+import loadresultsfromdisk as loadresultsfromdisk # There is a name conflict somewhere
+from results import results, resultsdakota, solution
 
 
@@ -16,5 +17,4 @@
 
     TODO:
-    - Run Dakota test to check that updates from 6/26 are working
     - Add checks to Popen
     """
@@ -34,10 +34,10 @@
     qmuoutfile = str(md.miscellaneous.name) + '.qmu.out'
     [method, dresp_out, scm, pcm, srcm, prcm] = dakota_out_parse(qmuoutfile)
-    dakotaresults = struct()
+    dakotaresults = resultsdakota()
     dakotaresults.dresp_out = dresp_out
-    dakotaresults.scm = scm
-    dakotaresults.pcm = pcm
-    dakotaresults.srcm = srcm
-    dakotaresults.prcm = prcm
+    dakotaresults.scm       = scm
+    dakotaresults.pcm       = pcm
+    dakotaresults.srcm      = srcm
+    dakotaresults.prcm      = prcm
 
     if isfile('dakota_tabular.dat'):
@@ -50,18 +50,20 @@
             dakotaresults.modelresults = []
             md2 = deepcopy(md)
+            md2.results = results()
             md2.qmu.isdakota = 0
             for i in range(md2.qmu.method.params.samples):
                 outbin_name = '{}.outbin.{}'.format(md2.miscellaneous.name, (i + 1))
-                print('reading qmu file {}'.format(outbin_name))
-                md2 = loadresults.loadresultsfromdisk(md2, outbin_name)
-                dakotaresults.modelresults.append(md2.results)
+                print('Reading qmu file {}'.format(outbin_name))
+                md2 = loadresultsfromdisk.loadresultsfromdisk(md2, outbin_name)
+                dakotaresults.modelresults.append(deepcopy(md2.results))
+            del md2
 
     if md.qmu.statistics.method[0]['name'] != 'None':
         md.qmu.isdakota = 0
-        md = loadresults.loadresultsfromdisk(md, '{}.stats'.format(md.miscellaneous.name))
+        md = loadresultsfromdisk.loadresultsfromdisk(md, '{}.stats'.format(md.miscellaneous.name))
         md.qmu.isdakota = 1
 
     # put dakotaresults in their right location.
-    md.results.dakota = dakotaresults
+    md.results.dakota = deepcopy(dakotaresults)
 
     # move all the individual function evalutations into zip files
Index: /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.py	(revision 25725)
+++ /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.py	(revision 25726)
@@ -1,5 +1,3 @@
 import os
-
-import numpy as np # Remove: used temporarily for debugging
 
 from helpers import fieldnames, struct
@@ -31,8 +29,8 @@
             err_msg += '   Please check for error messages above or in the outlog   \n'
             err_msg += '============================================================\n'
+            print(err_msg)
+            return
 
-            raise OSError(err_msg)
-
-        # Initialize md.results if not a results structure yet
+        # Initialize md.results if not a structure yet
         if not isinstance(md.results, results):
             md.results = results()
Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py	(revision 25725)
+++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py	(revision 25726)
@@ -4,5 +4,5 @@
 import numpy as np
 
-from results import resultselement
+from results import solution
 
 
@@ -29,6 +29,5 @@
     loadres = ReadDataDimensions(fid)
     while loadres:
-
-        #Get time and step
+        # Get time and step
         if loadres['step'] > len(saveres):
             for i in range(len(saveres), loadres['step'] - 1):
@@ -38,8 +37,8 @@
         setattr(saveres[loadres['step'] - 1], 'time', loadres['time'])
 
-    #Add result
+        # Add result
         setattr(saveres[loadres['step'] - 1], loadres['fieldname'], float('NaN'))
 
-    #read next result
+        # Read next result
         loadres = ReadDataDimensions(fid)
 
@@ -65,11 +64,11 @@
         setattr(saveres[loadres['step'] - 1], 'time', loadres['time'])
 
-    #Add result
+        # Add result
         setattr(saveres[loadres['step'] - 1], loadres['fieldname'], loadres['field'])
 
-    #read next result
+        # Read next result
         loadres = ReadData(fid, md)
 
-    #close file
+    # Close file
     fid.close()
 
@@ -108,5 +107,6 @@
 
     # Ok, now construct structure
-    results = resultselement()
+    results = solution()
+
     for i in range(numresults):
         result = allresults[i]
@@ -118,5 +118,4 @@
             setattr(results[index], 'time', result['time'])
         setattr(results[index], result['fieldname'], result['field'])
-
     return results
 # }}}
@@ -322,2 +321,12 @@
     return saveres
 # }}}
+
+def addfieldtorecord(a, descr): #{{{
+    if a.dtype.fields is None:
+        raise ValueError('\'a\' must be a structured numpy array')
+    b = np.empty(a.shape, dtype=a.dtype.descr + descr)
+    for name in a.dtype.names:
+        b[name] = a[name]
+
+    return b
+#}}}