Index: /issm/trunk-jpl/src/m/array/MatlabArray.py =================================================================== --- /issm/trunk-jpl/src/m/array/MatlabArray.py (revision 23095) +++ /issm/trunk-jpl/src/m/array/MatlabArray.py (revision 23095) @@ -0,0 +1,276 @@ +from copy import deepcopy +import numpy as np +from MatlabFuncs import * + +#move this later +from helpers import * + +def allempty(cin): + ''' + function to return an empty cell array if all array elements are empty + cout=allempty(cin) +''' + for i in cin: + if not isempty(i): + cout = cin + return cout + return [] + +def allequal(ain,aval): + ''' + function to return an empty array if all array elements are + equal to the given value, which may also be empty but not nan. + + (note that by definition, nan is not equal to nan; this could + be changed by using isequalwithequalnans.) + + aout=allequal(ain,aval) +''' + if type(ain) != type(aval): + print allequal.__doc__ + raise RuntimeError("ain and aval must be of the same type") + + if type(ain) == list: + ain = np.array(ain) + if type(ain) == np.ndarray: + ain = ain.flatten() + + for i in ain: + if i != aval: + if type(ain) == str: + return '' + else: + return [] + return ain + +def array_numel(*args): + ''' + function to find a number of elements from a list of arrays. + + asize=array_numel(varargin) + + see array_size to check the number and shape of elements, if + multiple indices will be used. +''' + anum = 0 + inum = 0 + for arg in args: + if type(arg) == str: + inum = len(arg) + else: + inum = np.size(arg) + + if inum != 0: + if anum == 0: + anum = inum + else: + if inum != anum: + raise RuntimeError('Inputs had inconsistent number of elements') + return anum + +def array_size(*args): + ''' + function to find an array size from a list of arrays. + + asize=array_size(varargin) + + see array_numel to check only the number of elements, if + single indices will be used. + all arguments are assumed to be 1 or 2 dimensional + + Note: to call on all elements of an array use: array_size(*x) + in Matlab this would be array_size(x{1:end}) + + Note: to get all elements in a linear array use: array_size(np.array(x).flatten()[0:]) + in Matlab this would be array_size(x{1:end}) + because Matlab allows direct 1D access of nD arrays +''' + asize = (0,0) + isize = (0,0) + for arg in args: + if type(arg) == str: + isize = (1,1) #cellstr in matlab makes this happen + else: + isize = np.shape(arg) + if isize == (): #arg is a single value, ex. 0.3, 5, False, etc + isize = (1,1) + + if isize != (0,0): + if asize == (0,0): + asize = isize + else: + if isize != asize: + raise RuntimeError('Inputs had inconsistent shapes') + + #numpy gives (y,) if x = 1, must be reversed to match matlab syntax in this case + if len(asize) == 1: + asize = (1,asize[0]) + + return asize + +def str2int(astr,cfl='first',asint=True): + ''' + function to find and read the first or last positive integer + in a character string. cfl={'first','f','last','l'}; default: 'first' + returns 0 if astr has no positive integers + + Setting asint=False returns a list of strings + eg. ['1','2','3'] from '123' + + aint=str2int(astr,cfl) +''' + aint = [] + + num = '1234567890' + + if type(cfl) != str or type(astr) != str or len(cfl) == 0 or len(astr) == 0: + raise TypeError('str2int(astr,cfl): both arguments must be strings with length > 0') + + # find last positive int + if cfl[0] in ['l','L']: + for i in reversed(astr): + if i in num: + aint.append(i) + else: + if len(aint) > 0: + # aint is backwards since we were iterating backwards + aint = list(reversed(aint)) + if asint: + # convert list(str) to int + aint = int(reduce(lambda x,y: x+y,aint)) + break + + elif cfl[0] in ['f','F']: + for i in astr: + if i in num: + aint.append(i) + else: + if len(aint) > 0: + if asint: + # convert list(str) to int + aint = int(reduce(lambda x,y: x+y,aint)) + break + + # return 0 if aint is still [] (no integers found) + return aint or 0 + +def string_dim(a,idim,*args): + ''' + function to return the string dimension of an array element + + function sdim=string_dim(a,idim,varargin) + + such that: given the array/matrix a, + idim is the linear index of an element in a, + return the x/y/z/w/... coordinates of idim in n dimensions + + ex. a = [1 2 3 + 4 5 6] + + idim = 4 + (a[4] == 5; counted as [1,4,2,5,3,6] linearly in matlab) + + x = string_dim(a,4) -> '[1,1]' + + a[x] == a[4] == a[1,1] == 5 + + example use: exec('print a'+string_dim(a,4)) -> print a[1,1] +''' + sdmin = '' + if type(a) == list: + a = np.array(a) + if type(a) != np.ndarray: + raise TypeError('string_dim(a,idim,*args): a must be a numpy array . Try passing np.array(a) instead') + + if np.size(a) == 0 and idim == 0: + return sdim + + if idim >= np.size(a): + raise RuntimeError('string_dim(a,idim,*args): index idim exceeds number of elements in a') + + #check for column or row vector + for iarg in args: + if strcmpi(iarg,'vector'): + if a.ndim == 2 and (np.shape(a,1) == 1 or np.shape(a,2) == 1): + return '('+str(idim)+')' + + #transpose to compensate for differences in linear indexing in + # matlab vs in python (y/x + linear vs x/y) + a = a.T + + #general case + asize = np.shape(a) + i = np.zeros((np.shape(asize))) + aprod = np.prod(asize) + idim = idim - 1 + index = np.zeros((len(asize))) + + #calculate indices base 0 + for i in range(len(asize)): + aprod=aprod/asize[i] + index[i]=np.floor(idim/aprod) + idim=idim-index[i]*aprod + + #assemble string for output + sdim ='[' + for i in range(len(asize)-1): + sdim += str(int(index[i])) + ',' + + sdim += str(int(index[-1])) + ']' + + # happens due to how python does indexing, this response in matlab is just '' + if sdim == '[-1]': + return '' + + return sdim + +def string_vec(a): + ''' + function to return the string of a vector + + function svec=string_vec(a) +''' + return str(a) + +def struc_class(sclass,cstr,name): + ''' + function to find the structural fields of a specified class + + sclasso=struc_class(sclass,cstr,variable_name) + + such that: + sclasso.variable_name == sclass (hard copy) + + if variable_name == '' + sclasso.cstr == sclass (hard copy) +''' + try: + # I tried other methods, but this is, unfortunately, the best behaving by far + exec 'from '+cstr+' import *' + except: + raise RuntimeError('MatlabArray.struc_class Class Error: class "'+cstr+'" does not exist') + + sclasso = struct() + + if isinstance(sclass,eval(cstr)): + # if we were given no name, call it by its class name + if name != '': + setattr(sclasso, name, deepcopy(sclass)) + else: + setattr(sclasso, cstr, deepcopy(sclass)) + else: + raise RuntimeError('MatlabArray.struc_class Match Error: provided object of type "'+str(type(sclass))+'" does not match provided string; object should be of type '+cstr) + + #may need this later depending on how src/m/classes/qmu works out + + #if len(vars(sclass)) == 0: + #return Lstruct() + + #else: + #fnames = fieldnames(sclass) + #for f in fnames: + #if isinstance(vars(sclass)[f],eval(cstr)): + #exec('sclasso.%s = vars(sclass)[f]')%(f) + #vars(sclasso)[f] = vars(sclass)[f] + + return sclasso Index: /issm/trunk-jpl/src/m/classes/balancethickness.py =================================================================== --- /issm/trunk-jpl/src/m/classes/balancethickness.py (revision 23094) +++ /issm/trunk-jpl/src/m/classes/balancethickness.py (revision 23095) @@ -17,4 +17,6 @@ self.omega = float('NaN') + self.slopex = float('NaN') + self.slopey = float('NaN') #set defaults @@ -56,4 +58,6 @@ WriteData(fid,prefix,'object',self,'fieldname','thickening_rate','format','DoubleMat','mattype',1,'scale',1./yts) WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','omega','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'fieldname','slopex','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'fieldname','slopey','format','DoubleMat','mattype',1) + WriteData(fid,prefix,'object',self,'fieldname','omega','format','DoubleMat','mattype',1) # }}} Index: /issm/trunk-jpl/src/m/classes/clusters/generic.py =================================================================== --- /issm/trunk-jpl/src/m/classes/clusters/generic.py (revision 23094) +++ /issm/trunk-jpl/src/m/classes/clusters/generic.py (revision 23095) @@ -72,6 +72,6 @@ executable='issm.exe'; if isdakota: - version=IssmConfig('_DAKOTA_VERSION_')[0:2] - version=float(version) + version=IssmConfig('_DAKOTA_VERSION_') + version=float(version[0]) if version>=6: executable='issm_dakota.exe' Index: /issm/trunk-jpl/src/m/classes/fourierlove.py =================================================================== --- /issm/trunk-jpl/src/m/classes/fourierlove.py (revision 23094) +++ /issm/trunk-jpl/src/m/classes/fourierlove.py (revision 23095) @@ -12,14 +12,14 @@ def __init__(self): # {{{ - self.nfreq = float('NaN'); - self.frequencies = float('NaN'); - self.sh_nmax = float('NaN'); - self.sh_nmin = float('NaN'); - self.g0 = float('NaN'); - self.r0 = float('NaN'); - self.mu0 = float('NaN'); - self.allow_layer_deletion = float('NaN'); - self.love_kernels = float('NaN'); - self.forcing_type = float('NaN'); + self.nfreq = 0; + self.frequencies = 0; + self.sh_nmax = 0; + self.sh_nmin = 0; + self.g0 = 0; + self.r0 = 0; + self.mu0 = 0; + self.allow_layer_deletion = 0; + self.love_kernels = 0; + self.forcing_type = 0; #set defaults Index: /issm/trunk-jpl/src/m/classes/inversion.py =================================================================== --- /issm/trunk-jpl/src/m/classes/inversion.py (revision 23094) +++ /issm/trunk-jpl/src/m/classes/inversion.py (revision 23095) @@ -163,4 +163,5 @@ WriteData(fid,prefix,'object',self,'fieldname','iscontrol','format','Boolean') WriteData(fid,prefix,'object',self,'fieldname','incomplete_adjoint','format','Boolean') + WriteData(fid,prefix,'object',self,'fieldname','vel_obs','format','DoubleMat','mattype',1,'scale',1./yts) if not self.iscontrol: return Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 23094) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 23095) @@ -218,5 +218,6 @@ return self # }}} - def extract(md,area): # {{{ + #@staticmethod + def extract(self,area): # {{{ """ extract - extract a model according to an Argus contour or flag list @@ -240,5 +241,5 @@ #copy model - md1=copy.deepcopy(md) + md1=copy.deepcopy(self) #get elements that are inside area @@ -366,5 +367,5 @@ #Edges - if md.mesh.domaintype()=='2Dhorizontal': + if md1.mesh.domaintype()=='2Dhorizontal': if np.ndim(md2.mesh.edges)>1 and np.size(md2.mesh.edges,axis=1)>1: #do not use ~isnan because there are some np.nans... #renumber first two columns Index: /issm/trunk-jpl/src/m/classes/qmu.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu.py (revision 23094) +++ /issm/trunk-jpl/src/m/classes/qmu.py (revision 23095) @@ -1,3 +1,5 @@ import numpy as np +from MatlabFuncs import * +from IssmConfig import * from project3d import project3d from collections import OrderedDict @@ -5,4 +7,5 @@ from checkfield import checkfield from WriteData import WriteData +from helpers import * class qmu(object): @@ -16,6 +19,6 @@ def __init__(self): # {{{ self.isdakota = 0 - self.variables = OrderedDict() - self.responses = OrderedDict() + self.variables = OrderedStruct() + self.responses = OrderedStruct() self.method = OrderedDict() self.params = OrderedDict() @@ -117,12 +120,29 @@ return - if not md.qmu.params.evaluation_concurrency==1: - md.checkmessage("concurrency should be set to 1 when running dakota in library mode") - if md.qmu.partition: - if not np.size(md.qmu.partition)==md.mesh.numberofvertices: - md.checkmessage("user supplied partition for qmu analysis should have size md.mesh.numberofvertices x 1") - if not min(md.qmu.partition)==0: + version=IssmConfig('_DAKOTA_VERSION_') + version=float(version[0]) + + if version < 6: + if not md.qmu.params.evaluation_concurrency==1: + md.checkmessage("concurrency should be set to 1 when running dakota in library mode") + else: + if not strcmpi(self.params.evaluation_scheduling,'master'): + md.checkmessage('evaluation_scheduling in qmu.params should be set to "master"') + + if md.cluster.np <= 1: + md.checkmessage('in parallel library mode, Dakota needs to run on at least 2 cpus, 1 cpu for the master, 1 cpu for the slave. Modify md.cluser.np accordingly.') + + if self.params.processors_per_evaluation < 1: + md.checkmessage('in parallel library mode, Dakota needs to run at least one slave on one cpu (md.qmu.params.processors_per_evaluation >=1)!') + + if np.mod(md.cluster.np-1,self.params.processors_per_evaluation): + md.checkmessage('in parallel library mode, the requirement is for md.cluster.np = md.qmu.params.processors_per_evaluation * number_of_slaves, where number_of_slaves will automatically be determined by Dakota. Modify md.cluster.np accordingly') + + if np.size(md.qmu.partition) > 0: + if np.size(md.qmu.partition)!=md.mesh.numberofvertices and np.size(md.qmu.partition) != md.mesh.numberofelements: + md.checkmessage("user supplied partition for qmu analysis should have size (md.mesh.numberofvertices x 1) or (md.mesh.numberofelements x 1)") + if not min(md.qmu.partition.flatten())==0: md.checkmessage("partition vector not indexed from 0 on") - if max(md.qmu.partition)>=md.qmu.numberofpartitions: + if max(md.qmu.partition.flatten())>=md.qmu.numberofpartitions: md.checkmessage("for qmu analysis, partitioning vector cannot go over npart, number of partition areas") @@ -139,5 +159,5 @@ WriteData(fid,prefix,'object',self,'fieldname','variabledescriptors','format','StringArray') WriteData(fid,prefix,'object',self,'fieldname','responsedescriptors','format','StringArray') - if not self.mass_flux_segments: + if not isempty(self.mass_flux_segments): WriteData(fid,prefix,'data',self.mass_flux_segments,'name','md.qmu.mass_flux_segments','format','MatArray'); flag=True; Index: /issm/trunk-jpl/src/m/classes/qmu/@dakota_method/dakota_method.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/@dakota_method/dakota_method.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/@dakota_method/dakota_method.py (revision 23095) @@ -0,0 +1,894 @@ +#move this later +from helpers import * + +from MatlabFuncs import * +import numpy as np + +class dakota_method(object): + ''' + definition for the dakota_method class. + + [dm]=dakota_method(method) + + where the required input is: + method (char, beginning of method name) + + and the output properties and defaults are: + method (char, full method name, '') + type (char, type of method, '') + variables (cell array, applicable variable types, []) + lcspec (cell array, linear constraint specs, []) + responses (cell array, applicable response types, []) + ghspec (cell array, gradient and hessian specs, []) + params (structure, method-depent parameters, []) + + this class is used to guide the writing of a dakota input + file for the specified dakota_method. + + note that zero arguments constructs a default instance one + argument of the class copies the instance and one argument + with enough characters to match a unique method constructs + a new instance of that method. + + "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 np.information + to foreign countries or providing access to foreign persons." + ''' + + def __init__(self,*args): + self.method ='' + self.type ='' + self.variables=[] + self.lcspec =[] + self.responses=[] + self.ghspec =[] + #properites + self.params =struct() + + @staticmethod + def dakota_method(*args): + + dm = dakota_method() + + # return a default object + if len(args) == 0: + return dm + + # copy the object or create the object from the input + elif len(args) == 1: + method = args[0] + + #given argument was a method, copy it + if isinstance(method,dakota_method): + #dm=method + object=method + for field in object.iterkeys(): + if field in vars(dm): + setattr(dm,field,object[field]) + return dm + + #given argument was a way of constructing a method + else: + mlist=[ + 'dot_bfgs', + 'dot_frcg', + 'dot_mmfd', + 'dot_slp', + 'dot_sqp', + 'npsol_sqp', + 'conmin_frcg', + 'conmin_mfd', + 'optpp_cg', + 'optpp_q_newton', + 'optpp_fd_newton', + 'optpp_newton', + 'optpp_pds', + 'asynch_pattern_search', + 'coliny_cobyla', + 'coliny_direct', + 'coliny_ea', + 'coliny_pattern_search', + 'coliny_solis_wets', + 'ncsu_direct', + 'surrogate_based_local', + 'surrogate_based_global', + 'moga', + 'soga', + 'nl2sol', + 'nlssol_sqp', + 'optpp_g_newton', + 'nond_sampling', + 'nond_local_reliability', + 'nond_global_reliability', + 'nond_polynomial_chaos', + 'nond_stoch_collocation', + 'nond_evidence', + 'dace', + 'fsu_quasi_mc', + 'fsu_cvt', + 'vector_parameter_study', + 'list_parameter_study', + 'centered_parameter_study', + 'multidim_parameter_study', + 'bayes_calibration'] + + mlist2=[] + for i in range(len(mlist)): + if strncmpi(method,mlist[i],len(method)): + mlist2.append(mlist[i]) + + # check for a unique match in the list of methods + + l = len(mlist2) + if l == 0: + raise RuntimeError('Unrecognized method: '+str(method)+'.') + elif l == 1: + dm.method=mlist2[0] + else: + raise RuntimeError('Non-unique method: '+str(method)+' matches '+string_cell(mlist2)) + + # assign the default values for the method + #switch dm.method + if dm.method in ['dot_bfgs', + 'dot_frcg']: + dm.type ='dot' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.constraint_tolerance=False + dm.params.output=False + dm.params.speculative=False + dm.params.scaling=False + dm.params.optimization_type='minimize' + elif dm.method in ['dot_mmfd', + 'dot_slp', + 'dot_sqp']: + dm.type ='dot' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =['linear_inequality_constraint', + 'linear_equality_constraint'] + dm.responses=['objective_function', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.constraint_tolerance=False + dm.params.output=False + dm.params.speculative=False + dm.params.scaling=False + dm.params.optimization_type='minimize' + + elif dm.method == 'npsol_sqp': + dm.type ='npsol' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =['linear_inequality_constraint', + 'linear_equality_constraint'] + dm.responses=['objective_function', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.constraint_tolerance=False + dm.params.output=False + dm.params.speculative=False + dm.params.scaling=False + dm.params.verify_level=-1 + dm.params.function_precision=1.0e-10 + dm.params.linesearch_tolerance=0.9 + + elif dm.method == 'conmin_frcg': + dm.type ='conmin' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.constraint_tolerance=False + dm.params.output=False + dm.params.speculative=False + dm.params.scaling=False + elif dm.method == 'conmin_mfd': + dm.type ='conmin' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =['linear_inequality_constraint', + 'linear_equality_constraint'] + dm.responses=['objective_function', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.constraint_tolerance=False + dm.params.output=False + dm.params.speculative=False + dm.params.scaling=False + + elif dm.method == 'optpp_cg': + dm.type ='optpp' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.output=False + dm.params.speculative=False + dm.params.scaling=False + dm.params.max_step=1000. + dm.params.gradient_tolerance=1.0e-4 + elif dm.method in ['optpp_q_newton', + 'optpp_fd_newton', + 'optpp_newton']: + dm.type ='optpp' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =['linear_inequality_constraint', + 'linear_equality_constraint'] + dm.responses=['objective_function', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.output=False + dm.params.speculative=False + dm.params.scaling=False + dm.params.value_based_line_search=False + dm.params.gradient_based_line_search=False + dm.params.trust_region=False + dm.params.tr_pds=False + dm.params.max_step=1000. + dm.params.gradient_tolerance=1.0e-4 + dm.params.merit_function='argaez_tapia' + dm.params.central_path=dm.params.merit_function + dm.params.steplength_to_boundary=0.99995 + dm.params.centering_parameter=0.2 + elif dm.method == 'optpp_pds': + dm.type ='optpp' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.output=False + dm.params.speculative=False + dm.params.scaling=False + dm.params.search_scheme_size=32 + + elif dm.method == 'asynch_pattern_search': + dm.type ='apps' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =['linear_inequality_constraint', + 'linear_equality_constraint'] + dm.responses=['objective_function', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_function_evaluations=False + dm.params.constraint_tolerance=False + dm.params.output=False + dm.params.scaling=False + dm.params.initial_delta=1.0 + dm.params.threshold_delta=0.01 + dm.params.contraction_factor=0.5 + dm.params.solution_target=False + dm.params.synchronization='nonblocking' + dm.params.merit_function='merit2_smooth' + dm.params.constraint_penalty=1.0 + dm.params.smoothing_factor=1.0 + + elif dm.method == 'coliny_cobyla': + dm.type ='coliny' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.output=False + dm.params.scaling=False + dm.params.show_misc_options=False + dm.params.misc_options=[] + dm.params.solution_accuracy=-np.inf + dm.params.initial_delta=[] + dm.params.threshold_delta=[] + elif dm.method == 'coliny_direct': + dm.type ='coliny' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.output=False + dm.params.scaling=False + dm.params.show_misc_options=False + dm.params.misc_options=[] + dm.params.solution_accuracy=-np.inf + dm.params.division='major_dimension' + dm.params.global_balance_parameter=0.0 + dm.params.local_balance_parameter=1.0e-8 + dm.params.max_boxsize_limit=0.0 + dm.params.min_boxsize_limit=0.0001 + dm.params.constraint_penalty=1000.0 + elif dm.method == 'coliny_ea': + dm.type ='coliny' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.output=False + dm.params.scaling=False + dm.params.show_misc_options=False + dm.params.misc_options=[] + dm.params.solution_accuracy=-np.inf + dm.params.seed=False + dm.params.population_size=50 + dm.params.initialization_type='unique_random' + dm.params.fitness_type='linear_rank' + dm.params.replacement_type='elitist' + dm.params.random=[] + dm.params.chc=[] + dm.params.elitist=[] + dm.params.new_solutions_generated='population_size - replacement_size' + dm.params.crossover_type='two_point' + dm.params.crossover_rate=0.8 + dm.params.mutation_type='offset_normal' + dm.params.mutation_scale=0.1 + dm.params.mutation_range=1 + dm.params.dimension_ratio=1.0 + dm.params.mutation_rate=1.0 + dm.params.non_adaptive=False + elif dm.method == 'coliny_pattern_search': + dm.type ='coliny' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.output=False + dm.params.scaling=False + dm.params.show_misc_options=False + dm.params.misc_options=[] + dm.params.solution_accuracy=-np.inf + dm.params.stochastic=False + dm.params.seed=False + dm.params.initial_delta=[] + dm.params.threshold_delta=[] + dm.params.constraint_penalty=1.0 + dm.params.constant_penalty=False + dm.params.pattern_basis='coordinate' + dm.params.total_pattern_size=False + dm.params.no_expansion=False + dm.params.expand_after_success=1 + dm.params.contraction_factor=0.5 + dm.params.synchronization='nonblocking' + dm.params.exploratory_moves='basic_pattern' + elif dm.method == 'coliny_solis_wets': + dm.type ='coliny' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.output=False + dm.params.scaling=False + dm.params.show_misc_options=False + dm.params.misc_options=[] + dm.params.solution_accuracy=-np.inf + dm.params.seed=False + dm.params.initial_delta=[] + dm.params.threshold_delta=[] + dm.params.no_expansion=False + dm.params.expand_after_success=5 + dm.params.contract_after_failure=3 + dm.params.contraction_factor=0.5 + dm.params.constraint_penalty=1.0 + dm.params.constant_penalty=False + + elif dm.method == 'ncsu_direct': + dm.type ='ncsu' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =['linear_inequality_constraint', + 'linear_equality_constraint'] # ? + dm.responses=['objective_function', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] # ? + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.scaling=False + dm.params.solution_accuracy=0. + dm.params.min_boxsize_limit=1.0e-8 + dm.params.vol_boxsize_limit=1.0e-8 + + # if dm.method in ['surrogate_based_local', + # 'surrogate_based_global']: + + elif dm.method == 'moga': + dm.type ='jega' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =['linear_inequality_constraint', + 'linear_equality_constraint'] + dm.responses=['objective_function', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.output=False + dm.params.scaling=False + dm.params.seed=False + dm.params.log_file='JEGAGlobal.log' + dm.params.population_size=50 + dm.params.print_each_pop=False + # according to documentation, uses method-indepent control + # dm.params.output='normal' + dm.params.initialization_type='unique_random' + dm.params.mutation_type='replace_uniform' + dm.params.mutation_scale=0.15 + dm.params.mutation_rate=0.08 + dm.params.replacement_type='' + dm.params.below_limit=6 + dm.params.shrinkage_percentage=0.9 + dm.params.crossover_type='shuffle_random' + dm.params.multi_point_binary=[] + dm.params.multi_point_parameterized_binary=[] + dm.params.multi_point_real=[] + dm.params.shuffle_random=[] + dm.params.num_parents=2 + dm.params.num_offspring=2 + dm.params.crossover_rate=0.8 + dm.params.fitness_type='' + dm.params.niching_type=False + dm.params.radial=[0.01] + dm.params.distance=[0.01] + dm.params.metric_tracker=False + dm.params.percent_change=0.1 + dm.params.num_generations=10 + dm.params.postprocessor_type=False + dm.params.orthogonal_distance=[0.01] + elif dm.method == 'soga': + dm.type ='jega' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =['linear_inequality_constraint', + 'linear_equality_constraint'] + dm.responses=['objective_function', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.output=False + dm.params.scaling=False + dm.params.seed=False + dm.params.log_file='JEGAGlobal.log' + dm.params.population_size=50 + dm.params.print_each_pop=False + dm.params.output='normal' + dm.params.initialization_type='unique_random' + dm.params.mutation_type='replace_uniform' + dm.params.mutation_scale=0.15 + dm.params.mutation_rate=0.08 + dm.params.replacement_type='' + dm.params.below_limit=6 + dm.params.shrinkage_percentage=0.9 + dm.params.crossover_type='shuffle_random' + dm.params.multi_point_binary=[] + dm.params.multi_point_parameterized_binary=[] + dm.params.multi_point_real=[] + dm.params.shuffle_random=[] + dm.params.num_parents=2 + dm.params.num_offspring=2 + dm.params.crossover_rate=0.8 + dm.params.fitness_type='merit_function' + dm.params.constraint_penalty=1.0 + dm.params.replacement_type='' + dm.params.convergence_type=False + dm.params.num_generations=10 + dm.params.percent_change=0.1 + + elif dm.method == 'nl2sol': + dm.type ='lsq' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['least_squares_term'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.output=False + dm.params.scaling=False + dm.params.function_precision=1.0e-10 + dm.params.absolute_conv_tol=-1. + dm.params.x_conv_tol=-1. + dm.params.singular_conv_tol=-1. + dm.params.singular_radius=-1. + dm.params.False_conv_tol=-1. + dm.params.initial_trust_radius=-1. + dm.params.covariance=0 + dm.params.regression_stressbalances=False + elif dm.method == 'nlssol_sqp': + dm.type ='lsq' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =['linear_inequality_constraint', + 'linear_equality_constraint'] + dm.responses=['least_squares_term', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.constraint_tolerance=False + dm.params.output=False + dm.params.speculative=False + dm.params.scaling=False + dm.params.verify_level=-1 + dm.params.function_precision=1.0e-10 + dm.params.linesearch_tolerance=0.9 + elif dm.method == 'optpp_g_newton': + dm.type ='lsq' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =['linear_inequality_constraint', + 'linear_equality_constraint'] + dm.responses=['least_squares_term', + 'nonlinear_inequality_constraint', + 'nonlinear_equality_constraint'] + dm.ghspec =['grad'] + dm.params.max_iterations=False + dm.params.max_function_evaluations=False + dm.params.convergence_tolerance=False + dm.params.output=False + dm.params.speculative=False + dm.params.scaling=False + dm.params.value_based_line_search=False + dm.params.gradient_based_line_search=False + dm.params.trust_region=False + dm.params.tr_pds=False + dm.params.max_step=1000. + dm.params.gradient_tolerance=1.0e-4 + dm.params.merit_function='argaez_tapia' + dm.params.central_path=dm.params.merit_function + dm.params.steplength_to_boundary=0.99995 + dm.params.centering_parameter=0.2 + + elif dm.method == 'nond_sampling': + dm.type ='nond' + dm.variables=['normal_uncertain', + 'uniform_uncertain', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['response_function'] + dm.ghspec =[] + # not documented, but apparently works + dm.params.output=False + dm.params.seed=False + dm.params.fixed_seed=False + dm.params.rng=False + dm.params.samples=False + dm.params.sample_type='lhs' + dm.params.all_variables=False + dm.params.variance_based_decomp=False + dm.params.previous_samples=0 + elif dm.method == 'nond_local_reliability': + dm.type ='nond' + dm.variables=['normal_uncertain', + 'uniform_uncertain', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['response_function'] + dm.ghspec =['grad'] + # not documented, but may work + dm.params.output=False + dm.params.max_iterations=False + dm.params.convergence_tolerance=False + dm.params.mpp_search=False + dm.params.sqp=False + dm.params.nip=False + dm.params.integration='first_order' + dm.params.refinement=False + dm.params.samples=0 + dm.params.seed=False + elif dm.method == 'nond_global_reliability': + dm.type ='nond' + dm.variables=['normal_uncertain', + 'uniform_uncertain', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['response_function'] + dm.ghspec =['grad'] + # not documented, but may work + dm.params.output=False + dm.params.x_gaussian_process=False + dm.params.u_gaussian_process=False + dm.params.all_variables=False + dm.params.seed=False + elif dm.method == 'nond_polynomial_chaos': + dm.type ='nond' + dm.variables=['normal_uncertain', + 'uniform_uncertain', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['response_function'] + dm.ghspec =['grad'] + # not documented, but may work + dm.params.output=False + dm.params.expansion_order=[] + dm.params.expansion_terms=[] + dm.params.quadrature_order=[] + dm.params.sparse_grid_level=[] + dm.params.expansion_samples=[] + dm.params.incremental_lhs=False + dm.params.collocation_points=[] + dm.params.collocation_ratio=[] + dm.params.reuse_samples=False + dm.params.expansion_import_file='' + dm.params.seed=False + dm.params.fixed_seed=False + dm.params.samples=0 + dm.params.sample_type='lhs' + dm.params.all_variables=False + elif dm.method == 'nond_stoch_collocation': + dm.type ='nond' + dm.variables=['normal_uncertain', + 'uniform_uncertain', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['response_function'] + dm.ghspec =['grad'] + # not documented, but may work + dm.params.output=False + dm.params.quadrature_order=[] + dm.params.sparse_grid_level=[] + dm.params.seed=False + dm.params.fixed_seed=False + dm.params.samples=0 + dm.params.sample_type='lhs' + dm.params.all_variables=False + elif dm.method == 'nond_evidence': + dm.type ='nond' + dm.variables=['normal_uncertain', + 'uniform_uncertain', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['response_function'] + dm.ghspec =['grad'] + # not documented, but may work + dm.params.output=False + dm.params.seed=False + dm.params.samples=10000 + + elif dm.method == 'dace': + dm.type ='dace' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function', + 'response_function'] + dm.ghspec =[] + dm.params.grid=False + dm.params.random=False + dm.params.oas=False + dm.params.lhs=False + dm.params.oa_lhs=False + dm.params.box_behnken=False + dm.params.central_composite=False + dm.params.seed=False + dm.params.fixed_seed=False + dm.params.samples=False + dm.params.symbols=False + dm.params.quality_metrics=False + dm.params.variance_based_decomp=False + elif dm.method == 'fsu_quasi_mc': + dm.type ='dace' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function', + 'response_function'] + dm.ghspec =[] + dm.params.halton=False + dm.params.hammersley=False + dm.params.samples=0 + dm.params.sequence_start=[0] + dm.params.sequence_leap=[1] + dm.params.prime_base=False + dm.params.fixed_sequence=False + dm.params.latinize=False + dm.params.variance_based_decomp=False + dm.params.quality_metrics=False + elif dm.method == 'fsu_cvt': + dm.type ='dace' + dm.variables=['continuous_design', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function', + 'response_function'] + dm.ghspec =[] + dm.params.seed=False + dm.params.fixed_seed=False + dm.params.samples=0 + dm.params.num_trials=10000 + dm.params.trial_type='random' + dm.params.latinize=False + dm.params.variance_based_decomp=False + dm.params.quality_metrics=False + + elif dm.method == 'vector_parameter_study': + dm.type ='param' + dm.variables=['continuous_design', + 'normal_uncertain', + 'uniform_uncertain', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function', + 'response_function'] + dm.ghspec =[] + dm.params.output=False + dm.params.final_point=[] + dm.params.step_length=[] + dm.params.num_steps=[] + dm.params.step_vector=[] + dm.params.num_steps=[] + elif dm.method == 'list_parameter_study': + dm.type ='param' + dm.variables=['continuous_design', + 'normal_uncertain', + 'uniform_uncertain', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function', + 'response_function'] + dm.ghspec =[] + dm.params.output=False + dm.params.list_of_points=[] + elif dm.method == 'centered_parameter_study': + dm.type ='param' + dm.variables=['continuous_design', + 'normal_uncertain', + 'uniform_uncertain', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function', + 'response_function'] + dm.ghspec =[] + dm.params.output=False + dm.params.percent_delta=[] + dm.params.deltas_per_variable=[] + elif dm.method == 'multidim_parameter_study': + dm.type ='param' + dm.variables=['continuous_design', + 'normal_uncertain', + 'uniform_uncertain', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function', + 'response_function'] + dm.ghspec =[] + dm.params.output=False + dm.params.partitions=[] + elif dm.method == 'bayes_calibration': + dm.type ='bayes' + dm.variables=['continuous_design', + 'normal_uncertain', + 'uniform_uncertain', + 'continuous_state'] + dm.lcspec =[] + dm.responses=['objective_function', + 'response_function', + 'calibration_function'] + dm.ghspec =[] + dm.params.queso=False + dm.params.dream=False + dm.params.gpmsa=False + dm.params.samples=0 + dm.params.seed=False + dm.params.output=False + + dm.params.metropolis_hastings=False + + dm.params.proposal_covariance=False + dm.params.diagonal=False + dm.params.values=[] + + else: + raise RuntimeError('Unimplemented method: '+str(dm.method)+'.') + + # if more than one argument, issue warning + else: + print 'Warning: dakota_method:extra_arg: Extra arguments for object of class '+str(type(dm))+'.' + return dm + + def __repr__(dm): + + # display the object + string = '\nclass dakota_method object = \n' + string += ' method: '+str(dm.method) + '\n' + string += ' type: '+str(dm.type) + '\n' + string += ' variables: '+str(dm.variables) + '\n' + string += ' lcspec: '+str(dm.lcspec) + '\n' + string += ' responses: '+str(dm.responses) + '\n' + string += ' ghspec: '+str(dm.ghspec) + '\n' + + # display the parameters within the object + + fnames=fieldnames(dm.params) + #get rid of stuff we aren't using + try: + fnames.remove('__module__') + except ValueError: + pass + + maxlen=0 + for i in range(len(fnames)): + maxlen=max(maxlen,len(fnames[i])) + + for i in fnames: + string += ' params.{:{space}s}: {}\n'.format(str(i),str(dm.params.__dict__[i]),space=maxlen+1) + #params.x : y + #with maxlen+1 spaces between x and : + return string + Index: /issm/trunk-jpl/src/m/classes/qmu/@dakota_method/dmeth_params_set.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/@dakota_method/dmeth_params_set.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/@dakota_method/dmeth_params_set.py (revision 23095) @@ -0,0 +1,27 @@ +from helpers import * +from dakota_method import * + +def dmeth_params_set(dm,*args): +# +# set parameters of a dakota_method object. +# +# dm=dmeth_params_set(dm,*args) +# + + if not isinstance(dm,dakota_method): + raise RuntimeError('Provided object is a \''+str(type(dm))+'\' class object, not \'dakota_method\'') + + # loop through each parameter field in the input list + for i in range(0,len(args),2): + if isfield(dm.params,args[i]): + #vars(dresp)[fnames[i]] + exec('dm.params.%s = args[i+1]')%(args[i]) + #vars(dm.params)[args[i]]=args[i+1] + else: + print 'WARNING: dmeth_params_set:unknown_param No parameter \''+str(args[i])+'\' for dakota_method \''+str(dm.method)+'\'.' + + return dm + + + + Index: /issm/trunk-jpl/src/m/classes/qmu/@dakota_method/dmeth_params_write.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/@dakota_method/dmeth_params_write.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/@dakota_method/dmeth_params_write.py (revision 23095) @@ -0,0 +1,587 @@ +from dakota_method import * +from MatlabFuncs import * +from IssmConfig import * + +#move this later: +from helpers import * + +# +# write the parameters from a dakota_method object. +# +# []=dmeth_params_write(dm,fid,sbeg) +# +def dmeth_params_write(dm,fid,sbeg='\t '): + + if not isinstance(dm,dakota_method): + raise RuntimeError('Object '+str(dm)+' is a '+type(dm)+' class object, not .') + + if sbeg == None or sbeg =='': + sbeg='\t ' + + # perform some error checking, but leave the rest to dakota. + # unfortunately this prevents merely looping through the fields + # of the parameters structure. + + # write method-indepent controls + + # param_write(fid,sbeg,'id_method',' = ','\n',dm.params) + # param_write(fid,sbeg,'model_pointer',' = ','\n',dm.params) + + # write method-depent controls + + #switch dm.type + if dm.type == 'dot': + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params) + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params) + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'constraint_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'output',' ','\n',dm.params) + param_write(fid,sbeg,'speculative','','\n',dm.params) + param_write(fid,sbeg,'scaling','','\n',dm.params) + #switch dm.method + if dm.method in ['dot_bfgs', + 'dot_frcg', + 'dot_mmfd', + 'dot_slp', + 'dot_sqp']: + param_write(fid,sbeg,'optimization_type',' = ','\n',dm.params) + + else: + raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + + + elif dm.type == 'npsol': + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params) + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params) + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'constraint_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'output',' ','\n',dm.params) + param_write(fid,sbeg,'speculative','','\n',dm.params) + param_write(fid,sbeg,'scaling','','\n',dm.params) + #switch dm.method + if dm.method == 'npsol_sqp': + param_write(fid,sbeg,'verify_level',' = ','\n',dm.params) + param_write(fid,sbeg,'function_precision',' = ','\n',dm.params) + param_write(fid,sbeg,'linesearch_tolerance',' = ','\n',dm.params) + + else: + raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + + + elif dm.type == 'conmin': + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params) + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params) + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'constraint_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'output',' ','\n',dm.params) + param_write(fid,sbeg,'speculative','','\n',dm.params) + param_write(fid,sbeg,'scaling','','\n',dm.params) + #switch dm.method + if dm.method in ['conmin_frcg', + 'conmin_mfd']: + pass + else: + raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + + + elif dm.type == 'optpp': + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params) + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params) + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'output',' ','\n',dm.params) + param_write(fid,sbeg,'speculative','','\n',dm.params) + param_write(fid,sbeg,'scaling','','\n',dm.params) + #switch dm.method + if dm.method == 'optpp_cg': + param_write(fid,sbeg,'max_step',' = ','\n',dm.params) + param_write(fid,sbeg,'gradient_tolerance',' = ','\n',dm.params) + + elif dm.method in ['optpp_q_newton', + 'optpp_fd_newton', + 'optpp_newton']: + if (dm.params.value_based_line_search + + dm.params.gradient_based_line_search + + dm.params.trust_region + + dm.params.tr_pds > 1): + raise RuntimeError('#s'' method must have only one algorithm.', + dm.method) + + param_write(fid,sbeg,'value_based_line_search','','\n',dm.params) + param_write(fid,sbeg,'gradient_based_line_search','','\n',dm.params) + param_write(fid,sbeg,'trust_region','','\n',dm.params) + param_write(fid,sbeg,'tr_pds','','\n',dm.params) + param_write(fid,sbeg,'max_step',' = ','\n',dm.params) + param_write(fid,sbeg,'gradient_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'merit_function',' = ','\n',dm.params) + param_write(fid,sbeg,'central_path',' = ','\n',dm.params) + param_write(fid,sbeg,'steplength_to_boundary',' = ','\n',dm.params) + param_write(fid,sbeg,'centering_parameter',' = ','\n',dm.params) + + elif dm.method == 'optpp_pds': + param_write(fid,sbeg,'search_scheme_size',' = ','\n',dm.params) + + else: + raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + + + elif dm.type == 'apps': + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params) + param_write(fid,sbeg,'constraint_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'output',' ','\n',dm.params) + param_write(fid,sbeg,'scaling','','\n',dm.params) + #switch dm.method + if dm.method == 'asynch_pattern_search': + param_write(fid,sbeg,'initial_delta',' = ','\n',dm.params) + param_write(fid,sbeg,'threshold_delta',' = ','\n',dm.params) + param_write(fid,sbeg,'contraction_factor',' = ','\n',dm.params) + param_write(fid,sbeg,'solution_target',' = ','\n',dm.params) + param_write(fid,sbeg,'synchronization',' = ','\n',dm.params) + param_write(fid,sbeg,'merit_function',' = ','\n',dm.params) + param_write(fid,sbeg,'constraint_penalty',' = ','\n',dm.params) + param_write(fid,sbeg,'smoothing_factor',' = ','\n',dm.params) + + else: + raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + + + elif dm.type == 'coliny': + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params) + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params) + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'output',' ','\n',dm.params) + param_write(fid,sbeg,'scaling','','\n',dm.params) + + param_write(fid,sbeg,'show_misc_options','','\n',dm.params) + param_write(fid,sbeg,'misc_options',' = ','\n',dm.params) + param_write(fid,sbeg,'solution_accuracy',' = ','\n',dm.params) + #switch dm.method + if dm.method == 'coliny_cobyla': + param_write(fid,sbeg,'initial_delta',' = ','\n',dm.params) + param_write(fid,sbeg,'threshold_delta',' = ','\n',dm.params) + + elif dm.method == 'coliny_direct': + param_write(fid,sbeg,'division',' = ','\n',dm.params) + param_write(fid,sbeg,'global_balance_parameter',' = ','\n',dm.params) + param_write(fid,sbeg,'local_balance_parameter',' = ','\n',dm.params) + param_write(fid,sbeg,'max_boxsize_limit',' = ','\n',dm.params) + param_write(fid,sbeg,'min_boxsize_limit',' = ','\n',dm.params) + param_write(fid,sbeg,'constraint_penalty',' = ','\n',dm.params) + + elif dm.method == 'coliny_ea': + param_write(fid,sbeg,'seed',' = ','\n',dm.params) + param_write(fid,sbeg,'population_size',' = ','\n',dm.params) + param_write(fid,sbeg,'initialization_type',' = ','\n',dm.params) + param_write(fid,sbeg,'fitness_type',' = ','\n',dm.params) + param_write(fid,sbeg,'replacement_type',' = ','\n',dm.params) + param_write(fid,sbeg,'random',' = ','\n',dm.params) + param_write(fid,sbeg,'chc',' = ','\n',dm.params) + param_write(fid,sbeg,'elitist',' = ','\n',dm.params) + param_write(fid,sbeg,'new_solutions_generated',' = ','\n',dm.params) + param_write(fid,sbeg,'crossover_type',' = ','\n',dm.params) + param_write(fid,sbeg,'crossover_rate',' = ','\n',dm.params) + param_write(fid,sbeg,'mutation_type',' = ','\n',dm.params) + param_write(fid,sbeg,'mutation_scale',' = ','\n',dm.params) + param_write(fid,sbeg,'mutation_range',' = ','\n',dm.params) + param_write(fid,sbeg,'dimension_ratio',' = ','\n',dm.params) + param_write(fid,sbeg,'mutation_rate',' = ','\n',dm.params) + param_write(fid,sbeg,'non_adaptive','','\n',dm.params) + + elif dm.method == 'coliny_pattern_search': + param_write(fid,sbeg,'stochastic','','\n',dm.params) + param_write(fid,sbeg,'seed',' = ','\n',dm.params) + param_write(fid,sbeg,'initial_delta',' = ','\n',dm.params) + param_write(fid,sbeg,'threshold_delta',' = ','\n',dm.params) + param_write(fid,sbeg,'constraint_penalty',' = ','\n',dm.params) + param_write(fid,sbeg,'constant_penalty','','\n',dm.params) + param_write(fid,sbeg,'pattern_basis',' = ','\n',dm.params) + param_write(fid,sbeg,'total_pattern_size',' = ','\n',dm.params) + param_write(fid,sbeg,'no_expansion','','\n',dm.params) + param_write(fid,sbeg,'expand_after_success',' = ','\n',dm.params) + param_write(fid,sbeg,'contraction_factor',' = ','\n',dm.params) + param_write(fid,sbeg,'synchronization',' = ','\n',dm.params) + param_write(fid,sbeg,'exploratory_moves',' = ','\n',dm.params) + + elif dm.method == 'coliny_solis_wets': + param_write(fid,sbeg,'seed',' = ','\n',dm.params) + param_write(fid,sbeg,'initial_delta',' = ','\n',dm.params) + param_write(fid,sbeg,'threshold_delta',' = ','\n',dm.params) + param_write(fid,sbeg,'no_expansion','','\n',dm.params) + param_write(fid,sbeg,'expand_after_success',' = ','\n',dm.params) + param_write(fid,sbeg,'contract_after_failure',' = ','\n',dm.params) + param_write(fid,sbeg,'contraction_factor',' = ','\n',dm.params) + param_write(fid,sbeg,'constraint_penalty',' = ','\n',dm.params) + param_write(fid,sbeg,'constant_penalty','','\n',dm.params) + + else: + raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + + + elif dm.type == 'ncsu': + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params) + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params) + param_write(fid,sbeg,'scaling','','\n',dm.params) + #switch dm.method + if dm.method == 'ncsu_direct': + param_write(fid,sbeg,'solution_accuracy',' = ','\n',dm.params) + param_write(fid,sbeg,'min_boxsize_limit',' = ','\n',dm.params) + param_write(fid,sbeg,'vol_boxsize_limit',' = ','\n',dm.params) + + else: + raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + + + elif dm.type == 'jega': + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params) + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params) + param_write(fid,sbeg,'output',' ','\n',dm.params) + param_write(fid,sbeg,'scaling','','\n',dm.params) + + param_write(fid,sbeg,'seed',' = ','\n',dm.params) + param_write(fid,sbeg,'log_file',' = ','\n',dm.params) + param_write(fid,sbeg,'population_size',' = ','\n',dm.params) + param_write(fid,sbeg,'print_each_pop','','\n',dm.params) + param_write(fid,sbeg,'output',' = ','\n',dm.params) + param_write(fid,sbeg,'initialization_type',' = ','\n',dm.params) + param_write(fid,sbeg,'mutation_type',' = ','\n',dm.params) + param_write(fid,sbeg,'mutation_scale',' = ','\n',dm.params) + param_write(fid,sbeg,'mutation_rate',' = ','\n',dm.params) + param_write(fid,sbeg,'replacement_type',' = ','\n',dm.params) + param_write(fid,sbeg,'below_limit',' = ','\n',dm.params) + param_write(fid,sbeg,'shrinkage_percentage',' = ','\n',dm.params) + param_write(fid,sbeg,'crossover_type',' = ','\n',dm.params) + param_write(fid,sbeg,'multi_point_binary',' = ','\n',dm.params) + param_write(fid,sbeg,'multi_point_parameterized_binary',' = ','\n',dm.params) + param_write(fid,sbeg,'multi_point_real',' = ','\n',dm.params) + param_write(fid,sbeg,'shuffle_random',' = ','\n',dm.params) + param_write(fid,sbeg,'num_parents',' = ','\n',dm.params) + param_write(fid,sbeg,'num_offspring',' = ','\n',dm.params) + param_write(fid,sbeg,'crossover_rate',' = ','\n',dm.params) + + #switch dm.method + if dm.method == 'moga': + param_write(fid,sbeg,'fitness_type',' = ','\n',dm.params) + param_write(fid,sbeg,'niching_type',' = ','\n',dm.params) + if not isempty(dm.params.radial) and not isempty(dm.params.distance): + raise RuntimeError('#s'' method must have only one niching distance.', + dm.method) + + param_write(fid,sbeg,'radial',' = ','\n',dm.params) + param_write(fid,sbeg,'distance',' = ','\n',dm.params) + param_write(fid,sbeg,'metric_tracker','','\n',dm.params) + param_write(fid,sbeg,'percent_change',' = ','\n',dm.params) + param_write(fid,sbeg,'num_generations',' = ','\n',dm.params) + param_write(fid,sbeg,'postprocessor_type',' = ','\n',dm.params) + param_write(fid,sbeg,'orthogonal_distance',' = ','\n',dm.params) + + elif dm.method == 'soga': + param_write(fid,sbeg,'fitness_type',' = ','\n',dm.params) + param_write(fid,sbeg,'constraint_penalty',' = ','\n',dm.params) + param_write(fid,sbeg,'replacement_type',' = ','\n',dm.params) + param_write(fid,sbeg,'convergence_type',' = ','\n',dm.params) + param_write(fid,sbeg,'num_generations',' = ','\n',dm.params) + param_write(fid,sbeg,'percent_change',' = ','\n',dm.params) + + else: + raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + + + elif dm.type == 'lsq': + #switch dm.method + if dm.method == 'nl2sol': + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params) + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params) + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'output',' ','\n',dm.params) + param_write(fid,sbeg,'scaling','','\n',dm.params) + + param_write(fid,sbeg,'function_precision',' = ','\n',dm.params) + param_write(fid,sbeg,'absolute_conv_tol',' = ','\n',dm.params) + param_write(fid,sbeg,'x_conv_tol',' = ','\n',dm.params) + param_write(fid,sbeg,'singular_conv_tol',' = ','\n',dm.params) + param_write(fid,sbeg,'singular_radius',' = ','\n',dm.params) + param_write(fid,sbeg,'false_conv_tol',' = ','\n',dm.params) + param_write(fid,sbeg,'initial_trust_radius',' = ','\n',dm.params) + param_write(fid,sbeg,'covariance',' = ','\n',dm.params) + param_write(fid,sbeg,'regression_stressbalances','','\n',dm.params) + + elif dm.method == 'nlssol_sqp': + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params) + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params) + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'constraint_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'output',' ','\n',dm.params) + param_write(fid,sbeg,'speculative','','\n',dm.params) + param_write(fid,sbeg,'scaling','','\n',dm.params) + + param_write(fid,sbeg,'verify_level',' = ','\n',dm.params) + param_write(fid,sbeg,'function_precision',' = ','\n',dm.params) + param_write(fid,sbeg,'linesearch_tolerance',' = ','\n',dm.params) + + elif dm.method == 'optpp_g_newton': + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params) + param_write(fid,sbeg,'max_function_evaluations',' = ','\n',dm.params) + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'output',' ','\n',dm.params) + param_write(fid,sbeg,'speculative','','\n',dm.params) + param_write(fid,sbeg,'scaling','','\n',dm.params) + + if (dm.params.value_based_line_search + + dm.params.gradient_based_line_search + + dm.params.trust_region + + dm.params.tr_pds > 1): + raise RuntimeError('#s'' method must have only one algorithm.', + dm.method) + + param_write(fid,sbeg,'value_based_line_search','','\n',dm.params) + param_write(fid,sbeg,'gradient_based_line_search','','\n',dm.params) + param_write(fid,sbeg,'trust_region','','\n',dm.params) + param_write(fid,sbeg,'tr_pds','','\n',dm.params) + param_write(fid,sbeg,'max_step',' = ','\n',dm.params) + param_write(fid,sbeg,'gradient_tolerance',' = ','\n',dm.params) + param_write(fid,sbeg,'merit_function',' = ','\n',dm.params) + param_write(fid,sbeg,'central_path',' = ','\n',dm.params) + param_write(fid,sbeg,'steplength_to_boundary',' = ','\n',dm.params) + param_write(fid,sbeg,'centering_parameter',' = ','\n',dm.params) + + else: + raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + + + elif dm.type == 'nond': + #switch dm.method + if dm.method == 'nond_sampling': + param_write(fid,sbeg,'seed',' = ','\n',dm.params) + param_write(fid,sbeg,'fixed_seed','','\n',dm.params) + dver = str(IssmConfig('_DAKOTA_VERSION_')[0]) + if ((int(dver[0])==4 and int(dver[2])>2) or int(dver[0])>4): + param_write(fid,sbeg,'rng',' ','\n',dm.params) + + param_write(fid,sbeg,'samples',' = ','\n',dm.params) + param_write(fid,sbeg,'sample_type',' ','\n',dm.params) + param_write(fid,sbeg,'all_variables','','\n',dm.params) + param_write(fid,sbeg,'variance_based_decomp','','\n',dm.params) + if strcmp(dm.params.sample_type,'incremental_random') or strcmp(dm.params.sample_type,'incremental_lhs'): + param_write(fid,sbeg,'previous_samples',' = ','\n',dm.params) + + param_write(fid,sbeg,'output',' ','\n',dm.params) + + elif dm.method == 'nond_local_reliability': + param_write(fid,sbeg,'max_iterations',' = ','\n',dm.params) + param_write(fid,sbeg,'convergence_tolerance',' = ','\n',dm.params) + + param_write(fid,sbeg,'mpp_search',' = ','\n',dm.params) + if type(dm.params.mpp_search) == str: + if (dm.params.sqp + + dm.params.nip > 1): + raise RuntimeError('#s'' method must have only one algorithm.', + dm.method) + + param_write(fid,sbeg,'sqp','','\n',dm.params) + param_write(fid,sbeg,'nip','','\n',dm.params) + param_write(fid,sbeg,'integration',' ','\n',dm.params) + param_write(fid,sbeg,'refinement',' = ','\n',dm.params) + if type(dm.params.refinement) == str: + param_write(fid,sbeg,'samples',' = ','\n',dm.params) + param_write(fid,sbeg,'seed',' = ','\n',dm.params) + + + param_write(fid,sbeg,'output',' ','\n',dm.params) + + elif dm.method == 'nond_global_reliability': + if (dm.params.x_gaussian_process + dm.params.u_gaussian_process != 1): + raise RuntimeError('#s'' method must have one and only one algorithm.', + dm.method) + + param_write(fid,sbeg,'x_gaussian_process','','\n',dm.params) + param_write(fid,sbeg,'u_gaussian_process','','\n',dm.params) + param_write(fid,sbeg,'all_variables','','\n',dm.params) + param_write(fid,sbeg,'seed',' = ','\n',dm.params) + + elif dm.method == 'nond_polynomial_chaos': + param_write(fid,sbeg,'expansion_order',' = ','\n',dm.params) + param_write(fid,sbeg,'expansion_terms',' = ','\n',dm.params) + param_write(fid,sbeg,'quadrature_order',' = ','\n',dm.params) + param_write(fid,sbeg,'sparse_grid_level',' = ','\n',dm.params) + param_write(fid,sbeg,'expansion_samples',' = ','\n',dm.params) + param_write(fid,sbeg,'incremental_lhs','','\n',dm.params) + param_write(fid,sbeg,'collocation_points',' = ','\n',dm.params) + param_write(fid,sbeg,'collocation_ratio',' = ','\n',dm.params) + param_write(fid,sbeg,'reuse_samples','','\n',dm.params) + param_write(fid,sbeg,'expansion_import_file',' = ','\n',dm.params) + param_write(fid,sbeg,'seed',' = ','\n',dm.params) + param_write(fid,sbeg,'fixed_seed','','\n',dm.params) + param_write(fid,sbeg,'samples',' = ','\n',dm.params) + param_write(fid,sbeg,'sample_type',' = ','\n',dm.params) + param_write(fid,sbeg,'all_variables','','\n',dm.params) + + elif dm.method == 'nond_stoch_collocation': + param_write(fid,sbeg,'quadrature_order',' = ','\n',dm.params) + param_write(fid,sbeg,'sparse_grid_level',' = ','\n',dm.params) + param_write(fid,sbeg,'seed',' = ','\n',dm.params) + param_write(fid,sbeg,'fixed_seed','','\n',dm.params) + param_write(fid,sbeg,'samples',' = ','\n',dm.params) + param_write(fid,sbeg,'sample_type',' = ','\n',dm.params) + param_write(fid,sbeg,'all_variables','','\n',dm.params) + + elif dm.method == 'nond_evidence': + param_write(fid,sbeg,'seed',' = ','\n',dm.params) + param_write(fid,sbeg,'samples',' = ','\n',dm.params) + + else: + raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + + + elif dm.type == 'dace': + #switch dm.method + if dm.method == 'dace': + if (dm.params.grid + dm.params.random + dm.params.oas + dm.params.lhs + dm.params.oa_lhs + dm.params.box_behnken + dm.params.central_composite != 1): + raise RuntimeError('#s'' method must have one and only one algorithm.', + dm.method) + + param_write(fid,sbeg,'grid','','\n',dm.params) + param_write(fid,sbeg,'random','','\n',dm.params) + param_write(fid,sbeg,'oas','','\n',dm.params) + param_write(fid,sbeg,'lhs','','\n',dm.params) + param_write(fid,sbeg,'oa_lhs','','\n',dm.params) + param_write(fid,sbeg,'box_behnken','','\n',dm.params) + param_write(fid,sbeg,'central_composite','','\n',dm.params) + param_write(fid,sbeg,'seed',' = ','\n',dm.params) + param_write(fid,sbeg,'fixed_seed','','\n',dm.params) + param_write(fid,sbeg,'samples',' = ','\n',dm.params) + param_write(fid,sbeg,'symbols',' = ','\n',dm.params) + param_write(fid,sbeg,'quality_metrics','','\n',dm.params) + param_write(fid,sbeg,'variance_based_decomp','','\n',dm.params) + + elif dm.method == 'fsu_quasi_mc': + if (dm.params.halton + dm.params.hammersley != 1): + raise RuntimeError('#s'' method must have one and only one sequence type.',dm.method) + + param_write(fid,sbeg,'halton','','\n',dm.params) + param_write(fid,sbeg,'hammersley','','\n',dm.params) + param_write(fid,sbeg,'samples',' = ','\n',dm.params) + param_write(fid,sbeg,'sequence_start',' = ','\n',dm.params) + param_write(fid,sbeg,'sequence_leap',' = ','\n',dm.params) + param_write(fid,sbeg,'prime_base',' = ','\n',dm.params) + param_write(fid,sbeg,'fixed_sequence','','\n',dm.params) + param_write(fid,sbeg,'latinize','','\n',dm.params) + param_write(fid,sbeg,'variance_based_decomp','','\n',dm.params) + param_write(fid,sbeg,'quality_metrics','','\n',dm.params) + + elif dm.method == 'fsu_cvt': + param_write(fid,sbeg,'seed',' = ','\n',dm.params) + param_write(fid,sbeg,'fixed_seed','','\n',dm.params) + param_write(fid,sbeg,'samples',' = ','\n',dm.params) + param_write(fid,sbeg,'num_trials',' = ','\n',dm.params) + param_write(fid,sbeg,'trial_type',' = ','\n',dm.params) + param_write(fid,sbeg,'latinize','','\n',dm.params) + param_write(fid,sbeg,'variance_based_decomp','','\n',dm.params) + param_write(fid,sbeg,'quality_metrics','','\n',dm.params) + + else: + raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + + + elif dm.type == 'param': + param_write(fid,sbeg,'output',' ','\n',dm.params) + #switch dm.method + if dm.method == 'vector_parameter_study': + if not np.logical_xor(isempty(dm.params.final_point),isempty(dm.params.step_vector)): + raise RuntimeError(str(dm.method)+' method must have one and only one specification.') + + if not isempty(dm.params.final_point): + param_write(fid,sbeg,'final_point',' = ','\n',dm.params) + param_write(fid,sbeg,'step_length',' = ','\n',dm.params) + param_write(fid,sbeg,'num_steps',' = ','\n',dm.params) + elif not isempty(dm.params.step_vector): + param_write(fid,sbeg,'step_vector',' = ','\n',dm.params) + param_write(fid,sbeg,'num_steps',' = ','\n',dm.params) + + + elif dm.method == 'list_parameter_study': + param_write(fid,sbeg,'list_of_points',' = ','\n',dm.params) + + elif dm.method == 'centered_parameter_study': + param_write(fid,sbeg,'percent_delta',' = ','\n',dm.params) + param_write(fid,sbeg,'deltas_per_variable',' = ','\n',dm.params) + + elif dm.method == 'multidim_parameter_study': + param_write(fid,sbeg,'partitions',' = ','\n',dm.params) + + else: + raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + + + elif dm.type == 'bayes': + #switch dm.method + if dm.method == 'bayes_calibration': + # if (dm.params.queso + + # dm.params.dream + + # dm.params.gpmsa ~= 1) + # raise RuntimeError('''#s'' method must have one and only one bayes type. YOU SUCK', + # dm.method) + # + param_write(fid,sbeg,'queso','','\n',dm.params) + param_write(fid,sbeg,'dream','','\n',dm.params) + param_write(fid,sbeg,'gpmsa','','\n',dm.params) + param_write(fid,sbeg,'samples',' = ','\n',dm.params) + param_write(fid,sbeg,'seed',' = ','\n',dm.params) + param_write(fid,sbeg,'output',' =','\n',dm.params) + param_write(fid,sbeg,'metropolis_hastings','','\n',dm.params) + param_write(fid,sbeg,'proposal_covariance','','\n',dm.params) + param_write(fid,sbeg,'diagonal','','\n',dm.params) + param_write(fid,sbeg,'values',' = ','\n',dm.params) + + else: + raise RuntimeError('Unrecognized '+dm.type+' method: '+dm.method+'.') + + +## function to write a structure of parameters + +def param_struc_write(fidi,sbeg,smid,s,params): + + # loop through each parameter field in the structure + + fnames=fieldnames(params) + + for i in range(np.size(fnames)): + param_write(fidi,sbeg,fnames[i],smid,s,params) + + return + +## function to write a parameter + +def param_write(fidi,sbeg,pname,smid,s,params): + + # check for errors + + if not isfield(params,pname): + warning('param_write:param_not_found', + 'Parameter '+str(pname)+' not found in '+params+'.') + return + elif type(vars(params)[pname]) == bool and not vars(params)[pname]: + return + elif isempty(vars(params)[pname]): + print 'Warning: param_write:param_empty: Parameter '+pname+' requires input of type '+type(vars(params)[pname])+'.' + return + + + # construct the parameter string based on type + + if type(vars(params)[pname]) == bool: + fidi.write(sbeg+str(pname)+s) + elif type(vars(params)[pname]) in [int,float]: + fidi.write(sbeg+str(pname)+smid+str(vars(params)[pname])+s) + elif type(vars(params)[pname]) == list: + fidi.write(sbeg+str(pname)+smid+str(vars(params)[pname][0])) + for i in range(1,np.size(vars(params)[pname])): + fidi.write(' '+str(vars(params)[pname][i])) + + fidi.write(s) + elif type(vars(params)[pname]) == str: + fidi.write(sbeg+str(pname)+smid+str(vars(params)[pname])+s) + else: + print 'Warning: param_write:param_unrecog: Parameter '+pname+' is of unrecognized type '+type(vars(params)[pname])+'.' + return + + + Index: /issm/trunk-jpl/src/m/classes/qmu/calibration_function.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/calibration_function.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/calibration_function.py (revision 23095) @@ -0,0 +1,143 @@ +import numpy as np +from rlist_write import * +from MatlabArray import * + +class calibration_function(object): + ''' + definition for the calibration_function class. + + [cf] = calibration_function.calibration_function(args) + cf = calibration_function() + + where the required args are: + descriptor (char, description, '') + and the optional args and defaults are: + scale_type (char, scaling type, 'none') + scale (double, scaling factor, 1.) + weight (double, weighting factor, 1.) + + note that zero arguments constructs a default instance, one + argument of the class copies the instance, and one or more + arguments constructs a new instance from the arguments. +''' + def __init__(self): + self.descriptor = '' + self.scale_type = 'none' + self.scale = 1. + self.weight = 1. + + @staticmethod + def calibration_function(*args): + nargin = len(args) + + # create a default object + if nargin == 0: + return calibration_function() + + # copy the object or create the object from the input + else: + if (nargin == 1) and isinstance(args[0],calibration_function): + cf = args[0] + else: + asizec = array_size(*args[0:min(nargin,4)]) + cf = [calibration_function() for i in range(asizec[0]) for j in range(asizec[1])] + + for i in range(np.size(cf)): + if (np.size(args[0]) > 1): + cf[i].descriptor = args[0][i] + else: + cf[i].descriptor = str(args[0])+string_dim(cf,i,'vector') + + if nargin >= 2: + for i in range(np.size(cf)): + cf[i].scale_type = str(args[1]) + if nargin >= 3: + for i in range(np.size(cf)): + cf[i].scale = args[2] + if nargin >= 4: + for i in range(np.size(cf)): + cf[i].weight = args[3] + if nargin > 4: + print 'WARNING: calibration_function:extra_arg: Extra arguments for object of class '+str(type(cf))+'.' + + return cf + + def __repr__(self): + # display the object + string = '\n' + string += 'class "calibration_function" object = \n' + string += ' descriptor: '+str(self.descriptor) + '\n' + string += ' scale_type: '+str(self.scale_type) + '\n' + string += ' scale: '+str(self.scale) + '\n' + string += ' weight: '+str(self.weight) + '\n' + + return string + + # from here on, cf is either a single, or a 1d vector of, calibration_function + + @staticmethod + def prop_desc(cf,dstr): + if type(cf) not in [list,np.ndarray]: + if cf.descriptor != '' or type(cf.descriptor) != str: + desc = str(cf.descriptor) + elif dstr != '': + desc = str(dstr) + else: + desc = 'cf' + return desc + + desc = ['' for i in range(np.size(cf))] + for i in range(np.size(cf)): + if cf[i].descriptor != '' or type(cf[i].descriptor) != str: + desc[i] = str(cf[i].descriptor) + elif dstr != '': + desc[i] = str(dstr)+str(string_dim(cf,i,'vector')) + else: + desc[i] = 'cf'+str(string_dim(cf,i,'vector')) + + desc = allempty(desc) + + return desc + + @staticmethod + def prop_stype(cf): + stype='' + return stype + + @staticmethod + def prop_weight(cf): + weight=[] + return weight + + @staticmethod + def prop_lower(cf): + lower=[] + return lower + + @staticmethod + def prop_upper(cf): + upper=[] + return upper + + @staticmethod + def prop_target(cf): + target=[] + return target + + @staticmethod + def prop_scale(cf): + scale=[] + return scale + + @staticmethod + def dakota_write(fidi,dresp,rdesc): + # collect only the responses of the appropriate class + cf = [struc_class(i,'calibration_function','cf') for i in dresp] + + # write responses + rdesc = rlist_write(fidi,'calibration_terms','calibration_function',cf,rdesc) + return rdesc + + @staticmethod + def dakota_rlev_write(fidi,dresp,params): + return Index: /issm/trunk-jpl/src/m/classes/qmu/continuous_design.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/continuous_design.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/continuous_design.py (revision 23095) @@ -0,0 +1,231 @@ +import numpy as np +from vlist_write import * +from MatlabArray import * + +class continuous_design(object): + ''' + definition for the continuous_design class. + + [cdv] = continuous_design.continuous_design(args) + cdv = continuous_design() + + where the required args are: + descriptor (char, description, '') + initpt (double, initial point, 0.) + and the optional args and defaults are: + lower (double, lower bound, -Inf) + upper (double, upper bound, Inf) + scale_type (char, scaling type, 'none') + scale (double, scaling factor, 1.) + + note that zero arguments constructs a default instance, one + argument of the class copies the instance, and two or more + arguments constructs a new instance from the arguments. +''' + + def __init__(self): + self.descriptor = '' + self.initpt = 0. + self.lower = -np.inf + self.upper = np.inf + self.scale_type = 'none' + self.scale = 1. + + @staticmethod + def continuous_design(*args): + nargin = len(args) + + # create a default object + if nargin == 0: + return continuous_design() + + # copy the object + if nargin == 1: + if isinstance(args[0],continuous_design): + cdv = args[0] + else: + raise RuntimeError('Object is a '+str(type(args[0]))+' class object, not "continuous_design".') + + # create the object from the input + else: + shapec = array_size(*args[0:min(nargin,6)]) + cdv = [continuous_design() for i in range(shapec[0]) for j in range(shapec[1])] + # in case cdv doesn't use array-like args + #cdv = continuous_design() + + for i in range(np.size(cdv)): + if (np.size(args[0]) > 1): + cdv[i].descriptor = args[0][i] + else: + cdv[i].descriptor = str(args[0])+string_dim(cdv,i,'vector') + + if (nargin >= 2): + for i in range(np.size(cdv)): + if (np.size(args[1]) > 1): + cdv[i].initpt = args[1][i] + else: + cdv[i].initpt = args[1] + + if (nargin >= 3): + for i in range(np.size(cdv)): + if (np.size(args[2]) > 1): + cdv[i].lower = args[2][i] + else: + cdv[i].lower = args[2] + + if (nargin >= 4): + for i in range(np.size(cdv)): + if (np.size(args[3]) > 1): + cdv[i].upper = args[3][i] + else: + cdv[i].upper = args[3] + + if (nargin >= 5): + for i in range(np.size(cdv)): + if (np.size(args[4]) > 1): + cdv[i].scale_type = args[4][i] + else: + cdv[i].scale_type = str(args[4]) + + if (nargin >= 6): + for i in range(np.size(cdv)): + if (np.size(args[5]) > 1): + cdv[i].scale = args[5][i] + else: + cdv[i].scale = args[5] + + if (nargin > 6): + print 'WARNING: continuous_design:extra_arg: Extra arguments for object of class '+str(type(cdv))+'.' + + return cdv + + + def __repr__(self): + # display the object + string = '\n' + string += 'class "continuous_design" object = \n' + string += ' descriptor: ' +str(self.descriptor) + '\n' + string += ' initpt: ' +str(self.initpt) + '\n' + string += ' lower: ' +str(self.lower) + '\n' + string += ' upper: ' +str(self.upper) + '\n' + string += ' scale_type: ' +str(self.scale_type) + '\n' + string += ' scale: ' +str(self.scale) + '\n' + + return string + + @staticmethod + def prop_desc(cdv,dstr): + if type(cdv) not in [list,np.ndarray]: + cdv = [cdv] + # in case cdv doesn't use array-like args + #if cdv.descriptor != '' or type(cdv.descriptor) != str: + #desc = str(cdv.descriptor) + #elif dstr != '': + #desc = str(dstr) + #else: + #desc = 'cdv' + #return desc + + desc = ['' for i in range(np.size(cdv))] + for i in range(np.size(cdv)): + if cdv[i].descriptor != '' or type(cdv[i].descriptor) != str: + desc[i] = str(cdv[i].descriptor) + elif dstr != '': + desc[i] = str(dstr)+str(string_dim(cdv,i,'vector')) + else: + desc[i] = 'cdv'+str(string_dim(cdv,i,'vector')) + + desc = allempty(desc) + + return desc + + @staticmethod + def prop_initpt(cdv): + if type(cdv) not in [list,np.ndarray]: + return cdv.initpt + + initpt = np.zeros(np.size(cdv)) + for i in range(np.size(cdv)): + initpt[i] = cdv[i].initpt + + initpt = allequal(initpt,0.) + + return initpt + + @staticmethod + def prop_lower(cdv): + if type(cdv) not in [list,np.ndarray]: + return cdv.lower + + lower = np.zeros(np.size(cdv)) + for i in range(np.size(cdv)): + lower[i] = cdv[i].lower + + lower = allequal(lower,-np.inf) + + return lower + + @staticmethod + def prop_upper(cdv): + if type(cdv) not in [list,np.ndarray]: + return cdv.upper + + upper = np.zeros(np.size(cdv)) + for i in range(np.size(cdv)): + upper[i] = cdv[i].upper + + upper = allequal(upper, np.inf) + + return upper + + @staticmethod + def prop_mean(cdv): + mean=[] + return mean + + @staticmethod + def prop_stddev(cdv): + stddev=[] + return sttdev + + @staticmethod + def prop_initst(cdv): + initst=[] + return initst + + @staticmethod + def prop_stype(cdv): + if type(cdv) not in [list,np.ndarray]: + return str(cdv.scale_type) + + stype = np.empty(np.size(cdv)) + stype.fill(0.0) + for i in range(np.size(cdv)): + stype[i] = str(cdv[i].scale_type) + + stype = allequal(stype,'none') + + return stype + + @staticmethod + def prop_scale(cdv): + if type(cdv) not in [list,np.ndarray]: + return cdv.scale + + scale = np.zeros(np.size(cdv)) + for i in range(np.size(cdv)): + scale[i] = cdv[i].scale + + scale = allequal(scale,1.) + + return scale + + + @staticmethod + def dakota_write(fidi,dvar): + # collect only the variables of the appropriate class + cdv = [struc_class(i,'continuous_design','cdv') for i in dvar] + + # write variables + vlist_write(fidi,'continuous_design','cdv',cdv) + Index: /issm/trunk-jpl/src/m/classes/qmu/continuous_state.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/continuous_state.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/continuous_state.py (revision 23095) @@ -0,0 +1,182 @@ +import numpy as np +from vlist_write import * +from MatlabArray import * + +class continuous_state(object): + ''' + definition for the continuous_state class. + + [csv] = continuous_state.continuous_state(args) + csv = continuous_state() + + where the required args are: + descriptor (char, description, '') + initst (double, initial state, 0.) + and the optional args and defaults are: + lower (double, lower bound, -Inf) + upper (double, upper bound, Inf) + + note that zero arguments constructs a default instance, one + argument of the class copies the instance, and two or more + arguments constructs a new instance from the arguments. +''' + + def __init__(self): + self.descriptor = '' + self.initst = 0. + self.lower = -np.inf + self.upper = np.inf + + @staticmethod + def continuous_state(*args): + nargin = len(args) + + # create a default object + if nargin == 0: + return continuous_state() + + # copy the object + if nargin == 1: + if isinstance(args[0],continuous_state): + csv = args[0] + else: + raise RuntimeError('Object is a '+str(type(args[0]))+' class object, not "continuous_state".') + + # create the object from the input + else: + shapec = np.shape(*args[0:min(nargin,4)]) + csv = [continuous_state() for i in range(shapec[0]) for j in range(shapec[1])] + + for i in range(np.size(csv)): + if (np.size(args[0]) > 1): + csv[i].descriptor = args[0][i] + else: + csv[i].descriptor = str(args[0])+string_dim(csv,i,'vector') + + if (nargin >= 2): + for i in range(np.size(csv)): + if (np.size(args[1]) > 1): + csv[i].initst = args[1][i] + else: + csv[i].initst = args[1] + + if (nargin >= 3): + for i in range(np.size(csv)): + if (np.size(args[2]) > 1): + csv[i].lower = args[2][i] + else: + csv[i].lower = args[2] + + if (nargin >= 4): + for i in range(np.size(csv)): + if (np.size(args[3]) > 1): + csv[i].upper = args[3][i] + else: + csv[i].upper = args[3] + + if (nargin > 4): + print 'continuous_state:extra_arg','Extra arguments for object of class '+str(type(csv))+'.' + + return csv + + def __repr__(self): + # display the object + string = '\n' + string += 'class "continuous_state" object = \n' + string += ' descriptor: ' +str(self.descriptor) + '\n' + string += ' initst: ' +str(self.initst) + '\n' + string += ' lower: ' +str(self.lower) + '\n' + string += ' upper: ' +str(self.upper) + '\n' + + return string + + @staticmethod + def prop_desc(csv,dstr): + if type(csv) not in [list,np.ndarray]: + csv = [csv] + + desc = ['' for i in range(np.size(csv))] + for i in range(np.size(csv)): + if csv[i].descriptor != '' or type(cdv[i].descriptor) != str: + desc[i] = str(csv[i].descriptor) + elif dstr != '': + desc[i] = str(dstr)+str(string_dim(csv,i,'vector')) + else: + desc[i] = 'csv'+str(string_dim(csv,i,'vector')) + + desc = allempty(desc) + + return desc + + @staticmethod + def prop_initpt(csv): + initpt=[] + return initpt + + @staticmethod + def prop_lower(csv): + if type(csv) not in [list,np.ndarray]: + return csv.lower + + lower = np.zeros(np.size(csv)) + for i in range(np.size(csv)): + lower[i] = csv[i].lower + + lower = allequal(lower,-np.inf) + + return lower + + @staticmethod + def prop_upper(csv): + if type(csv) not in [list,np.ndarray]: + return csv.upper + + upper = np.zeros(np.size(csv)) + for i in range(np.size(csv)): + upper[i] = csv[i].upper + + upper = allequal(upper, np.inf) + + return upper + + @staticmethod + def prop_mean(csv): + mean=[] + return mean + + @staticmethod + def prop_stddev(csv): + stddev=[] + return sttdev + + @staticmethod + def prop_initst(csv): + if type(csv) not in [list,np.ndarray]: + return csv.initst + + initst = np.zeros(np.size(csv)) + for i in range(np.size(csv)): + initst[i] = csv[i].initst + + initst = allequal(initst,0.) + + return initst + + @staticmethod + def prop_stype(csv): + stype='' + return stype + + @staticmethod + def prop_scale(csv): + scale=[] + return scale + + @staticmethod + def dakota_write(fidi,dvar): + # collect only the variables of the appropriate class + csv = [struc_class(i,'continuous_state','csv') for i in dvar] + + # write variables + vlist_write(fidi,'continuous_state','csv',csv) + Index: /issm/trunk-jpl/src/m/classes/qmu/least_squares_term.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/least_squares_term.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/least_squares_term.py (revision 23095) @@ -0,0 +1,171 @@ +import numpy as np +from rlist_write import * +from MatlabArray import * + +class least_squares_term(object): + ''' + definition for the least_squares_term class. + + [lst] = least_squares_term.least_squares_term(args) + lst = least_squares_term() + + where the required args are: + descriptor (char, description, '') + and the optional args and defaults are: + scale_type (char, scaling type, 'none') + scale (double, scaling factor, 1.) + weight (double, weighting factor, 1.) + + note that zero arguments constructs a default instance, one + argument of the class copies the instance, and one or more + arguments constructs a new instance from the arguments. +''' + def __init__(self): + self.descriptor = '' + self.scale_type = 'none' + self.scale = 1. + self.weight = 1. + + @staticmethod + def least_squares_term(*args): + nargin = len(args) + + #create a default object + if nargin == 0: + return least_squares_term() + + #copy the object or create the object from the input + else: + if (nargin == 1) and isinstance(args[0],least_squares_term): + lst = args[0] + else: + asizec = np.shape(*args[0:min(nargin,4)]) + lst = [least_squares_term() for i in range(asizec[0]) for j in range(asizec[1])] + + for i in range(np.size(lst)): + if (np.size(args[0]) > 1): + lst[i].descriptor = args[0][i] + else: + lst[i].descriptor = str(args[0])+string_dim(lst,i,'vector') + + if (nargin >= 2): + for i in range(np.size(lst)): + if (np.size(args[1]) > 1): + lst[i].scale_type = args[1][i] + else: + lst[i].scale_type = str(args[1]) + + if (nargin >= 3): + for i in range(np.size(lst)): + if (np.size(args[2]) > 1): + lst[i].scale = args[2][i] + else: + lst[i].scale = args[2] + + if (nargin >= 4): + for i in range(np.size(lst)): + if (np.size(args[3]) > 1): + lst[i].weight = args[3][i] + else: + lst[i].weight = args[3] + + if (nargin > 4): + print 'WARNING: least_squares_term:extra_arg Extra arguments for object of class '+str(type(lst))+'.' + + return lst + + def __repr__(self): + # display the object + string = '\n' + string += 'class "least_squares_term" object = \n' + string += ' descriptor: '+str(self.descriptor) + '\n' + string += ' scale_type: '+str(self.scale_type) + '\n' + string += ' scale: '+str(self.scale) + '\n' + string += ' weight: '+str(self.weight) + '\n' + + return string + + @staticmethod + def prop_desc(lst,dstr): + if type(lst) not in [list,np.ndarray]: + lst = [lst] + + desc = ['' for i in range(np.size(lst))] + for i in range(np.size(lst)): + if lst[i].descriptor != '' or type(cdv[i].descriptor) != str: + desc[i] = str(lst[i].descriptor) + elif dstr != '': + desc[i] = str(dstr)+str(string_dim(lst,i,'vector')) + else: + desc[i] = 'lst'+str(string_dim(lst,i,'vector')) + + desc = allempty(desc) + + return desc + + @staticmethod + def prop_stype(lst): + if type(lst) not in [list,np.ndarray]: + return str(lst.scale_type) + + stype = ['' for i in range(np.size(lst))] + for i in range(np.size(lst)): + stype[i] = str(lst[i].scale_type) + + stype = allequal(stype,'none') + + return stype + + @staticmethod + def prop_scale(lst): + if type(lst) not in [list,np.ndarray]: + return lst.scale + + scale = np.zeros(np.size(lst)) + for i in range(np.size(lst)): + scale[i] = lst[i].scale + + scale = allequal(scale,1.) + + return scale + + @staticmethod + def prop_weight(lst): + if type(lst) not in [list,np.ndarray]: + return lst.weight + + weight = np.zeros(np.size(lst)) + for i in range(np.size(lst)): + weight[i] = lst[i].weight + + weight = allequal(weight,1.) + + return weight + + @staticmethod + def prop_lower(lst): + lower=[] + return lower + + @staticmethod + def prop_upper(lst): + upper=[] + return upper + + @staticmethod + def prop_target(lst): + target=[] + return target + + @staticmethod + def dakota_write(fidi,dresp,rdesc): + #collect only the responses of the appropriate class + lst = [struc_class(i,'least_squares_term','lst') for i in dresp] + + #write responses + rdesc = rlist_write(fidi,'least_squares_terms','least_squares_term',lst,rdesc) + return rdesc + + @staticmethod + def dakota_rlev_write(fidi,dresp,params): + return Index: /issm/trunk-jpl/src/m/classes/qmu/linear_equality_constraint.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/linear_equality_constraint.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/linear_equality_constraint.py (revision 23095) @@ -0,0 +1,165 @@ +import numpy as np +from lclist_write import * +from MatlabArray import * + +class linear_equality_constraint: + ''' + constructor for the linear_equality_constraint class. + + [lec] = linear_equality_constraint.linear_equality_constraint(args) + lec = linear_equality_constraint() + + where the required args are: + matrix (double row, variable coefficients, float('NaN')) + target (double vector, target values, 0.) + and the optional args and defaults are: + scale_type (char, scaling type, 'none') + scale (double, scaling factor, 1.) + + note that zero arguments constructs a default instance, one + argument of the class copies the instance, and two or more + arguments constructs a new instance from the arguments. +''' + def __init__(self): + self.matrix = np.array([[float('NaN')]]) + self.target = 0. + self.scale_type = 'none' + self.scale = 1. + + @staticmethod + def linear_equality_constraint(*args): + nargin = len(args) + + # create a default object + if nargin == 0: + return linear_equality_constraint() + + # copy the object + elif nargin == 1: + if isinstance(args[0],linear_equality_constraint): + lec = args[0] + else: + raise RuntimeError('Object is a '+str(type(args[0]))+' class object, not "linear_equality_constraint"') + + # create the object from the input + else: + if (np.shape(args[0],1) == array_np.size(args[1:min(nargin,4)]) or np.shape(args[0],1) == 1): + asizec = np.shape(args[1:min(nargin,4)]) + elif (array_np.size(args[1:min(nargin,4)]) == 1): + asizec = [np.shape(args[0],1),1] + else: + raise RuntimeError('Matrix for object of class '+str(type(lec))+' has inconsistent number of rows.') + + lec = [linear_equality_constraint() for i in range(asizec[0]) for j in range(asizec[1])] + + for i in range(np.size(lec)): + if (np.shape(args[0])[0] > 1): + lec[i].matrix = args[0][i,:] + else: + lec[i].matrix = args[0] + + if (nargin >= 2): + for i in range(np.size(lec)): + if (np.size(args[1]) > 1): + lec[i].target = args[1][i] + else: + lec[i].target = args[1] + + if (nargin >= 3): + for i in range(np.size(lec)): + if (np.size(args[2]) > 1): + lec[i].scale_type = args[2][i] + else: + lec[i].scale_type = str(args[2]) + + if (nargin >= 4): + for i in range(np.size(lec)): + if (np.size(args[3]) > 1): + lec[i].scale = args[3][i] + else: + lec[i].scale = args[3] + + if (nargin > 4): + print 'WARNING: linear_equality_constraint:extra_arg: Extra arguments for object of class '+str(type(lec))+'.' + + return lec + + + def __repr__(self): + # display the object + string = '\n' + string += 'class "linear_equality_constraint" object = \n' + string += ' matrix: ' +str(self.matrix) + '\n' + string += ' target: ' +str(self.target) + '\n' + string += ' scale_type: ' +str(self.scale_type) + '\n' + string += ' scale: ' +str(self.scale) + '\n' + + return string + + @staticmethod + def prop_matrix(lec): + if type(lec) not in [list,np.ndarray]: + return lec.matrix + + matrix = np.zeros(np.size(lec)) + for i in range(np.size(lec)): + matrix[i,0:np.shape(lec[i].matrix)[1]] = lec[i].matrix[0,:] + + return matrix + + @staticmethod + def prop_lower(lec): + lower=[] + return lower + + @staticmethod + def prop_upper(lec): + upper=[] + return upper + + @staticmethod + def prop_target(lec): + if type(lec) not in [list,np.ndarray]: + return lec.target + + target = np.zeros(np.shape(lec)) + for i in range(np.size(lec)): + target[i] = lec[i].target + + target = allequal(target,0.) + + return target + + @staticmethod + def prop_stype(lec): + if type(lec) not in [list,np.ndarray]: + return lec.scale_type + + stype = ['' for i in range(np.size(lec))] + for i in range(np.size(lec)): + stype[i] = str(lec[i].scale_type) + + stype = allequal(stype,'none') + + return stype + + @staticmethod + def prop_scale(lec): + if type(lec) not in [list,np.ndarray]: + return lec.scale + + scale = np.zeros(np.shape(lec)) + for i in range(np.size(lec)): + scale[i] = lec[i].scale + + scale = allequal(scale,1.) + + return scale + + @staticmethod + def dakota_write(fidi,dvar): + # collect only the variables of the appropriate class + lec = [struc_type(i,'linear_equality_constraint','lec') for i in dvar] + + # write constraints + lclist_write(fidi,'linear_equality_constraints','linear_equality',lec) Index: /issm/trunk-jpl/src/m/classes/qmu/linear_inequality_constraint.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/linear_inequality_constraint.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/linear_inequality_constraint.py (revision 23095) @@ -0,0 +1,187 @@ +import numpy as np +from lclist_write import * +from MatlabArray import * + +class linear_inequality_constraint: + ''' + constructor for the linear_inequality_constraint class. + + [lic] = linear_inequality_constraint.linear_inequality_constraint(args) + lic = linear_inequality_constraint() + + where the required args are: + matrix (double row, variable coefficients, float('NaN')) + lower (double vector, lower bounds, -np.Inf) + upper (double vector, upper bounds, 0.) + and the optional args and defaults are: + scale_type (char, scaling type, 'none') + scale (double, scaling factor, 1.) + + note that zero arguments constructs a default instance, one + argument of the class copies the instance, and three or more + arguments constructs a new instance from the arguments. +''' + def __init__(self): + self.matrix = np.array([[float('NaN')]]) + self.lower = -np.Inf + self.upper = 0. + self.scale_type = 'none' + self.scale = 1. + + @staticmethod + def linear_inequality_constraint(*args): + nargin = len(args) + + # create a default object + if nargin == 0: + return linear_inequality_constraint() + + # copy the object + if nargin == 1: + if isinstance(args[0],linear_inequality_constraint): + lic = args[0] + else: + raise RuntimeError('Object is a '+str(type(args[0]))+' class object, not "linear_inequality_constraint".') + + # not enough arguments + if nargin == 2: + raise RuntimeError('Construction of linear_inequality_constraint class object requires at least 3 inputs.') + + # create the object from the input + else: + if (np.shape(args[0],1) == array_numel(args[1:min(nargin,5)]) or np.shape(args[0],1) == 1): + asizec = array_size(args[1:min(nargin,5)]) + elif (array_numel(args[1:min(nargin,5)]) == 1): + asizec = [array_size(args[0],1),1] + else: + raise RuntimeError('Matrix for object of class '+str(type(lic))+' has inconsistent number of rows.') + + lic = [linear_inequality_constraint() for i in range(asizec[0]) for j in range(asizec[1])] + + for i in range(np.size(lic)): + if (np.shape(args[0],1) > 1): + lic[i].matrix = args[0][i,:] + else: + lic[i].matrix = args[0] + + if (nargin >= 2): + for i in range(np.size(lic)): + if (np.size(args[1]) > 1): + lic[i].lower = args[1][i] + else: + lic[i].lower = args[1] + + if (nargin >= 3): + for i in range(np.size(lic)): + if (np.size(args[2]) > 1): + lic[i].upper = args[2][i] + else: + lic[i].upper = args[2] + + if (nargin >= 4): + for i in range(np.size(lic)): + if (np.size(args[3]) > 1): + lic[i].scale_type = args[3][i] + else: + lic[i].scale_type = str(args[3]) + + if (nargin >= 5): + for i in range(np.size(lic)): + if (np.size(args[4]) > 1): + lic[i].scale = args[4][i] + else: + lic[i].scale = args[4] + + if (nargin > 5): + print 'WARNING: linear_inequality_constraint:extra_arg: Extra arguments for object of class '+str(type(lic))+'.' + + return lic + + + def __repr__(self): + # display the object + string = '\n' + string += 'class "linear_inequality_constraint" object = \n' + string += ' matrix: ' +str(string_vec(self.matrix)) + '\n' + string += ' lower: ' +str(self.lower) + '\n' + string += ' upper: ' +str(self.upper) + '\n' + string += ' scale_type: ' +str(self.scale_type) + '\n' + string += ' scale: ' +str(self.scale) + '\n' + + return string + + @staticmethod + def prop_matrix(lic): + if type(lic) not in [list,np.ndarray]: + return lic.matrix + + matrix = np.zeros(np.size(lic)) + for i in range(np.size(lic)): + matrix[i,0:np.shape(lic[i].matrix)[1]] = lic[i].matrix[0,:] + + return matrix + + @staticmethod + def prop_lower(lic): + if type(lic) not in [list,np.ndarray]: + return lic.lower + + lower = np.zeros(np.shape(lic)) + for i in range(np.size(lic)): + lower[i] = lic[i].lower + + lower = allequal(lower,-np.Inf) + + return lower + + @staticmethod + def prop_upper(lic): + if type(lic) not in [list,np.ndarray]: + return lic.upper + + upper = np.zeros(np.shape(lic)) + for i in range(np.size(lic)): + upper[i] = lic[i].upper + + upper = allequal(upper,0.) + + return upper + + @staticmethod + def prop_target(lic): + target=[] + return target + + @staticmethod + def prop_stype(lic): + if type(lic) not in [list,np.ndarray]: + return lic.scale_type + + stype = ['' for i in range(np.size(lic))] + for i in range(np.size(lic)): + stype[i] = str(lic[i].scale_type) + + stype = allequal(stype,'none') + + return stype + + @staticmethod + def prop_scale(lic): + if type(lic) not in [list,np.ndarray]: + return lic.scale + + scale = np.zeros(np.shape(lic)) + for i in range(np.size(lic)): + scale[i] = lic[i].scale + + scale = allequal(scale,1.) + + return scale + + @staticmethod + def dakota_write(fidi,dvar): + # collect only the variables of the appropriate class + lic = [struc_class(i,'linear_inequality_constraint','lic') for i in dvar] + + # write constraints + lclist_write(fidi,'linear_inequality_constraints','linear_inequality',lic) Index: /issm/trunk-jpl/src/m/classes/qmu/nonlinear_equality_constraint.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/nonlinear_equality_constraint.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/nonlinear_equality_constraint.py (revision 23095) @@ -0,0 +1,179 @@ +import numpy as np +from rlist_write import * +from MatlabArray import * + +class nonlinear_equality_constraint: + ''' + constructor for the nonlinear_equality_constraint class. + + [nec] = nonlinear_equality_constraint.nonlinear_equality_constraint(args) + nec = nonlinear_equality_constraint() + + where the required args are: + descriptor (char, description, '') + target (double vector, target values, 0.) + and the optional args and defaults are: + scale_type (char, scaling type, 'none') + scale (double, scaling factor, 1.) + + note that zero arguments constructs a default instance, one + argument of the class copies the instance, and two or more + arguments constructs a new instance from the arguments. +''' + def __init__(self): + self.descriptor = '' + self.target = 0. + self.scale_type = 'none' + self.scale = 1. + + @staticmethod + def nonlinear_equality_constraint(*args): + nargin = len(args) + + # create a default object + if nargin == 0: + return nonlinear_equality_constraint() + + # copy the object + elif nargin == 1: + if isinstance(args[0],nonlinear_equality_constraint): + nec = args[0] + else: + raise RuntimeError('Object is a '+str(type(args[0]))+' class object, not "nonlinear_equality_constraint"') + + # create the object from the input + else: + asizec = array_size(*args[0:min(nargin,4)]) + nec = [nonlinear_equality_constraint() for i in range(asizec[0]) for j in range(asizec[1])] + + for i in range(np.size(nec)): + if (np.shape(args[0])[0] > 1): + nec[i].descriptor = args[0][i] + else: + nec[i].descriptor = str(args[0]) + if (np.size(args[1]) > 1): + nec[i].target = args[1][i] + else: + nec[i].target = args[1] + + if (nargin >= 3): + for i in range(np.size(nec)): + if (np.size(args[2]) > 1): + nec[i].scale_type = args[2][i] + else: + nec[i].scale_type = str(args[2]) + + if (nargin >= 4): + for i in range(np.size(nec)): + if (np.size(args[3]) > 1): + nec[i].scale =args[3][i] + else: + nec[i].scale =args[3] + + if (nargin > 4): + print 'WARNING: nonlinear_equality_constraint:extra_arg: Extra arguments for object of class '+str(type(nec))+'.' + + return nec + + + def __repr__(self): + # display the object + string = '\n' + string += 'class "nonlinear_equality_constraint" object = \n' + string += ' descriptor: ' +str(self.descriptor) + '\n' + string += ' target: ' +str(self.target) + '\n' + string += ' scale_type: ' +str(self.scale_type) + '\n' + string += ' scale: ' +str(self.scale) + '\n' + + return string + + @staticmethod + def prop_desc(nec,dstr): + if type(nec) not in [list,np.ndarray]: + if nec.descriptor != '' or type(nec.descriptor) != str: + desc = str(nec.descriptor) + elif dstr != '': + desc = str(dstr) + else: + desc = 'nec' + return desc + + desc=['' for i in range(np.size(nec))] + for i in range(np.size(nec)): + if nec[i].descriptor != '' or type(nec[i].descriptor) != str: + desc[i] = str(nec[i].descriptor) + elif dstr != '': + desc[i] = str(dstr)+str(string_dim(nec,i,'vector')) + else: + desc[i] = 'nec'+str(string_dim(nec,i,'vector')) + + desc=allempty(desc) + + return desc + + @staticmethod + def prop_lower(nec): + lower=[] + return lower + + @staticmethod + def prop_upper(nec): + upper=[] + return upper + + @staticmethod + def prop_weight(nec): + weight=[] + return weight + + @staticmethod + def prop_target(nec): + if type(nec) not in [list,np.ndarray]: + return nec.target + + target = np.zeros(np.shape(nec)) + for i in range(np.size(nec)): + target[i] = nec[i].target + + target = allequal(target,0.) + + return target + + @staticmethod + def prop_stype(nec): + if type(nec) not in [list,np.ndarray]: + return nec.scale_type + + stype = ['' for i in range(np.size(nec))] + for i in range(np.size(nec)): + stype[i] = str(nec[i].scale_type) + + stype = allequal(stype,'none') + + return stype + + @staticmethod + def prop_scale(nec): + if type(nec) not in [list,np.ndarray]: + return nec.scale + + scale = np.zeros(np.shape(nec)) + for i in range(np.size(nec)): + scale[i] = nec[i].scale + + scale = allequal(scale,1.) + + return scale + + @staticmethod + def dakota_write(fidi,dresp,rdesc): + # colnect only the variables of the appropriate class + nec = [struc_type(i,'nonlinear_equality_constraint','nec') for i in dresp] + + # write constraints + rdesc = rlist_write(fidi,'nonlinear_equality_constraints','nonlinear_equality',nec,rdesc) + return rdesc + + @staticmethod + def dakota_rlev_write(fidi,dresp,params): + return Index: /issm/trunk-jpl/src/m/classes/qmu/nonlinear_inequality_constraint.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/nonlinear_inequality_constraint.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/nonlinear_inequality_constraint.py (revision 23095) @@ -0,0 +1,198 @@ +import numpy as np +from rlist_write import * +from MatlabArray import * + +class nonlinear_inequality_constraint: + ''' + constructor for the nonlinear_inequality_constraint class. + + [nic] = nonlinear_inequality_constraint.nonlinear_inequality_constraint(args) + nic = nonlinear_inequality_constraint() + + where the required args are: + descriptor (char, description, '') + lower (double, lower bound, -np.inf) + upper (double, upper bound, 0.) + and the optional args and defaults are: + scale_type (char, scaling type, 'none') + scale (double, scaling factor, 1.) + + note that zero arguments constructs a default instance, one + argument of the class copies the instance, and three or more + arguments constructs a new instance from the arguments. +''' + def __init__(self): + self.descriptor = '' + self.lower = -np.inf + self.upper = 0. + self.scale_type = 'none' + self.scale = 1. + + @staticmethod + def nonlinear_inequality_constraint(*args): + nargin = len(args) + + # create a default object + if nargin == 0: + return nonlinear_inequality_constraint() + + # copy the object + if nargin == 1: + if isinstance(args[0],nonlinear_inequality_constraint): + nic = args[0] + else: + raise RuntimeError('Object is a '+str(type(args[0]))+' class object, not "nonlinear_inequality_constraint".') + + # not enough arguments + if nargin == 2: + raise RuntimeError('Construction of nonlinear_inequality_constraint class object requires at least 3 inputs.') + + # create the object from the input + else: + asizec = array_size(*args[0:min(nargin,3)]) + nic = [nonlinear_inequality_constraint() for i in range(asizec[0]) for j in range(asizec[1])] + + for i in range(np.size(nic)): + if (np.size(args[0]) > 1): + nic[i].descriptor = args[0][i]; + else: + nic[i].descriptor = str(args[0])+string_dim(nic,i,'vector') + if (np.size(args[1]) > 1): + nic[i].lower = args[1][i] + else: + nic[i].lower = args[1] + + if (np.size(args[2]) > 1): + nic[i].upper = args[2][i] + else: + nic[i].upper = args[2] + + if (nargin >= 4): + for i in range(np.size(nic)): + if (np.size(args[3]) > 1): + nic[i].scale_type = args[3][i] + else: + nic[i].scale_type = str(args[3]) + + if (nargin >= 5): + for i in range(np.size(nic)): + if (np.size(args[4]) > 1): + nic[i].upper = args[4][i] + else: + nic[i].upper = args[4] + + if (nargin > 5): + print 'WARNING: nonlinear_inequality_constraint:extra_arg: Extra arguments for object of class '+str(type(nic))+'.' + + return nic + + def __repr__(self): + # display the object + string = '\n' + string += 'class "nonlinear_inequality_constraint" object = \n' + string += ' descriptor: ' +str(self.descriptor) + '\n' + string += ' lower: ' +str(self.lower) + '\n' + string += ' upper: ' +str(self.upper) + '\n' + string += ' scale_type: ' +str(self.scale_type) + '\n' + string += ' scale: ' +str(self.scale) + '\n' + + return string + + @staticmethod + def prop_desc(nic,dstr): + if type(nic) not in [list,np.ndarray]: + if nic.descriptor != '' or type(nic.descriptor) != str: + desc = str(nic.descriptor) + elif dstr != '': + desc = str(dstr) + else: + desc = 'nic' + return desc + + desc = ['' for i in range(np.size(nic))] + for i in range(np.size(nic)): + if nic[i].descriptor != '' or type(nic[i].descriptor) != str: + desc[i] = str(nic[i].descriptor) + elif dstr != '': + desc[i] = str(dstr)+str(string_dim(nic,i,'vector')) + else: + desc[i] = 'nic'+str(string_dim(nic,i,'vector')) + + desc = allempty(desc) + + return desc + + @staticmethod + def prop_stype(nic): + if type(nic) not in [list,np.ndarray]: + return nic.scale_type + + stype = ['' for i in range(np.size(nic))] + for i in range(np.size(nic)): + stype[i] = str(nic[i].scale_type) + + stype = allequal(stype,'none') + + return stype + + @staticmethod + def prop_scale(nic): + if type(nic) not in [list,np.ndarray]: + return nic.scale + + scale = np.zeros(np.shape(nic)) + for i in range(np.size(nic)): + scale[i] = nic[i].scale + + scale = allequal(scale,1.) + + return scale + + @staticmethod + def prop_weight(nic): + weight=[] + return weight + + @staticmethod + def prop_lower(nic): + if type(nic) not in [list,np.ndarray]: + return nic.lower + + lower = np.zeros(np.shape(nic)) + for i in range(np.size(nic)): + lower[i] = nic[i].lower + + lower = allequal(lower,-np.inf) + + return lower + + @staticmethod + def prop_upper(nic): + if type(nic) not in [list,np.ndarray]: + return nic.upper + + upper = np.zeros(np.shape(nic)) + for i in range(np.size(nic)): + upper[i] = nic[i].upper + + upper = allequal(upper,0.) + + return upper + + @staticmethod + def prop_target(nic): + target=[] + return target + + @staticmethod + def dakota_write(fidi,dresp,rdesc): + # collect only the variables of the appropriate class + nic = [struc_class(i,'nonlinear_inequality_constraint','nic') for i in dresp] + + # write constraints + rdesc = rlist_write(fidi,'nonlinear_inequality_constrants','nonlinear_inequality',nic,rdesc) + return rdesc + + @staticmethod + def dakota_rlev_write(fidi,dresp,params): + return Index: /issm/trunk-jpl/src/m/classes/qmu/normal_uncertain.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/normal_uncertain.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/normal_uncertain.py (revision 23095) @@ -0,0 +1,181 @@ +import numpy as np +from vlist_write import * +from MatlabArray import * + +class normal_uncertain(object): + ''' + definition for the normal_uncertain class. + + [nuv] = normal_uncertain.normal_uncertain(args) + nuv = normal_uncertain() + + where the required args are: + descriptor (str, description, '') + mean (float, mean, float('NaN')) + stddev (float, standard deviation, float('NaN')) + and the optional args and defaults are: + lower (float, lower bound, -np.Inf) + upper (float, upper bound, np.Inf) + + note that zero arguments constructs a default instance, one + argument of the class copies the instance, and three or more + arguments constructs a new instance from the arguments. +''' + def __init__(self): + self.descriptor = '' + self.mean = float('NaN') + self.stddev = float('NaN') + self.lower =-np.Inf + self.upper = np.Inf + + @staticmethod + def normal_uncertain(*args): + nargin = len(args) + + # create a default object + if nargin == 0: + return normal_uncertain() + + # copy the object + elif nargin == 1: + if isinstance(args[0],normal_uncertain): + nuv = args[0] + else: + raise RuntimeError('Object '+str(args[0])+' is a '+str(type(args[0]))+' class object, not "normal_uncertain".') + + # not enough arguments + elif nargin == 2: + raise RuntimeError('Construction of "normal_uncertain" class object requires at least 3 inputs.') + + # create the object from the input + else: + # lines differ here in other classes/tests; see asizec problem in notes + nuv=normal_uncertain() + + nuv.descriptor = str(args[0]) + nuv.mean = args[1] + nuv.stddev = args[2] + if nargin >= 4: + nuv.lower = args[3] + if nargin >= 5: + nuv.upper = args[4] + if nargin > 5: + print 'WARNING: normal_uncertain:extra_arg: Extra arguments for object of class '+str(type(nuv))+'.' + + return [nuv] + + def __repr__(self): + # display an individual object + string = '\n' + string += 'class "normal_uncertain" object = \n' + string += ' descriptor: '+str(self.descriptor) + '\n' + string += ' mean: '+str(self.mean) + '\n' + string += ' stddev: '+str(self.stddev) + '\n' + string += ' lower: '+str(self.lower) + '\n' + string += ' upper: '+str(self.upper) + '\n' + + return string + + # from here on, nuv is either a single, or a 1d vector of, normal_uncertain + + @staticmethod + def prop_desc(nuv,dstr): + if type(nuv) not in [list,np.ndarray]: + if nuv.descriptor != '' or type(nuv.descriptor) != str: + desc = str(nuv.descriptor) + elif dstr != '': + desc = str(dstr) + else: + desc = 'nuv' + return desc + + desc = ['' for i in range(np.size(nuv))] + for i in range(np.size(nuv)): + if nuv[i].descriptor != '' or type(nuv[i].descriptor) != str: + desc[i] = str(nuv[i].descriptor) + elif dstr != '': + desc[i] = str(dstr)+str(string_dim(nuv,i,'vector')) + else: + desc[i] = 'nuv'+str(string_dim(nuv,i,'vector')) + + desc = allempty(desc) + + return desc + + @staticmethod + def prop_initpt(nuv): + initpt=[] + return initpt + + @staticmethod + def prop_lower(nuv): + if type(nuv) not in [list,np.ndarray]: + return nuv.lower + + lower = np.zeros(np.size(nuv)) + for i in range(np.size(nuv)): + lower[i] = nuv[i].lower + + lower = allequal(lower,-np.inf) + + return lower + + @staticmethod + def prop_upper(nuv): + if type(nuv) not in [list,np.ndarray]: + return nuv.upper + + upper = np.zeros(np.size(nuv)) + for i in range(np.size(nuv)): + upper[i] = nuv[i].upper + + upper = allequal(upper,-np.inf) + + return upper + + @staticmethod + def prop_mean(nuv): + if type(nuv) not in [list,np.ndarray]: + return nuv.mean + + mean = np.zeros(np.size(nuv)) + for i in range(np.size(nuv)): + mean[i] = nuv[i].mean + + return mean + + @staticmethod + def prop_stddev(nuv): + if type(nuv) not in [list,np.ndarray]: + return nuv.stddev + + stddev = np.zeros(np.size(nuv)) + for i in range(np.size(nuv)): + stddev[i] = nuv[i].stddev + + return stddev + + @staticmethod + def prop_initst(nuv): + initst=[] + return initst + + @staticmethod + def prop_stype(nuv): + stype=[] + return stype + + @staticmethod + def prop_scale(nuv): + scale=[] + return scale + + @staticmethod + def dakota_write(fidi,dvar): + # collect only the variables of the appropriate class + nuv = [struc_class(i,'normal_uncertain','nuv') for i in dvar] + + # possible namespace pollution, the above import seems not to work + from vlist_write import * + # write variables + vlist_write(fidi,'normal_uncertain','nuv',nuv) Index: /issm/trunk-jpl/src/m/classes/qmu/objective_function.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/objective_function.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/objective_function.py (revision 23095) @@ -0,0 +1,178 @@ +import numpy as np +from rlist_write import * +from MatlabArray import * + +class objective_function(object): + ''' + definition for the objective_function class. + + [of] = objective_function.objective_function(args) + of = objective_function() + + where the required args are: + descriptor (char, description, '') + and the optional args and defaults are: + scale_type (char, scaling type, 'none') + scale (double, scaling factor, 1.) + weight (double, weighting factor, 1.) + + note that zero arguments constructs a default instance, one + argument of the class copies the instance, and one or more + arguments constructs a new instance from the arguments. +''' + def __init__(self): + self.descriptor = '' + self.scale_type = 'none' + self.scale = 1. + self.weight = 1. + + @staticmethod + def objective_function(*args): + nargin = len(args) + + # create a default object + if nargin == 0: + return objective_function() + + # copy the object or create the object from the input + else: + if (nargin == 1) and isinstance(args[0],objective_function): + of = args[0] + else: + shapec = array_size(*args[0:min(nargin,4)]) + of = [objective_function() for i in range(shapec[0]) for j in range(shapec[1])] + + for i in range(np.size(of)): + if (np.size(args[0]) > 1): + of[i].descriptor = args[0][i] + else: + of[i].descriptor = str(args[0])+string_dim(of,i,'vector') + + if (nargin >= 2): + for i in range(np.size(of)): + if (np.size(args[1]) > 1): + of[i].scale_type = args[1][i] + else: + of[i].scale_type = str(args[1]) + + if (nargin >= 3): + for i in range(np.size(of)): + if (np.size(args[2]) > 1): + of[i].scale = args[2][i] + else: + of[i].scale = args[2] + + if (nargin >= 4): + for i in range(np.size(of)): + if (np.size(args[3]) > 1): + of[i].weight = args[3][i] + else: + of[i].weight = args[3] + + if (nargin > 4): + print 'WARNING: objective_function:extra_arg Extra arguments for object of class '+str(type(of))+'.' + + return of + + + def __repr__(self): + # display the object + string = '\n' + string += 'class "objective_function" object = \n' + string += ' descriptor: ' +str(self.descriptor) + '\n' + string += ' scale_type: ' +str(self.scale_type) + '\n' + string += ' scale: ' +str(self.scale) + '\n' + string += ' weight: ' +str(self.weight) + '\n' + + return string + + @staticmethod + def prop_desc(of,dstr): + if type(of) not in [list,np.ndarray]: + if of.descriptor != '' or type(of.descriptor) != str: + desc = str(of.descriptor) + elif dstr != '': + desc = str(dstr) + else: + desc = 'of' + return desc + + desc = ['' for i in range(np.size(of))] + for i in range(np.size(of)): + if of[i].descriptor != '' or type(of[i].descriptor) != str: + desc[i] = str(of[i].descriptor) + elif dstr != '': + desc[i] = str(dstr)+str(string_dim(of,i,'vector')) + else: + desc[i] = 'of'+str(string_dim(of,i,'vector')) + + desc = allempty(desc) + + return desc + + @staticmethod + def prop_lower(of): + lower=[] + return lower + + @staticmethod + def prop_upper(of): + upper=[] + return upper + + @staticmethod + def prop_target(of): + target=[] + return target + + @staticmethod + def prop_weight(of): + if type(of) not in [list,np.ndarray]: + return of.weight + + weight = np.zeros(np.shape(of)) + for i in range(np.size(of)): + weight[i] = of[i].weight + + weight = allequal(weight,1.) + + return weight + + @staticmethod + def prop_stype(of): + if type(of) not in [list,np.ndarray]: + return of.scale_type + + stype = ['' for i in range(np.size(of))] + for i in range(np.size(of)): + stype[i] = str(of[i].scale_type) + + stype = allequal(stype,'none') + + return stype + + @staticmethod + def prop_scale(of): + if type(of) not in [list,np.ndarray]: + return of.scale + + scale = np.zeros(np.shape(of)) + for i in range(np.size(of)): + scale[i] = of[i].scale + + scale = allequal(scale,1.) + + return scale + + @staticmethod + def dakota_write(fidi,dresp,rdesc): + # coloft only the variables of the appropriate class + of = [struc_class(i,'objective_functions','of') for i in dresp] + + # write constraints + rdesc = rlist_write(fidi,'objective_functions','objective_function',of,rdesc) + return rdesc + + @staticmethod + def dakota_rlev_write(fidi,dresp,params): + return Index: /issm/trunk-jpl/src/m/classes/qmu/qmu_classes.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/qmu_classes.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/qmu_classes.py (revision 23095) @@ -0,0 +1,30 @@ +# nuv +from normal_uncertain import * +# uuv +from uniform_uncertain import * + +# lic +from linear_inequality_constraint import * +# lec +from linear_equality_constraint import * + +#nic +from nonlinear_inequality_constraint import * +#nec +from nonlinear_equality_constraint import * + +# csv +from continuous_design import * +# cdv +from continuous_state import * + +# lst +from least_squares_term import * + +# of +from objective_function import * +# rf +from response_function import * + +# cf +from calibration_function import * Index: /issm/trunk-jpl/src/m/classes/qmu/response_function.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/response_function.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/response_function.py (revision 23095) @@ -0,0 +1,191 @@ +import numpy as np +from rlist_write import * +from rlev_write import * +from MatlabArray import * + +#move this later +from helpers import * + +class response_function(object): + ''' + definition for the response_function class. + + [rf] = response_function.response_function(args) + rf = response_function() + + where the required args are: + descriptor (char, description, '') + and the optional args and defaults are: + respl (double vector, response levels, []) + probl (double vector, probability levels, []) + rell (double vector, reliability levels, []) + grell (double vector, gen. reliability levels, []) + + note that zero arguments constructs a default instance, one + argument of the class copies the instance, and one or more + arguments constructs a new instance from the arguments. +''' + + def __init__(self): + self.descriptor = '' + self.respl = [] + self.probl = [] + self.rell = [] + self.grell = [] + + @staticmethod + def response_function(*args): + + nargin = len(args) + # create a default object + if nargin == 0: + return response_function() + + # copy the object or create the object from the input + else: + if nargin == 1 and isinstance(args[0],response_function): + rf = args[0] + else: + asizec = array_size(*args[0:min(nargin,1)]) + rf = [response_function() for i in range(asizec[0]) for j in range(asizec[1])] + + for i in range(np.size(rf)): + if (np.size(args[0]) > 1): + rf[i].descriptor = args[0][i] + else: + rf[i].descriptor = str(args[0])+string_dim(rf,i,'vector') + + if nargin >= 2: + for i in range(np.size(rf)): + rf[i].respl = args[1] + + if nargin >= 3: + for i in range(np.size(rf)): + rf[i].probl = args[2] + + if nargin >= 4: + for i in range(np.size(rf)): + rf[i].rell = args[3] + + if nargin >= 5: + for i in range(np.size(rf)): + rf[i].grell = args[4] + + if nargin > 5: + print 'WARNING: response_function:extra_arg: Extra arguments for object of class '+str(type(rf))+'.' + + return rf + + + def __repr__(self): + # display the object + string = '\n' + string += 'class "response_function" object = \n' + string += ' descriptor: ' +str(self.descriptor) + '\n' + string += ' respl: ' +str(self.respl) + '\n' + string += ' probl: ' +str(self.probl) + '\n' + string += ' rell: ' +str(self.rell) + '\n' + string += ' grell: ' +str(self.grell) + '\n' + + return string + + def __len__(self): + return max(len(self.respl),len(self.probl),len(self.rell),len(self.grell)) + + # from here on, rf is either a single, or a 1d vector of, response_function + + @staticmethod + def prop_desc(rf,dstr): + # response_function is always a vector, or should be, even with just 1 + if type(rf) not in [list,np.ndarray]: + rf = [rf] + + desc = ['' for i in range(np.size(rf))] + for i in range(np.size(rf)): + if rf[i].descriptor != '' or type(rf[i].descriptor) != str: + desc[i] = str(rf[i].descriptor) + elif dstr != '': + desc[i] = str(dstr)+str(string_dim(rf,i,'vector')) + else: + desc[i] = 'rf'+str(string_dim(rf,i,'vector')) + + desc = allempty(desc) + + return desc + + @staticmethod + def prop_stype(rf): + stype=[] + return stype + + @staticmethod + def prop_scale(rf): + scale=[] + return scale + + @staticmethod + def prop_weight(rf): + weight=[] + return weight + + @staticmethod + def prop_lower(rf): + lower=[] + return lower + + @staticmethod + def prop_upper(rf): + upper=[] + return upper + + @staticmethod + def prop_target(rf): + target=[] + return target + + @staticmethod + def prop_levels(rf): + # response_function is always a vector, or should be, even with just 1 + if type(rf) not in [list,np.ndarray]: + rf = [rf] + + respl = empty_nd_list(np.size(rf)) + + probl = empty_nd_list(np.size(rf)) + + rell = empty_nd_list(np.size(rf)) + + grell = empty_nd_list(np.size(rf)) + + for i in range(np.size(rf)): + respl[i] = rf[i].respl + probl[i] = rf[i].probl + rell [i] = rf[i].rell + grell[i] = rf[i].grell + + respl = allempty(respl) + probl = allempty(probl) + rell = allempty(rell) + grell = allempty(grell) + + return [respl,probl,rell,grell] + + @staticmethod + def dakota_write(fidi,dresp,rdesc): + # collect only the responses of the appropriate class + rf = [struc_class(vars(dresp)[i][j],'response_function','rf') for i in fieldnames(dresp) for j in range(len(vars(dresp)[i]))] + + #possible namespace pollution here + from rlist_write import rlist_write + # write responses + rdesc = rlist_write(fidi,'response_function','rf',rf,rdesc) + + return rdesc + + @staticmethod + def dakota_rlev_write(fidi,dresp,params): + # collect only the responses of the appropriate class + rf = [struc_class(vars(dresp)[i][j],'response_function','rf') for i in fieldnames(dresp) for j in range(len(vars(dresp)[i]))] + + # write response levels + rlev_write(fidi,rf,'response_function',params) Index: /issm/trunk-jpl/src/m/classes/qmu/uniform_uncertain.py =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu/uniform_uncertain.py (revision 23095) +++ /issm/trunk-jpl/src/m/classes/qmu/uniform_uncertain.py (revision 23095) @@ -0,0 +1,164 @@ +import numpy as np +from vlist_write import * +from MatlabArray import * + +class uniform_uncertain(object): + ''' + definition for the uniform_uncertain class. + + [uuv] = uniform_uncertain.uniform_uncertain(args) + uuv = uniform_uncertain() + + where the required args are: + descriptor (str, description, '') + lower (float, lower bound, -np.Inf) + upper (float, upper bound, np.Inf) + + note that zero arguments constructs a default instance, one + argument of the class copies the instance, and three or more + arguments constructs a new instance from the arguments. +''' + + def __init__(self): + self.descriptor = '' + self.lower = -np.Inf + self.upper = np.Inf + + @staticmethod + def uniform_uncertain(*args): + nargin = len(args) + + # create a default object + if nargin == 0: + return uniform_uncertain() + + # copy the object + elif nargin == 1: + if isinstance(args[0],uniform_uncertain): + uuv = args[0] + else: + raise RuntimeError('Object '+str(args[0])+' is a '+str(type(args[0]))+' class object, not "uniform_uncertain".') + + # not enough arguments + + elif nargin == 2: + raise RuntimeError('Construction of "uniform_uncertain" class object requires at least 3 inputs.') + + # create the object from the input + else: + # leaving this here in case it becomes important in the future + #asizec=array_size(*args[0:min(nargin,3)]) + #uuv = [uniform_uncertain() for i in range(asizec[0]) for j in range(asizec[1])] + + uuv = uniform_uncertain() + + uuv.descriptor = str(args[0]) + uuv.lower = args[1] + uuv.upper = args[2] + if (nargin > 3): + print 'WARNING: uniform_uncertain:extra_arg: Extra arguments for object of class '+type(uuv)+'.' + + return [uuv] + + + def __repr__(self): + # display an individual object + string = '\n' + string += 'class "uniform_uncertain" object = \n' + string += ' descriptor: ' + str(self.descriptor) + '\n' + string += ' lower: ' + str(self.lower) + '\n' + string += ' upper: ' + str(self.upper) + '\n' + + return string + + # from here on, uuv is either a single, or a 1d vector of, uniform_uncertain + + @staticmethod + def prop_desc(uuv,dstr): + if type(uuv) not in [list,np.ndarray]: + if uuv.descriptor != '' or type(uuv.descriptor) != str: + desc = str(uuv.descriptor) + elif dstr != '': + desc = str(dstr) + else: + desc = 'uuv' + return desc + + desc = ['' for i in range(np.size(uuv))] + for i in range(np.size(uuv)): + if uuv[i].descriptor != '' or type(uuv[i].descriptor) != str: + desc[i] = str(uuv[i].descriptor) + elif dstr != '': + desc[i] = str(dstr)+str(string_dim(uuv,i,'vector')) + else: + desc[i] = 'uuv'+str(string_dim(uuv,i,'vector')) + + desc = allempty(desc) + + return desc + + @staticmethod + def prop_initpt(uuv): + initpt=[] + return initpt + + @staticmethod + def prop_lower(uuv): + if type(uuv) not in [list,np.ndarray]: + return uuv.lower + + lower = np.zeros(np.size(uuv)) + for i in range(np.size(uuv)): + lower[i] = uuv[i].lower + + lower = allequal(lower,-np.Inf) + + return lower + + @staticmethod + def prop_upper(uuv): + if type(uuv) not in [list,np.ndarray]: + return uuv.upper + + upper = np.zeros(np.size(uuv)) + for i in range(np.size(uuv)): + upper[i] = uuv[i].upper + + upper = allequal(upper, np.Inf) + + return upper + + @staticmethod + def prop_mean(uuv): + mean=[] + return mean + + @staticmethod + def prop_stddev(uuv): + stddev=[] + return stddev + + @staticmethod + def prop_initst(uuv): + initst=[] + return initst + + @staticmethod + def prop_stype(uuv): + stype=[] + return stype + + @staticmethod + def prop_scale(uuv): + scale=[] + return scale + + @staticmethod + def dakota_write(fidi,dvar): + # collect only the variables of the appropriate class + uuv = [struc_class(i,'uniform_uncertain','uuv') for i in dvar] + + # possible namespace pollution, the above import seems not to work + from vlist_write import * + # write variables + vlist_write(fidi,'uniform_uncertain','uuv',uuv) Index: /issm/trunk-jpl/src/m/dev/ISSM.py =================================================================== --- /issm/trunk-jpl/src/m/dev/ISSM.py (revision 23095) +++ /issm/trunk-jpl/src/m/dev/ISSM.py (revision 23095) @@ -0,0 +1,71 @@ +print 'WARNING: EXPERIMENTAL FEATURE ISSM.py: universal Python ISSM import' + +#Most common imports +import numpy as np +import scipy.io as spio +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from IssmConfing import * + +#Secondary imports +import copy +from scipy.interpolate import interp1d +from operator import itemgetter +from generic import generic +from materials import * +from bamg import * +from SMBgemb import * +from calvingminthickness import * +from calvingvonmises import * +from bamgflowband import * +from paterson import * +from frictionsommers import * +from hydrologysommers import * +from transient import * +from mismipbasalforcings import * +from ComputeHessian import * +from ComputeMetric import * + +# qmu +from dakota_method import * +from qmu_classes import * +from partitioner import * +from dmeth_params_set import * +from dmeth_params_write import * + +#Helper functions +def python_help(): + '''Prints out key code fragments that may be useful to users''' + print 'Differences between Python and Matlab code:' + #... + +def find(to_find): + '''analagous to matlab's find function but requires separate and/or functions''' + return np.array(np.where(to_find)) + +def find_and(*args): + '''analagous to matlab's a & b functionality when used in conjunction with find(), + returns overlap across a and b + takes an arbitrary number of arguments of similar shape''' + result = args[0] + for arg in args[1:]: + if type(arg) != np.ndarray: + arg = np.array(arg) + result = np.intersect1d(result,arg) + return result + +def find_or(*args): + '''analagous to matlab's a | b functionality when used in conjunction with find(), + returns all unique values across a and b + takes an arbitrary number of arguments of similar shape''' + result = args[0] + for arg in args[1:]: + if type(arg) != np.ndarray: + arg = np.array(arg) + result = np.unique(np.concatenate((result,arg))) + return result Index: /issm/trunk-jpl/src/m/dev/devpath.py =================================================================== --- /issm/trunk-jpl/src/m/dev/devpath.py (revision 23094) +++ /issm/trunk-jpl/src/m/dev/devpath.py (revision 23095) @@ -33,6 +33,6 @@ #Manual imports for commonly used functions +from runme import runme #first because plotmodel may fail from plotmodel import plotmodel -from runme import runme #c = get_ipython().config Index: /issm/trunk-jpl/src/m/miscellaneous/normfit_issm.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/normfit_issm.m (revision 23094) +++ /issm/trunk-jpl/src/m/miscellaneous/normfit_issm.m (revision 23095) @@ -58,4 +58,3 @@ end end - end Index: /issm/trunk-jpl/src/m/miscellaneous/normfit_issm.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/normfit_issm.py (revision 23095) +++ /issm/trunk-jpl/src/m/miscellaneous/normfit_issm.py (revision 23095) @@ -0,0 +1,61 @@ +import numpy as np + +from scipy.stats import chi2,t + +# +# wrapper for normfit to avoid using the matlab statistics toolbox. +# +def normfit_issm(x,alpha=None): + + if alpha == None: + alpha=0.05 + + + # check for any NaN in any columns + if not np.isnan(x).any(): + + # explicitly calculate the moments + muhat =np.mean(x,0) + # numpy defaults to 0 delta degrees of freedom; matlab uses 1 + sigmahat=np.std(x,0,ddof=1) + + # no way to ask this in python, assume 4 outputs + #if (nargout>2): + + prob=1.-alpha/2. + + if (np.size(x,0) == 1): + # operate like matlab normfit, mean, std, etc. + n=np.size(x) + else: + n=np.size(x,0) + + muci =np.zeros((2,np.size(muhat ))) + sigmaci =np.zeros((2,np.size(sigmahat))) + + try: + muci[0,:] =muhat-t.ppf(prob,n-1)*sigmahat/np.sqrt(n) + muci[1,:] =muhat+t.ppf(prob,n-1)*sigmahat/np.sqrt(n) + sigmaci[0,:]=sigmahat*np.sqrt((n-1)/chi2.ppf(prob ,n-1)) + sigmaci[1,:]=sigmahat*np.sqrt((n-1)/chi2.ppf(1.-prob,n-1)) + except: + muci[0,:] =muhat + muci[1,:] =muhat + sigmaci[0,:]=sigmahat + sigmaci[1,:]=sigmahat + else: + # must loop over columns, since number of elements could be different + muhat =np.zeros((1,np.size(x,1))) + sigmahat=np.zeros((1,np.size(x,1))) + muci =np.zeros((2,np.size(x,1))) + sigmaci =np.zeros((2,np.size(x,1))) + + # remove any NaN and recursively call column + for j in range(np.shape(x,1)): + [muhat[j],sigmahat[j],muci[:,j],sigmaci[:,j]]=normfit_issm(x[not np.isnan(x[:,j]),j],alpha) + + return [muhat,sigmahat,muci,sigmaci] + + + + Index: /issm/trunk-jpl/src/m/miscellaneous/prctile_issm.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/prctile_issm.m (revision 23094) +++ /issm/trunk-jpl/src/m/miscellaneous/prctile_issm.m (revision 23095) @@ -8,5 +8,4 @@ catch me - if length(size(x)) > 2 error('Number of dimensions %d not implemented.',length(size(x))); @@ -50,5 +49,4 @@ % fill in high and low values - y(pxi(n),:)=repmat(x(n,:),nnz(p>xi(n)),1); @@ -83,4 +81,3 @@ end end - end Index: /issm/trunk-jpl/src/m/miscellaneous/prctile_issm.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/prctile_issm.py (revision 23095) +++ /issm/trunk-jpl/src/m/miscellaneous/prctile_issm.py (revision 23095) @@ -0,0 +1,65 @@ +import numpy as np +from scipy.interpolate import interp1d + +# +# wrapper for prctile to avoid using the matlab statistics toolbox. +# +def prctile_issm(x,p,dim): + + try: + raise RuntimeError('hello world') + + #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 + + 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) + + # 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 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/modules/Chaco.py =================================================================== --- /issm/trunk-jpl/src/m/modules/Chaco.py (revision 23095) +++ /issm/trunk-jpl/src/m/modules/Chaco.py (revision 23095) @@ -0,0 +1,20 @@ +from Chaco_python import Chaco_python + +def Chaco(A,vwgts,ewgts,x,y,z,options,nparts,goal): + '''CHACO + + Usage: + [assgn] = Chaco(A,vwgts,ewgts,x,y,z,options,nparts,goal); + + A: Input adjacency matrix + vwgts: weights for all vertices + ewgts: weights for all edges + x,y,z: coordinates for inertial method + options: architecture and partitioning options + nparts: number of parts options + goal: desired set sizes +''' + # Call mex module + assgn = Chaco_python(A,vwgts,ewgts,x,y,z,options,nparts,goal) + + return assgn Index: /issm/trunk-jpl/src/m/modules/MeshPartition.py =================================================================== --- /issm/trunk-jpl/src/m/modules/MeshPartition.py (revision 23095) +++ /issm/trunk-jpl/src/m/modules/MeshPartition.py (revision 23095) @@ -0,0 +1,42 @@ +from MeshPartition_python import MeshPartition_python + +from mesh3dprisms import * +from mesh2d import * +from mesh2dvertical import * + +def MeshPartition(md,numpartitions): + '''MESHPARTITION - Partition mesh according to the number of areas, using Metis library. + + Usage: + [element_partitioning,node_partitioning]=MeshPartition(md.mesh,numpartitions)") + + element_partitioning: Vector of partitioning area numbers, for every element. + node_partitioning: Vector of partitioning area numbers, for every node. +''' + if md == None or numpartitions == None: + print MeshPartition.__doc__ + raise RuntimeError('Wrong usage (see above)') + + #Get mesh info from md.mesh + numberofvertices = md.mesh.numberofvertices + numberofelements = md.mesh.numberofelements + elements = md.mesh.elements + numberofelements2d = 0 + numberofvertices2d = 0 + numberoflayers = 1 + elements2d = [] + if isinstance(md.mesh,mesh3dprisms): + elementtype = 'Penta' + numberofelements2d = md.mesh.numberofelements2d + numberofvertices2d = md.mesh.numberofvertices2d + numberoflayers = md.mesh.numberoflayers + elements2d = md.mesh.elements2d + elif isinstance(md.mesh,mesh2d): + elementtype = 'Tria' + elif isinstance(md.mesh,mesh2dvertical): + elementtype = 'Tria' + + #Call module + [element_partitioning, node_partitioning] = MeshPartition_python(numberofvertices, elements, numberofvertices2d, elements2d, numberoflayers, meshelementtype, numpartitions) + + return [element_partitioning, node_partitioning] Index: /issm/trunk-jpl/src/m/modules/Scotch.py =================================================================== --- /issm/trunk-jpl/src/m/modules/Scotch.py (revision 23095) +++ /issm/trunk-jpl/src/m/modules/Scotch.py (revision 23095) @@ -0,0 +1,12 @@ +from Scotch_python import Scotch_python + +def Scotch(*varargin): +'''SCOTCH - Scotch partitioner + + Usage: + maptab=Scotch(adjmat,vertlb,vertwt,edgewt,archtyp,archpar,Scotch-specific parameters) +''' + # Call mex module + maptab=Scotch_python(*varargin) + + return maptab Index: /issm/trunk-jpl/src/m/partition/adjacency.py =================================================================== --- /issm/trunk-jpl/src/m/partition/adjacency.py (revision 23095) +++ /issm/trunk-jpl/src/m/partition/adjacency.py (revision 23095) @@ -0,0 +1,33 @@ +import numpy as np +import MatlabFuncs as m +from NodeConnectivity import * +from GetAreas import * + +def adjacency(md): + #ADJACENCY - compute adjacency matrix, list of vertices and list of weights. + # + # function to create the adjacency matrix from the connectivity table. + # + # the required output is: + # md.adj_mat (double [sparse nv x nv], vertex adjacency matrix) + # md.qmu.vertex_weight (double [nv], vertex weights) + + indi=np.array([md.mesh.elements[:,0],md.mesh.elements[:,1],md.mesh.elements[:,2]]) + indj=np.array([md.mesh.elements[:,1],md.mesh.elements[:,2],md.mesh.elements[:,0]]) + values=np.ones(np.shape(indi)) + + md.qmu.adjacency=m.sparse(indi,indj,values,md.mesh.numberofvertices,md.mesh.numberofvertices) + md.qmu.adjacency=np.logical_or(md.qmu.adjacency, md.qmu.adjacency.T).astype(float) #change to reshape(-1,1) if needed + + #now, build vwgt: + areas=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y) + + #get node connectivity + md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices) + + connectivity=md.mesh.vertexconnectivity[0][:,0:-1] + pos=np.where(connectivity) + connectivity[pos]=areas[connectivity[pos]-1]/3. + md.qmu.vertex_weight=np.sum(connectivity,1) + + return md Index: /issm/trunk-jpl/src/m/partition/partitioner.py =================================================================== --- /issm/trunk-jpl/src/m/partition/partitioner.py (revision 23095) +++ /issm/trunk-jpl/src/m/partition/partitioner.py (revision 23095) @@ -0,0 +1,124 @@ +import numpy as np +from pairoptions import * +import MatlabFuncs as m +from adjacency import * +#from Chaco import * +#from Scotch import * +#from MeshPartition import * +from project3d import * + +def partitioner(md,*varargin): + help =''' +PARTITIONER - partition mesh + + List of options to partitioner: + + package: 'chaco', 'metis' + npart: number of partitions. + weighting: 'on' or 'off': default off + section: 1 by defaults(1=bisection, 2=quadrisection, 3=octasection) + recomputeadjacency: 'on' by default (set to 'off' to compute existing one) + Output: md.qmu.partition recover the partition vector + + Usage: + md=partitioner(md,'package','chaco','npart',100,'weighting','on') + ''' + + #get options: + options=pairoptions(*varargin) + + #set defaults + options.addfielddefault('package','chaco') + options.addfielddefault('npart',10) + options.addfielddefault('weighting','on') + options.addfielddefault('section',1) + options.addfielddefault('recomputeadjacency','on') + + #get package: + package=options.getfieldvalue('package') + npart=options.getfieldvalue('npart') + recomputeadjacency=options.getfieldvalue('recomputeadjacency') + + if(np.ndim(md.mesh)==3): + #partitioning essentially happens in 2D. So partition in 2D, then + #extrude the partition vector vertically. + md3d=md #save for later + md.mesh.elements=md.mesh.elements2d + md.mesh.x=md.mesh.x2d + md.mesh.y=md.mesh.y2d + md.mesh.numberofvertices=md.mesh.numberofvertices2d + md.mesh.numberofelements=md.mesh.numberofelements2d + md.qmu.vertex_weight=[] + md.mesh.vertexconnectivity=[] + recomputeadjacency='on' + + #adjacency matrix if needed: + if m.strcmpi(recomputeadjacency,'on'): + md=adjacency(md) + else: + print 'skipping adjacency matrix computation as requested in the options' + + if m.strcmpi(package,'chaco'): + raise RuntimeError('Chaco is not currently supported for this function') + + # default method (from chaco.m) + #method=np.array([1,1,0,0,1,1,50,0,.001,7654321]).reshape(-1,1) + #method[0]=3 # global method (3=inertial (geometric)) + #method[2]=0 # vertex weights (0=off, 1=on) + + #specify bisection + #method[5]=options.getfieldvalue('section')# ndims (1=bisection, 2=quadrisection, 3=octasection) + + #are we using weights? + #if m.strcmpi(options.getfieldvalue('weighting'),'on'): + #weights=np.floor(md.qmu.vertex_weight/min(md.qmu.vertex_weight)) + #method[2]=1 + #else: + #weights=[] + + # partition into nparts + #if isinstance(md.mesh,mesh2d): + #part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y,np.zeros((md.mesh.numberofvertices,)),method,npart,[]).T+1 #index partitions from 1 up. like metis. + #else: + #part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y,md.mesh.z,method,npart,[]).T+1 #index partitions from 1 up. like metis. + + + elif m.strcmpi(package,'scotch'): + raise RuntimeError('Scotch is not currently supported for this function') + + #are we using weights? + #if m.strcmpi(options.getfieldvalue('weighting'),'on'): + #weights=np.floor(md.qmu.vertex_weight/min(md.qmu.vertex_weight)) + #else: + #weights=[] + + #maptab=Scotch(md.qmu.adjacency,[],weights,[],'cmplt',[npart]) + + #part=maptab[:,1]+1#index partitions from 1 up. like metis. + + elif m.strcmpi(package,'linear'): + + if (npart == md.mesh.numberofelements) or (md.qmu.numberofpartitions == md.mesh.numberofelements): + part=np.arange(1,1+md.mesh.numberofelements,1) + print 'Linear partitioner requesting partitions on elements' + else: + part=np.arange(1,1+md.mesh.numberofvertices,1) + + elif m.strcmpi(package,'metis'): + raise RuntimeError('Metis/MeshPartition is not currently supported for this function') + #[element_partitioning,part]=MeshPartition(md,md.qmu.numberofpartitions) + + else: + print help + raise RuntimeError('partitioner error message: could not find '+str(package)+' partitioner') + + #extrude if we are in 3D: + if np.ndim(md.mesh)==3: + md3d.qmu.vertex_weight=md.qmu.vertex_weight + md3d.qmu.adjacency=md.qmu.adjacency + md=md3d + part=project3d(md,'vector',part.T,'type','node') + + md.qmu.partition=part + + return md Index: /issm/trunk-jpl/src/m/qmu/dakota_in_data.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/dakota_in_data.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/dakota_in_data.py (revision 23095) @@ -0,0 +1,90 @@ +#move this stuff elsewhere +from helpers import * + +from dakota_in_write import * +from dakota_in_params import * +from MatlabFuncs import * + +def dakota_in_data(dmeth,variables,responses,dparams,filei,*args): + ''' + define the data to write the dakota .in and .m files. + + []=dakota_in_data(dmeth,variables,responses,dparams,filei,*args) + + where the required input is: + dmeth (dakota_method, method class object) + variables (structure array, variable class objects) + responses (structure array, response class objects) + dparams (structure array, method-independent parameters) + filei (character, name of .in and .m files) + + params may be empty, in which case defaults will be used. + + the optional args are passed directly through to the + QmuUpdateFunctions brancher to be used by the analysis + package. for example, this could be model information. + + this function defines the data to write the dakota .in and + .m files. it is necessary for multiple reasons. first, + it collects the parameters and applies some defaults that + are unique to the environment. second, some analysis package + variables and/or responses may be treated differently by + dakota. for example, an analysis package variable may be + defined as an array, so the QmuSetupDesign brancher will + create dakota variables for each element of the array. + finally it calls the functions to write the .in and .m files. + this function is independent of the particular analysis + package. + + this data would typically be generated by a matlab script + for a specific model, using the method, variable, and + response class objects. +''' + + ## 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 in range(np.size(fnames)): + if not isfield(params,fnames[i]): + print 'WARNING: dakota_in_data:unknown_param: No parameter '+str(fnames[i])+' in default parameter set.' + + exec('params.%s = vars(dparams)[fnames[i]]')%(fnames[i]) + + # use matlab even though we are running python + if params.direct and params.analysis_driver == '': + params.analysis_driver='matlab' + + if strcmpi(params.analysis_driver,'matlab') and params.analysis_components == '': + [pathstr,name,ext] = fileparts(filei) + params.analysis_components=fullfile(pathstr,name+'.py') + + # merge method parameters, though they shouldn't be in dparams + # dmeth=dmeth_params_merge(dmeth,dparams) + + ## variables + fnames=fieldnames(variables) + + # works like matlab arbitrary structs/classes, remembers order of input attributes + dvar = OrderedStruct() + dresp = OrderedStruct() + + for i in range(len(fnames)): + # currently all variable types can just be copied + exec('dvar.%s = vars(variables)[fnames[i]]')%(fnames[i]) + + ## responses + fnames=fieldnames(responses) + + for i in range(len(fnames)): + # currently all response types can just be copied + exec('dresp.%s = vars(responses)[fnames[i]]')%(fnames[i]) + + + ## write files + #Write in file + dakota_in_write(dmeth,dvar,dresp,params,filei,*args) Index: /issm/trunk-jpl/src/m/qmu/dakota_in_params.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/dakota_in_params.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/dakota_in_params.py (revision 23095) @@ -0,0 +1,184 @@ +from dakota_in_data import * + +#move this later: +from helpers import * + +def dakota_in_params(params): + ''' + populate a Dakota parameter structure. + + params=dakota_in_params(params) + + where the optional input is: + params (structure array, method-independent parameters) + + and the output is the same. + + this function takes a structure of method-independent dakota + parameters, which may be empty, and adds default parameters + for those parameters which do not exist. + + the field names of the structure are identical to the dakota + parameter names (and are in fact used to write them to the + files). logical values are used for parameters which have + no associated data and are determined only by their presence + or absence. + + note that the method-dependent parameters are contained in + the dakota_method class object. +''' + if params == None: + help(dakota_in_params) + return + + ## process the input parameters + if len(fieldnames(params)) == 0: + params=struct() + + ## strategy section + if not isfield(params,'graphics'): + params.graphics=False + + if not isfield(params,'tabular_graphics_data'): + params.tabular_graphics_data=False + + # could use unique file name rather than 'dakota_tabular.dat' + if not isfield(params,'tabular_graphics_file'): + params.tabular_graphics_file=False + + ## method section + # nearly all method parameters are in the dakota_method class + # or result from the response level lists + if not isfield(params,'compute'): + params.compute='probabilities' + + if not isfield(params,'distribution'): + params.distribution='cumulative' + + ## model section + + ## interface section + if not isfield(params,'system'): + params.system=False + + if not isfield(params,'fork'): + params.fork=False + + if not isfield(params,'direct'): + params.direct=False + + # interface parallelism controls + if not isfield(params,'asynchronous'): + params.asynchronous=True + + if not isfield(params,'evaluation_concurrency'): + params.evaluation_concurrency=False + + if not isfield(params,'analysis_concurrency'): + params.analysis_concurrency=False + + if not isfield(params,'evaluation_servers'): + params.evaluation_servers=False + + if not isfield(params,'evaluation_self_scheduling'): + params.evaluation_self_scheduling=False + + if not isfield(params,'evaluation_static_scheduling'): + params.evaluation_static_scheduling=True + + if not isfield(params,'evaluation_scheduling'): + params.evaluation_scheduling=False + + if not isfield(params,'processors_per_evaluation'): + params.processors_per_evaluation=False + + if not isfield(params,'analysis_servers'): + params.analysis_servers=False + + if not isfield(params,'analysis_self_scheduling'): + params.analysis_self_scheduling=False + + if not isfield(params,'analysis_static_scheduling'): + params.analysis_static_scheduling=False + + # algebraic mappings + if not isfield(params,'algebraic_mappings'): + params.algebraic_mappings=False + + # simulation interface controls + if not isfield(params,'analysis_driver'): + params.analysis_driver='' + + if not isfield(params,'analysis_components'): + params.analysis_components='' + + if not isfield(params,'input_filter'): + params.input_filter='' + + if not isfield(params,'output_filter'): + params.output_filter='' + + if not isfield(params,'failure_capture'): + params.failure_capture='abort' + + if not isfield(params,'deactivate'): + params.deactivate='evaluation_cache restart_file' + + # system call or fork interface + if not isfield(params,'parameters_file'): + params.parameters_file='params.in' + + if not isfield(params,'results_file'): + params.results_file='results.out' + + if not isfield(params,'verbatim'): + params.verbatim=False + + if not isfield(params,'aprepro'): + params.aprepro=False + + if not isfield(params,'file_tag'): + params.file_tag=True + + if not isfield(params,'file_save'): + params.file_save=True + + # direct function interface + if not isfield(params,'processors_per_analysis'): + params.processors_per_analysis=False + + ## responses section + if not isfield(params,'numerical_gradients'): + params.numerical_gradients=False + + if not isfield(params,'method_source'): + params.method_source='dakota' + + if not isfield(params,'interval_type'): + params.interval_type='forward' + + if not isfield(params,'fd_gradient_step_size'): + params.fd_gradient_step_size=0.001 + + if not isfield(params,'analytic_gradients'): + params.analytic_gradients=False + + # mixed_gradients not fully implemented + if not isfield(params,'mixed_gradients'): + params.mixed_gradients=False + + if not isfield(params,'id_analytic_gradients'): + params.id_analytic_gradients=False + + if not isfield(params,'id_numerical_gradients'): + params.id_numerical_gradients=False + + # hessians not fully implemented + if not isfield(params,'numerical_hessians'): + params.numerical_hessians=True + + if not isfield(params,'hessian_gradient_step_size'): + params.hessian_gradient_step_size=0.001 + + return params + Index: /issm/trunk-jpl/src/m/qmu/dakota_in_write.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/dakota_in_write.m (revision 23094) +++ /issm/trunk-jpl/src/m/qmu/dakota_in_write.m (revision 23095) @@ -156,5 +156,4 @@ % write response levels - if strcmp(dmeth.type,'nond') for i=1:length(dmeth.responses) Index: /issm/trunk-jpl/src/m/qmu/dakota_in_write.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/dakota_in_write.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/dakota_in_write.py (revision 23095) @@ -0,0 +1,333 @@ +#Move this somewhere else later +from helpers import * + +from dakota_method import * +from dakota_in_params import * +from IssmConfig import * +from dmeth_params_write import * +from vector_write import * +from qmu_classes import * + +import itertools + +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. +''' + + # process the input parameters + if len(fieldnames(method)) == 0: + method=str(input('Method? ')) + + if type(method) == str: + dmeth=dakota_method(method) + elif isinstance(method,dakota_method): + dmeth=method + else: + raise RuntimeError('Method '+str(method)+' is unrecognized class '+str(type(method))+'. (should be either "str" or "dakota_method")') + + if len(filei) == 0: + filei=str(input('Dakota input file to write? ')) + + [pathstr,name,ext] = fileparts(filei) + if len(ext) == 0: + # fileparts only considers '.in' to be the extension, not '.qmu.in' + ext='.qmu.in' + + filei2=fullfile(pathstr,name+ext) + + print 'Opening Dakota input file \''+filei2 + '\'' + try: + with open(filei2,'w+') as fidi: + + if len(fieldnames(params)) == 0: + params=struct() + + params=dakota_in_params(params) + + # write the strategy section + if float(IssmConfig('_DAKOTA_VERSION_')[0]) < 6: + strategy_write(fidi,params) + else: + environment_write(fidi,params) + + # write the method section + method_write(fidi,dmeth,dresp,params) + + # write the model section + model_write(fidi) + + # write the variables section + variables_write(fidi,dmeth,dvar) + + # write the interface section + interface_write(fidi,params) + + # write the responses section + responses_write(fidi,dmeth,dresp,params) + + except IOError: + print filei2 + ' could not be opened.' + + print 'End of file successfully written.' + + +## function to write the strategy section of the file + +def strategy_write(fidi,params): + + print 'Writing strategy section of Dakota input file.' + + fidi.write('strategy,\n') + fidi.write('\tsingle_method\n\n') + param_write(fidi,'\t ','graphics','','\n',params) + param_write(fidi,'\t ','tabular_graphics_data','','\n',params) + param_write(fidi,'\t ','tabular_graphics_file',' ','\n',params) + fidi.write('\n') + + +## function to write the environment section of the file + +def environment_write(fidi,params): + + print 'Writing environment section of Dakota input file.' + + fidi.write('environment,\n') + param_write(fidi,'\t ','graphics','','\n',params) + param_write(fidi,'\t ','tabular_graphics_data','','\n',params) + param_write(fidi,'\t ','tabular_graphics_file',' ','\n',params) + fidi.write('\n') + + +## function to write the method section of the file + +def method_write(fidi,dmeth,dresp,params): + + print 'Writing method section of Dakota input file.' + + fidi.write('method,\n') + fidi.write('\t'+str(dmeth.method)+'\n') + + dmeth_params_write(dmeth,fidi) + + # write response levels + + if strcmp(dmeth.type,'nond'): + for i in range(len(dmeth.responses)): + str_name = dmeth.responses[i] + resp = eval("{}.{}()".format(str_name,str_name)) + resp.dakota_rlev_write(fidi,dresp,params) + + fidi.write('\n') + + +## function to write the model section of the file + +def model_write(fidi): + + print 'Writing model section of Dakota input file.' + + fidi.write('model,\n') + fidi.write('\tsingle\n\n') + + +## function to write the variables section of the file + +def variables_write(fidi,dmeth,dvar): + + print 'Writing variables section of Dakota input file.' + + fidi.write('variables,\n') + + # variables vary by method + fd = fieldnames(dvar) + types = [] + var = [] + for i in range(len(fd)): + i_type = eval('dvar.{}[0].__class__.__name__'.format(fd[i])) + j = dmeth.variables.index(i_type) + str_name = dmeth.variables[j] + + # organize so that multiple instances of the same qmu class + # (2 different variable instances of "normal_uncertain" for example) + # are in the same dakota_write call regardless of individual size; + # but that each class has its own dakota_write call + if str_name not in types: + types.append(str_name) + var.append(eval('dvar.{}'.format(fd[i]))) + else: + t = types.index(str_name) + var[t].extend(eval('dvar.{}'.format(fd[i]))) + + for t in range(len(types)): + v = eval('{}.{}()'.format(types[t],types[t])) + v.dakota_write(fidi,var[t]) + + # linear constraints vary by method + fc = dmeth.lcspec + + for i in range(len(dmeth.lcspec)): + str_name = dmeth.lcspec[i] + var = eval('{}.{}()'.format(str_name,str_name)) + # check that str_name is correct against matlab version which has no argument there + var.dakota_write(fidi,eval('dvar.{}[i]'.format(j)),str_name) + + fidi.write('\n') + + +## function to write the interface section of the file + +def interface_write(fidi,params): + + print 'Writing interface section of Dakota input file.' + + fidi.write('interface,\n') + + if (not params.system) and (not params.fork) and (not params.direct): + params.fork=True + elif params.system+params.fork+params.direct > 1: + raise RuntimeError('Too many interfaces selected.') + + if params.system or params.fork: + param_write(fidi,'\t','asynchronous','','\n',params) + param_write(fidi,'\t ','evaluation_concurrency',' = ','\n',params) + param_write(fidi,'\t ','analysis_concurrency',' = ','\n',params) + param_write(fidi,'\t ','evaluation_servers',' = ','\n',params) + param_write(fidi,'\t ','evaluation_self_scheduling','','\n',params) + param_write(fidi,'\t ','evaluation_static_scheduling','','\n',params) + param_write(fidi,'\t ','analysis_servers',' = ','\n',params) + param_write(fidi,'\t ','analysis_self_scheduling','','\n',params) + param_write(fidi,'\t ','analysis_static_scheduling','','\n',params) + param_write(fidi,'\t','algebraic_mappings',' = ','\n',params) + param_write(fidi,'\t','system','','\n',params) + param_write(fidi,'\t','fork','','\n',params) + param_write(fidi,'\t ','analysis_driver',' = \'','\'\n',params) + if len(params.input_filter) != 0: + param_write(fidi,'\t ','input_filter',' = ','\n',params) + + if len(params.output_filter) != 0: + param_write(fidi,'\t ','output_filter',' = ','\n',params) + + param_write(fidi,'\t ','failure_capture',' ','\n',params) + param_write(fidi,'\t ','deactivate',' ','\n',params) + param_write(fidi,'\t ','parameters_file',' = ','\n',params) + param_write(fidi,'\t ','results_file',' = ','\n',params) + param_write(fidi,'\t ','verbatim', '','\n',params) + param_write(fidi,'\t ','aprepro', '','\n',params) + param_write(fidi,'\t ','file_tag', '','\n',params) + param_write(fidi,'\t ','file_save','','\n',params) + elif params.direct: + # Error: asynchronous capability not yet supported in direct interfaces. + # Update: it is now possible to run in parallel in direct interfaces. + param_write(fidi,'\t','algebraic_mappings',' = ','\n',params) + param_write(fidi,'\t','direct','','\n',params) + param_write(fidi,'\t ','analysis_driver',' = \'','\'\n',params) + if IssmConfig('_DAKOTA_VERSION_') < 6: + param_write(fidi,'\t ','evaluation_static_scheduling','','\n',params) + else: + param_write(fidi,'\t ','evaluation_scheduling',' ','\n',params) + param_write(fidi,'\t ','processors_per_evaluation',' = ','\n',params) + if len(params.analysis_components) != 0: + [pathstr,name,ext] = fileparts(params.analysis_components) + if ext != '': + ext='.py' + + params.analysis_components=fullfile(pathstr,name+ext) + param_write(fidi,'\t ','analysis_components',' = \'','\'\n',params) + + if len(params.input_filter) != 0: + param_write(fidi,'\t ','input_filter',' = ','\n',params) + + if len(params.output_filter) != 0: + param_write(fidi,'\t ','output_filter',' = ','\n',params) + + param_write(fidi,'\t ','failure_capture',' ','\n',params) + param_write(fidi,'\t ','deactivate',' ','\n',params) + param_write(fidi,'\t ','processors_per_analysis',' = ','\n',params) + + fidi.write('\n') + + +## function to write the responses section of the file + +def responses_write(fidi,dmeth,dresp,params): + + print 'Writing responses section of Dakota input file.' + + fidi.write('responses,\n') + #fidi.write('calibration_terms = 1 \n') + + # functions, gradients, and hessians vary by method + + rdesc=[] + + for i in range(len(dmeth.responses)): + resp = eval(dmeth.responses[i]) + rdesc = resp.dakota_write(fidi,dresp,rdesc) + + # write accumulated response descriptors for all response classes + + if len(rdesc) != 0: + fidi.write('\tresponse_descriptors =\n') + vector_write(fidi,'\t ',rdesc,6,76) + + ghspec_write(fidi,params,dmeth.ghspec) + + fidi.write('\n') + + +## function to write gradient and hessian specifications + +def ghspec_write(fidi,params,ghspec): + + # gradients + if 'grad' in ghspec: + if (not params.numerical_gradients) and (not params.analytic_gradients): + params.numerical_gradients=True + elif (params.numerical_gradients+params.analytic_gradients > 1): + raise RuntimeError('Too many gradients selected.') + + if params.numerical_gradients: + param_write(fidi,'\t','numerical_gradients','','\n',params) + param_write(fidi,'\t ','method_source',' ','\n',params) + param_write(fidi,'\t ','interval_type',' ','\n',params) + param_write(fidi,'\t ','fd_gradient_step_size',' = ','\n',params) + elif params.analytic_gradients: + param_write(fidi,'\t','analytic_gradients','','\n',params) + # elif params.mixed_gradients + else: + fidi.write('\tno_gradients\n') + + # hessians (no implemented methods use them yet) + if 'hess' in ghspec: + raise RuntimeError('Hessians needed by method but not provided.') + else: + fidi.write('\tno_hessians\n') + + + Index: /issm/trunk-jpl/src/m/qmu/dakota_out_parse.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/dakota_out_parse.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/dakota_out_parse.py (revision 23095) @@ -0,0 +1,959 @@ +import numpy as np + +from os.path import isfile,getsize +import re + +from MatlabFuncs import * +from prctile_issm import * +from normfit_issm import * +#move this later +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 + +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. +''' + if filei == None: + help(dakota_out_parse) + return + + if not isfile(filei) or getsize(filei) == 0: + filei=str(input('Input file? ')) + + #fidi=fopen(sprintf('%s',filei),'r') + #try: + with open(filei, 'r') as fidi: + ## check the first line for the Dakota tabular output file + method=[] + fline=fidi.readline() + if getsize(filei) == 0 or fline == '': + raise RuntimeError('File '+filei+' is empty.') + + dresp=[] # of struct() + scm =struct() + pcm =struct() + srcm=struct() + prcm=struct() + + if strncmpi(fline,'%eval_id',8): + method='unknown' + dresp=dak_tab_out(fidi,fline) + return [method,dresp,scm,pcm,srcm,prcm] + else: + fidi.seek(0,0) + + ## loop through the file to find the Dakota method name + fline=findline(fidi,'method',True) + if fline == None: + #do nothing + pass + else: + if fline[6] == ',': + fline = fidi.readline() + [ntokens,tokens]=fltokens(fline) + method=tokens[0].strip() + 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 + counter = 0 + fline='' + nfeval=nfeval_read(fidi,fline) + + ## process each results section based on content of the file + while counter < 10: + # 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) + if counter >= 10: + break + if fline == '' or fline.isspace(): + counter += 1 + else: + counter = 0 + + #print fline + # ipos=ftell(fidi) + fline = fidi.readline() + if fline == '' or fline.isspace(): + pass + elif strncmp(fline,'<<<<< Function evaluation summary',33): + nfeval=nfeval_read(fidi,fline) + elif strncmp(fline,'Statistics based on ',20): + nsamp=nsamp_read(fidi,fline) + elif strncmp(fline,'Moments for each response function',34): + dresp=moments_read(fidi,dresp,fline) + elif strncmp(fline,'Moment-based statistics for each response function',50): + dresp=mbstats_read(fidi,dresp,fline) + elif strncmp(fline,'95% confidence intervals for each response function',51): + dresp=cis_read(fidi,dresp,fline) + elif strncmp(fline,'Probabilities for each response function',40) or strncmp(fline,'Level mappings for each response function',41): + dresp=cdfs_read(fidi,dresp,fline) + elif strncmp(fline,'Probability Density Function (PDF) histograms for each response function',72): + dresp=pdfs_read(fidi,dresp,fline) + elif strncmp(fline,'Simple Correlation Matrix',25): + scm=corrmat_read(fidi,'Simple Correlation Matrix',fline) + elif strncmp(fline,'Partial Correlation Matrix',26): + pcm=corrmat_read(fidi,'Partial Correlation Matrix',fline) + elif strncmp(fline,'Simple Rank Correlation Matrix',30): + srcm=corrmat_read(fidi,'Simple Rank Correlatio:n Matrix',fline) + elif strncmp(fline,'Partial Rank Correlation Matrix',31): + prcm=corrmat_read(fidi,'Partial Rank Correlation Matrix',fline) + elif strncmp(fline,'MV Statistics for ',18): + dresp=mvstats_read(fidi,dresp,fline) + elif strncmp(fline,'<<<<< Best ',11): + dresp=best_read(fidi,dresp,fline) + elif strncmp(fline,'The following lists volumetric uniformity measures',50): + dresp=vum_read(fidi,dresp,fline) + elif strncmp(fline,'<<<<< Iterator ',15) and (len(fline) > 26) and (' completed.' in fline[15:]): + method=itcomp_read(fidi,fline) + elif strncmp(fline,'-----',5): + pass + else: + 'Unexpected line: '+str(fline) + + # fidi.seek(ipos,0) + + ## loop through the file to verify the end + + # fline=findline(fidi,'<<<<< Single Method Strategy completed') + # if not ischar(fline) + # return + # + 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.') + + return [method,dresp,scm,pcm,srcm,prcm] + # ]]] + +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) + [ntokens,tokens]=fltokens(fline) + + # New file DAKOTA versions>6 + if strncmpi(fline,'%eval_id interface',18): + offset=2 + else: #DAKOTA versions<6 + offset=1 + + 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 + nrow=0 + while True: + + fline=fidi.readline() + + if fline == '' or fline.isspace(): + break + + if nrow > 0: + data = np.concatenate((data,[np.zeros(ntokens-offset)])) + + [ntokens,tokens]=fltokens(fline) + + # add row values to matrix (skipping eval_id) + + for i in range(ntokens-offset): + data[nrow,i] = tokens[i+offset] + + 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 + if (np.size(data,0) > 1): + #dmean =mean (data) + #dstddev=std (data,0) + [dmean,dstddev,dmeanci,dstddevci]=normfit_issm(data,0.05) + else: + dmean =np.zeros((1,np.size(data,1))) + dstddev =np.zeros((1,np.size(data,1))) + dmeanci =np.zeros((2,np.size(data,1))) + dstddevci=np.zeros((2,np.size(data,1))) + for i in range(np.size(data,1)): + [dmean[0,i],dstddev[0,i],dmeanci[:,i],dstddevci[:,i]]=normfit_issm(data[:,i],0.05) + + dmin =data.min(0) + dquart1=prctile_issm(data,25,0) + dmedian=np.median(data,0) + dquart3=prctile_issm(data,75,0) + dmax =data.max(0) + dmin95 =prctile_issm(data,5,0) + dmax95 =prctile_issm(data,95,0) + + dcorrel=np.corrcoef(data.T) + + # 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].sample =data[:,i] + dresp[i].mean =dmean[i] + dresp[i].stddev =dstddev[i] + dresp[i].meanci =dmeanci[:,i] + dresp[i].stddevci =dstddevci[:,i] + dresp[i].min =dmin[i] + dresp[i].quart1 =dquart1[i] + dresp[i].median =dmedian[i] + dresp[i].quart3 =dquart3[i] + dresp[i].max =dmax[i] + dresp[i].dmin95 =dmin95[i] + dresp[i].dmax95 =dmax95[i] + + # draw box plot + + # figure + # subplot(2,1,1) + # plot_boxplot(dresp) + + # draw normal probability plot + + # subplot(2,1,2) + # plot_normplot(dresp) + + return dresp + # ]]] + +def nfeval_read(fidi,fline): # [[[ +## function to find and read the number of function evaluations + + if fline == None or fline == '' or fline.isspace(): + fline=findline(fidi,'<<<<< Function evaluation summary') + nfeval = 0 + return + + [ntokens,tokens]=fltokens(fline) + nfeval=tokens[4] + print ' Dakota function evaluations='+str(int(nfeval))+'.' + + return nfeval + # ]]] + +def nsamp_read(fidi,fline): # [[[ +## function to find and read the number of samples + + if fline == None or fline == '' or fline.isspace(): + fline=findline(fidi,'Statistics based on ') + return + + [ntokens,tokens]=fltokens(fline) + nsamp=tokens[3] + print ' Dakota samples='+str(int(nsamp))+'.' + + return nsamp + # ]]] + +def moments_read(fidi,dresp,fline): # [[[ +## function to find and read the moments + + if fline == None or fline == '' or fline.isspace(): + fline=findline(fidi,'Moments for each response function') + return + + print 'Reading moments for response functions:' + + while True: + fline=fidi.readline() + if fline == '' or fline.isspace(): + break + + [ntokens,tokens]=fltokens(fline) + + # add new response function and moments + + dresp.append(struct()) + dresp[-1].descriptor=tokens[ 0] + print ' '+str(dresp[-1].descriptor) + dresp[-1].mean =tokens[ 3] + dresp[-1].stddev =tokens[ 6] + dresp[-1].coefvar =tokens[12] + + print ' Number of Dakota response functions='+str(len(dresp))+'.' + + return dresp + # ]]] + +def mbstats_read(fidi,dresp,fline): # [[[ +## function to find and read the moment-based statistics + + if fline == None or fline == '' or fline.isspace(): + fline=findline(fidi,'Moment-based statistics for each response function') + return + + print 'Reading moment-based statistics for response functions:' + + # skip column headings of moment-based statistics + + fline=fidi.readline() + + while True: + fline=fidi.readline() + if fline == '' or fline.isspace(): + break + + [ntokens,tokens]=fltokens(fline) + + # add new response function and moment-based statistics + + dresp.append(struct()) + dresp[-1].descriptor=tokens[ 0] + print ' '+str(dresp[-1].descriptor) + dresp[-1].mean =tokens[ 1] + dresp[-1].stddev =tokens[ 2] + dresp[-1].skewness =tokens[ 3] + dresp[-1].kurtosis =tokens[ 4] + + print ' Number of Dakota response functions='+str(len(dresp))+'.' + + return dresp + # ]]] + +def cis_read(fidi,dresp,fline): # [[[ +## function to find and read the confidence intervals + + if fline == None or fline == '' or fline.isspace(): + fline=findline(fidi,'95% confidence intervals for each response function') + return + + print 'Reading 95% confidence intervals for response functions:' + + while True: + fline=fidi.readline() + if fline == '' or fline.isspace(): + break + + [ntokens,tokens]=fltokens(fline) + + # check for column headings in Dakota 5.2 + + if (ntokens == 4): + fline=fidi.readline() + if fline == '' or fline.isspace(): + break + + [ntokens,tokens]=fltokens(fline) + + # find response function associated with confidence intervals + + idresp=-1 + for i in range(len(dresp)): + if strcmpi(tokens[0],dresp[i].descriptor): + idresp=i + break + + if idresp < 0: + idresp=len(dresp) + dresp.append(struct()) + dresp[idresp].descriptor=tokens[0] + print ' '+str(dresp[idresp].descriptor) + + # add confidence intervals to response functions + dresp[i].meanci = np.array([[np.nan],[np.nan]]) + dresp[i].stddevci = np.array([[np.nan],[np.nan]]) + + if (ntokens == 14): + dresp[i].meanci [0,0]=tokens[ 4] + dresp[i].meanci [1,0]=tokens[ 5] + dresp[i].stddevci[0,0]=tokens[11] + dresp[i].stddevci[1,0]=tokens[12] + else: + dresp[i].meanci [0,0]=tokens[ 1] + dresp[i].meanci [1,0]=tokens[ 2] + dresp[i].stddevci[0,0]=tokens[ 3] + dresp[i].stddevci[1,0]=tokens[ 4] + + print ' Number of Dakota response functions='+str(len(dresp))+'.' + + return dresp + # ]]] + +def cdfs_read(fidi,dresp,fline): # [[[ +## function to find and read the cdf's + + if fline == None or fline == '' or fline.isspace(): + fline=findline(fidi,'Probabilities for each response function') + if fline == None: + fline=findline(fidi,'Level mappings for each response function') + if fline == None: + return + + print 'Reading CDF''s for response functions:' + + while fline == '' or fline.isspace(): + fline=fidi.readline() + if fline == '' or fline.isspace(): + break + + # process header line of cdf + + while (fline != '' and not fline.isspace()): + [ntokens,tokens]=fltokens(fline) + + # find response function associated with cdf + # idresp is an index, so it can be 0, default to -1 + idresp=-1 + for i in range(len(dresp)): + if strcmpi(tokens[ 5],dresp[i].descriptor): + idresp=i + break + + + if idresp < 0: + idresp=len(dresp) + dresp.append(struct()) + dresp[idresp].descriptor=tokens[ 5] + print ' '+str(dresp(idresp).descriptor) + + + # skip column headings of cdf + + fline=fidi.readline() + fline=fidi.readline() + + # read and add cdf table to response function + + fline=fidi.readline() + icdf=0 + while (fline != '' and not fline.isspace()) and not strncmpi(fline,'Cumulative Distribution Function',32): + [ntokens,tokens]=fltokens(fline) + icdf=icdf+1 + dresp[idresp].cdf = np.zeros((icdf,4)) + dresp[idresp].cdf[icdf-1,0:4]=np.nan + # in later versions of Dakota, uncalculated columns are now blank + itoken=0 + for i in range(len(fline)/19): + if not isempty(fline[(i-1)*19:i*19]): + itoken=itoken+1 + dresp[idresp].cdf[icdf-1,i]=tokens[itoken] + + fline=fidi.readline() + + print ' Number of Dakota response functions='+str(len(dresp))+'.' + + return dresp + # ]]] + +def pdfs_read(fidi,dresp,fline): # [[[ +## function to find and read the pdf's + + if fline == None or fline == '' or fline.isspace(): + fline=findline(fidi,'Probability Density Function (PDF) histograms for each response function') + return + + print 'Reading PDF''s for response functions:' + + while (fline != '' and not fline.isspace()): + fline=fidi.readline() + if fline == '' or fline.isspace(): + break + + + # process header line of pdf + + while (fline != '' and not fline.isspace()): + [ntokens,tokens]=fltokens(fline) + + # find response function associated with pdf + # idresp is an index, so it can be 0, default to -1 + idresp=-1 + for i in range(len(dresp)): + if strcmpi(tokens[ 2],dresp[i].descriptor): + idresp=i + break + + + if idresp < 0: + idresp=len(dresp) + dresp.append(struct) + dresp[idresp].descriptor=tokens[ 2] + print ' '+str(dresp[idresp].descriptor) + + + # skip column headings of pdf + + fline=fidi.readline() + fline=fidi.readline() + + # read and add pdf table to response function + + fline=fidi.readline() + ipdf=0 + while (fline != '' and not fline.isspace()) and not strncmpi(fline,'PDF for', 7): + [ntokens,tokens]=fltokens(fline) + ipdf=ipdf+1 + dresp[idresp].pdf = np.zeros((ipdf,4)) + dresp[idresp].pdf[ipdf-1,0:3]=np.nan + for i in range(3): + dresp[idresp].pdf[ipdf-1,i]=tokens[i] + + fline=fidi.readline() + + print ' Number of Dakota response functions='+str(len(dresp))+'.' + + return dresp + # ]]] + +def corrmat_read(fidi,cmstr,fline): # [[[ +## function to find and read a correlation matrix + + if fline == None or fline == '' or fline.isspace(): + fline=findline(fidi,cmstr) + if fline == '' or fline.isspace(): + cmat=struct() + return + + print 'Reading ' +fline+ '.' + + cmat.title=fline + + while (fline != '' and not fline.isspace()): + fline=fidi.readline() + if fline == '' or fline.isspace(): + break + + # process column headings of matrix + + [ntokens,tokens]=fltokens(fline) + cmat.column= np.empty((1,ntokens)) + cmat.column.fill(0.0) + cmat.row = np.empty((1,1)) + cmat.row.fill(0.0) + cmat.matrix=np.zeros((1,ntokens)) + + for i in range(ntokens): + cmat.column[1,i]=str(tokens[i]) + + # process rows of matrix, reading until blank line + + nrow=0 + while True: + fline=fidi.readline() + if fline == '' or fline.isspace(): + break + + [ntokens,tokens]=fltokens(fline) + + # add row heading to matrix + + nrow=nrow+1 + cmat.row[nrow-1,0]=str(tokens[0]) + + # add row values to matrix + + for i in range(1,ntokens): + cmat.matrix[nrow-1,i-1]=tokens[i] + + return cmat + # ]]] + +def mvstats_read(fidi,dresp,fline): # [[[ +## function to find and read the MV statistics + + if fline == None or fline == '' or fline.isspace(): + fline=findline(fidi,'MV Statistics for ') + if fline == None: + return + + print 'Reading MV statistics for response functions:' + + ndresp=0 + + while (fline != '' and not fline.isspace()) and strncmpi(fline,'MV Statistics for ',18): + + # add new response function and moments + + [ntokens,tokens]=fltokens(fline) + dresp.append(struct()) + dresp[-1].descriptor=tokens[3] + print ' '+str(dresp[-1].descriptor) + fline=fidi.readline() + [ntokens,tokens]=fltokens(fline) + dresp[-1].mean =tokens[4] + fline=fidi.readline() + [ntokens,tokens]=fltokens(fline) + dresp[-1].stddev =tokens[6] + + # read and add importance factors to response function + + idvar=0 + fline=fidi.readline() + if fline == '' or fline.isspace(): + break + + # shape: [[0],[0],[0]...] + dresp[-1].var = [] + dresp[-1].impfac = [] + dresp[-1].sens = [] + + while (fline != '' and not fline.isspace()) and strncmpi(fline,' Importance Factor for variable ',33): + [ntokens,tokens]=fltokens(fline) + idvar=idvar+1 + dresp[-1].var.append(str(tokens[4])) + dresp[-1].impfac.append(tokens[6]) + if (ntokens >= 10): + dresp[-1].sens.append(tokens[9]) + else: + dresp[-1].sens.append(np.nan) + + + fline=fidi.readline() + + # if importance factors missing, skip to cdf + + if not idvar: + print ' Importance Factors not available.' + dresp[-1].var =[] + dresp[-1].impfac=[] + dresp[-1].sens =[] + while type(fline) == str and (fline != '' and not fline.isspace()) and not strncmpi(fline,'Cumulative Distribution Function',32) and not strncmpi(fline,'MV Statistics for ',18) and not strncmp(fline,'-',1): + fline=fidi.readline() + + # process header line of cdf + + icdf=0 + + # If there is a warning it MAY involve a lot of spaces; skip over them + if fline == '' or fline.isspace(): + fline=fidi.readline() + # Usually: "Warning: negligible standard deviation renders CDF results suspect." + if strncmpi(fline,'Warn',4): + fline = fidi.readline() + if fline == '' or fline.isspace(): + fline = fidi.readline() + + while (fline != '' and not fline.isspace()) and strncmpi(fline,'Cumulative Distribution Function',32): + [ntokens,tokens]=fltokens(fline) + + # find response function associated with cdf + # idresp is an index, so it can be 0, default to -1 + idresp=-1 + for i in range(len(dresp)): + if strcmpi(tokens[ 5],dresp[i].descriptor): + idresp=i + break + + if idresp < 0: + idresp=len(dresp) + dresp.append(struct()) + dresp[idresp].descriptor=tokens[ 5] + print ' '+str(dresp[idresp].descriptor) + + # skip column headings of cdf + fline=fidi.readline() + fline=fidi.readline() + + # read and add cdf table to response function + + fline=fidi.readline() + while (fline != '' and not fline.isspace()) and not strncmpi(fline,'MV Statistics for ',18) and not strncmp (fline,'-',1): + [ntokens,tokens]=fltokens(fline) + icdf=icdf+1 + dresp[idresp].cdf = np.zeros((icdf,4)) + dresp[idresp].cdf[icdf-1,0]=tokens[0] + dresp[idresp].cdf[icdf-1,1]=tokens[1] + if (ntokens == 4): + dresp[idresp].cdf[icdf-1,2]=tokens[2] + dresp[idresp].cdf[icdf-1,3]=tokens[3] + else: + dresp[idresp].cdf[icdf-1,2]=np.nan + dresp[idresp].cdf[icdf-1,3]=np.nan + + fline=fidi.readline() + + # if cdf missing, skip to end of response function + + if not icdf: + 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))+'.' + + return dresp + # ]]] + +def best_read(fidi,dresp,fline): # [[[ +## function to find and read the best evaluation + + if fline == None or fline == '' or fline.isspace(): + fline=findline(fidi,'<<<<< Best ') + if fline == None: + return + + if isempty(dresp): + dresp.append(struct()) + dresp[-1].best=struct() + + print 'Reading values for best function evaluation:' + + while (fline != '' and not fline.isspace()) and strncmpi(fline,'<<<<< Best ',11): + [ntokens,tokens]=fltokens(fline) + + # read and add best parameter(s) + + if strncmpi(str(tokens[2]),'parameter', 9): + print ' '+fline + + fline=fidi.readline() + dresp.best.param =[] + dresp.best.descriptor='' + + while (fline != '' and not fline.isspace()) and not strncmpi(fline,'<<<<< Best ',11): + [ntokens,tokens]=fltokens(fline) + dresp.best.param.append([0]) + dresp.best.param[-1]= tokens[0] + dresp.best.descriptor= str(tokens[1]) + fline=fidi.readline() + + # read and add best objective function(s) + + elif strncmpi(str(tokens[2]),'objective', 9) and strncmpi(str(tokens[3]),'function' , 8): + print ' '+fline + + fline=fidi.readline() + dresp.best.of=[] + + while (fline != '' and not fline.isspace()) and not strncmpi(fline,'<<<<< Best ',11): + [ntokens,tokens]=fltokens(fline) + dresp.best.of.append(0) + dresp.best.of[-1]=tokens[0] + fline=fidi.readline() + + # read and add best residual norms + + elif strncmpi(str(tokens[2]),'residual', 8) and strncmpi(str(tokens[3]),'norm' , 4): + print ' '+fline + dresp.best.norm = tokens[ 5] + dresp.best.hnormsq= tokens[10] + + fline=fidi.readline() + + while (fline != '' and not fline.isspace()) and not strncmpi(fline,'<<<<< Best ',11): + fline=fidi.readline() + + # read and add best residual term(s) + + elif strncmpi(str(tokens[2]),'residual', 8) and strncmpi(str(tokens[3]),'term' , 4): + print ' '+fline + + fline=fidi.readline() + dresp.best.res=[] + + while (fline != '' and not fline.isspace()) and not strncmpi(fline,'<<<<< Best ',11): + [ntokens,tokens]=fltokens(fline) + dresp.best.res.append(0) + dresp.best.res[-1]= tokens[0] + fline=fidi.readline() + + # read and add best constraint value(s) + + elif strncmpi(str(tokens[2]),'constraint',10) and strncmpi(str(tokens[3]),'value' , 5): + print ' '+fline + + fline=fidi.readline() + dresp.best.nc=[] + + while (fline != '' and not fline.isspace()) and not strncmpi(fline,'<<<<< Best ',11): + [ntokens,tokens]=fltokens(fline) + dresp.best.nc.append(0) + dresp.best.nc[-1]= tokens[0] + fline=fidi.readline() + + # read and add best data captured + + elif strncmpi(str(tokens[2]),'data' , 4) and strncmpi(str(tokens[3]),'captured', 8): + print ' '+fline + dresp.best.eval= tokens[7] + + fline=fidi.readline() + + while (fline != '' and not fline.isspace()) and not strncmpi(fline,'<<<<< Best ',11): + fline=fidi.readline() + + # read until next best or blank or end + else: + print ' '+fline+' (ignored)' + + fline=fidi.readline() + + while (fline != '' and not fline.isspace()) and not strncmpi(fline,'<<<<< Best ',11): + fline=fidi.readline() + + return dresp + # ]]] + +def vum_read(fidi,dresp,fline): # [[[ +## function to find and read the volumetric uniformity measures + + if fline == None or fline == '' or fline.isspace(): + fline=findline(fidi,'The following lists volumetric uniformity measures') + if fline == None: + return + + if isempty(dresp): + dresp.append(struct()) + dresp[-1].vum=[] + + print 'Reading measures for volumetric uniformity.' + + fline=fidi.readline() + fline=fidi.readline() + + while (fline != '' and not fline.isspace()): + [ntokens,tokens]=fltokens(fline) + check = tokens[0].lower() + if check == 'chi': + dresp.vum.chi=tokens[3] + elif check == 'd': + dresp.vum.d =tokens[3] + elif check == 'h': + dresp.vum.h =tokens[3] + elif check == 'tau': + dresp.vum.tau=tokens[3] + + fline=fidi.readline() + + return dresp + # ]]] + +def itcomp_read(fidi,fline): # [[[ +## function to find and read the iterator completion + + if fline == None or fline == '' or fline.isspace(): + while True: + fline=findline(fidi,'<<<<< Iterator ') + if fline == None: + return + + if (len(fline) > 26) and not (' completed.' in fline[15:]): + break + + [ntokens,tokens]=fltokens(fline) + method=tokens[2] + print 'Dakota iterator \''+str(method)+'\' completed.' + + return method + # ]]] + +def findline(fidi,string,goto_line=False): # [[[ +## function to find a file line starting with a specified string +## by default, return to previous position, before search +## if final argument is True, return such that fidi.readline() will read the line +## immediately after the searched line (or the first line if the search failed) + + ipos=fidi.tell() + npos = 0 + for fline in fidi: + npos += len(fline) + if (strncmpi(fline,string,len(string))): + if goto_line: + fidi.seek(npos,0) + else: + fidi.seek(ipos,0) + return fline + + # issue warning and reset file position + print 'Warning: findline:str_not_found: String '+str(string)+' not found in file.' + fidi.seek(ipos,0) + return None + # ]]] + +def fltokens(fline): # [[[ + ## function to parse a file line into tokens + if fline == None: + ntokens=-1 + tokens=[] + return [None,None] + + if fline == '' or fline.isspace(): + ntokens=0 + tokens=[] + return [None,None] + + # split wherever ' ' (space) or ':' occur + strings = re.split(':| ',fline) + # remove blank strings + strings = filter(lambda a: (a != '' and not a.isspace()),strings) + + ntokens=0 + tokens = ['' for i in range(len(strings))] + + # try to format substrings to float where possible and count tokens and ignore invalid values + for i in range(len(strings)): + if isempty(strings[i]): + continue + + # if the string is a number, make it a float, otherwise leave it alone + try: + tokens[ntokens] = float(strings[i]) + except ValueError: + tokens[ntokens] = strings[i] + + ntokens=ntokens+1 + + return [ntokens,tokens] + # ]]] + + Index: /issm/trunk-jpl/src/m/qmu/expandresponses.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/expandresponses.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/expandresponses.py (revision 23095) @@ -0,0 +1,19 @@ +from QmuSetupResponses import * +from helpers import * + +def expandresponses(md,responses): + #EXPANDRESPONSES - expand responses + + fnames=fieldnames(responses) + + # maintain order attributes were added + dresp = OrderedStruct() + + for k in fnames: + v = eval('responses.{}'.format(k)) + exec 'dresp.{} = type(v)()'.format(k) + for j in range(len(v)): + #call setupdesign + exec 'dresp.{}=QmuSetupResponses(md,dresp.{},v[j])'.format(k,k) + + return dresp Index: /issm/trunk-jpl/src/m/qmu/expandvariables.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/expandvariables.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/expandvariables.py (revision 23095) @@ -0,0 +1,28 @@ +from QmuSetupVariables import * +from helpers import * + +from qmu_classes import * + +def expandvariables(md,variables): + + fnames=fieldnames(variables) + + # maintain order attributes were added + dvar = OrderedStruct() + + for k in fnames: + v = eval('variables.{}'.format(k)) + + # for linear constraints, just copy + if isinstance(v,linear_inequality_constraint) or isinstance(v,linear_equality_constraint): + exec 'dvar.{} = v'.format(k) + + # for variables, call the setup function + else: + exec 'dvar.{} = type(v)()'.format(k) + for j in range(len(v)): + #call setupdesign + exec 'dvar.{}=QmuSetupVariables(md,dvar.{},v[j])'.format(k,k) + + + return dvar Index: /issm/trunk-jpl/src/m/qmu/helpers.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/helpers.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/helpers.py (revision 23095) @@ -0,0 +1,318 @@ +import numpy as np +from collections import OrderedDict + +class struct(object): + '''An empty struct that can be assigned arbitrary attributes''' + pass + +class Lstruct(list): + ''' +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'] + +Note that 'x' returns the array and x.__dict__ will only return + attributes other than the array + +List-based and struct-based behaviors work normally, however they are referenced + as if the other does not exist; len(x) corresponds only to + the list component of x, len(x.a) corresponds to x.a, x.__dict__ + corresponds only to the non-x-list attributes + +Example uses: + +x = Lstruct(1,2,3,4) -> [1,2,3,4] + x.a = 'hello' + len(x) -> 4 + x.append(5) + len(x) -> 5 + x[2] -> 3 + x.a -> 'hello' + print x -> [1,2,3,4,5] + x.__dict__ -> {'a': 'hello'} + x.b = [6,7,8,9] + x.b[-1] -> 9 + len(x.b) -> 4 + +Other valid constructors: + x = Lstruct(1,2,3,a='hello') -> x.a -> 'hello', x -> [1,2,3] + x = Lstruct(1,2,3)(a='hello') + x = Lstruct([1,2,3],x='hello') + x = Lstruct((1,2,3),a='hello') + +Credit: https://github.com/Vectorized/Python-Attribute-List +''' + + def __new__(self, *args, **kwargs): + return super(Lstruct, self).__new__(self, args, kwargs) + + def __init__(self, *args, **kwargs): + if len(args) == 1 and hasattr(args[0], '__iter__'): + list.__init__(self, args[0]) + else: + list.__init__(self, args) + self.__dict__.update(kwargs) + + def __call__(self, **kwargs): + self.__dict__.update(kwargs) + return self + +class OrderedStruct(object): + ''' +A form of dictionary-like structure that maintains the + ordering in which its fields/attributes and their + corresponding values were added. + +OrderedDict is a similar device, however this class + can be used as an "ordered struct/class" giving + it much more flexibility in practice. It is + also easier to work with fixed valued keys in-code. + +Eg: +OrderedDict: # a bit clumsy to use and look at + x['y'] = 5 + +OrderedStruct: # nicer to look at, and works the same way + x.y = 5 + OR + x['y'] = 5 # supports OrderedDict-style usage + +Supports: len(x), str(x), for-loop iteration. +Has methods: x.keys(), x.values(), x.items(), x.iterkeys() + +Usage: + x = OrderedStruct() + x.y = 5 + x.z = 6 + OR + x = OrderedStruct('y',5,'z',6) + + # note below that the output fields as iterables are always + # in the same order as the inputs + + x.keys() -> ['y','z'] + x.values() -> [5,6] + x.items() -> [('y',6),('z',6)] + x.__dict__ -> [('y',6),('z',6)] + vars(x) -> [('y',6),('z',6)] + + x.y -> 5 + x['y'] -> 5 + x.z -> 6 + x['z'] -> 6 + + for i in x: # same as x.items() + print i + -> + ('x',5) + ('y',6) + + Note: to access internal fields use dir(x) + (input fields will be included, but + are not technically internals) + ''' + + def __init__(self,*args): + '''Provided either nothing or a series of strings, construct a structure that will, +when accessed as a list, return its fields in the same order in which they were provided''' + + # keys and values + self._k = [] + self._v = [] + + if len(args) == 0: + return + + if len(args) % 2 != 0: + raise RuntimeError('OrderedStruct input error: OrderedStruct(*args) call must have an even number of inputs, in key/value pairs') + + for a,b in zip(args[0::2],args[1::2]): + exec('self.%s = b')%(a) + return + + def __repr__(self): + s = 'OrderedStruct:\n\t' + for a,b in zip(self._k,self._v): + s += str(a) + ' : ' + str(b) + '\n\t' + return s + + def __len__(self): + return len(self._k) + + def __getattr__(self, attr): + # called when __getattribute__ fails + try: + # check if in keys, then access + _k = object.__getattribute__(self, '_k') + _v = object.__getattribute__(self, '_v') + pos = _k.index(attr) + return _v[pos] + except ValueError: + # not in keys, not a valid attribute, raise error + raise AttributeError('Attribute "'+str(attr)+'" does not exist.') + + def __getattribute__(self, attr): + # re-route calls to vars(x) and x.__dict__ + if attr == '__dict__': + return OrderedDict(self.items()) + else: + return object.__getattribute__(self, attr) + + def __getitem__(self, key): + return self._v[self._k.index(key)] + + def __setattr__(self, name, value): + #super(OrderedStruct, self).__setattr__(name,value) + if name in ['_k','_v']: + object.__setattr__(self, name, value) + elif name not in self._k: + self._k.append(name) + self._v.append(value) + else: + self._v[self._k.index(name)] = value + + def __delattr__(self, key): + if name not in self._k: + raise AttributeError('Attribute "'+str(attr)+'" does not exist or is an internal field and therefore cannot be deleted safely.') + self.pop(key) + + def __iter__(self): + for a,b in zip(self._k,self._v): + yield(a,b) + + def iterkeys(self): + for k in self._k: + yield k + + def pop(self,key): + i = self._k.index(key) + k = self._k.pop(i) + v = self._v.pop(i) + #exec('del self.%s')%(key) + return (k,v) + + def keys(self): + return self._k + def values(self): + return self._v + def items(self): + return zip(self._k,self._v) + +def isempty(x): + '''returns true if object is +\-infinity, NaN, None, '', has length 0, or is an + array/matrix composed only of such components (includes mixtures of "empty" types)''' + + if type(x) in [list,np.ndarray,tuple]: + if np.size(x) == 0: + return True + + # if anything in that array/matrix is not empty, the whole thing is not empty + try: + x = np.concatenate(x) + except (ValueError): + pass + for i in x: + if not isempty(i): + return False + # the array isn't empty but is full of "empty" type objects, so return True + return True + + if x == None: + return True + if type(x) == str and x.lower() in ['','nan','none','inf','infinity','-inf','-infinity']: + return True + + # type may not be understood by numpy, in which case it definitely is NOT NaN or infinity + try: + if np.isnan(x) or np.isinf(x): + return True + except (TypeError): + pass + + # if all of that fails, then if is not empty + return False + +def fieldnames(x,ignore_internals = True): + '''returns a list of fields of x + ignore_internals ignores all fieldnames starting with '_' and is True by default''' + result = vars(x).keys() + + if ignore_internals: + result = filter(lambda i: i[0] != '_',result) + + return result + +def isfield(x,y,ignore_internals=True): + '''is y is a field of x? + ignore_internals ignores all fieldnames starting with '_' and is True by default''' + return str(y) in fieldnames(x,ignore_internals) + +def fileparts(x): + ''' + given: "path/path/.../file_name.ext" + returns: [path,file_name,ext] (list of strings)''' + try: + a = x[:x.rindex('/')] #path + b = x[x.rindex('/')+1:] #full filename + except ValueError: #no path provided + a = '' + b = x + try: + c,d = b.split('.') #file name, extension + except ValueError: #no extension provided + return [a,b,''] + return [a,c,'.'+d] + +def fullfile(*args): + ''' + use: + + fullfile(path, path, ... , file_name+ext) + returns: "path/path/.../file_name.ext" + + with all arguments as strings with no "/"s + + regarding extensions and the '.': + as final arguments ('file.doc') or ('file'+'.doc') will work; + ('final','.doc'), and the like, will not (you'd get 'final/.doc') +''' + result = str(args[0]) + for i in range(len(args[1:])): + # if last argument wasn't empty, add a '/' between it and the next argument + if len(args[i]) != 0: + result += '/'+str(args[i+1]) + else: + result += str(args[i+1]) + return result + +def findline(fidi,s): + ''' + returns full first line containing s (as a string), or None + Note: will include any newlines or tabs that occur in that line, + use str(findline(f,s)).strip() to remove these, str() in case result is None''' + for line in fidi: + if s in line: + return line + return None + +def empty_nd_list(shape,filler=0.,as_numpy_ndarray=False): + ''' +returns a python list of the size/shape given (shape must be int or tuple) + the list will be filled with the optional second argument + + filler is 0.0 by default + + as_numpy_ndarray will return the result as a numpy.ndarray and is False by default + + Note: the filler must be either None/np.nan/float('NaN'), float/double, or int + other numpy and float values such as +/- np.inf will also work + +use: + empty_nd_list((5,5), 0.0) # returns a 5x5 matrix of 0.0's + empty_nd_list(5, None) # returns a 5 long array of NaN +''' + result = np.empty(shape) + result.fill(filler) + if not as_numpy_ndarray: + return result.tolist() + return result + Index: /issm/trunk-jpl/src/m/qmu/importancefactors.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/importancefactors.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/importancefactors.py (revision 23095) @@ -0,0 +1,63 @@ +import numpy as np +from MatlabFuncs import * + +def importancefactors(md,variablename,responsename): + '''IMPORTANCEFACTORS - compute importance factors for a certain variable and response. + +Usage: + factors=importancefactors(md,variablename,responsename) + + Example: factors=importancefactors(md,'drag','max_vel') +''' + + variablenamelength=len(variablename) + + #go through all response functions and find the one corresponding to the correct responsename + responsefunctions=md.qmu.results.dresp_out + found=-1 + for i in range(len(responsefunctions)): + if strcmpi(responsefunctions[i].descriptor,responsename): + found=i + break + if found < 0: + raise RuntimeError('importancefactors error message: could not find correct response function') + + responsefunctions=responsefunctions[found] + nfun=np.size(responsefunctions.var) + + #Now recover response to the correct design variable + importancefactors=[] + count=0 + for i in range(nfun): + desvar=responsefunctions.var[i] + if strncmpi(desvar,variablename,variablenamelength): + importancefactors.append(responsefunctions.impfac[i]) + count=count+1 + + if count==0: + raise RuntimeError('importancefactors error message: either response does not exist, or importancefactors are empty') + + importancefactors = np.array(importancefactors) + + if count==1: #we have scalar + factors=importancefactors + return factors + else: + #distribute importance factor + factors=importancefactors[(md.qmu.partition.conj().T).flatten().astype(int)] + #md.qmu.partition was created to index "c" style + + #weight importancefactors by area + #if numel(factors)==md.mesh.numberofvertices, + # #get areas for each vertex. + # aire=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y) + # num_elements_by_node=md.nodeconnectivity(:,) + # grid_aire=zeros(md.mesh.numberofvertices,1) + # for i=1:md.mesh.numberofvertices, + # for j=1:num_elements_by_node(i), + # grid_aire(i)=grid_aire(i)+aire(md.nodeconnectivity(i,j)) + # + # + # factors=factors./grid_aire + # + return factors Index: /issm/trunk-jpl/src/m/qmu/lclist_write.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/lclist_write.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/lclist_write.py (revision 23095) @@ -0,0 +1,67 @@ +import numpy as np +#move this later +from helpers import * + +from vector_write import * + +from linear_equality_constraint import * +from linear_inequality_constraint import * + +def lclist_write(fidi,cstring,cstring2,dvar): + ''' +function to write linear constraint list +''' + if dvar == None: + return + + # put linear constraints into lists for writing + + nvar=0 + pmatrix=[] + plower =[] + pupper =[] + ptarget=[] + pstype =[] + pscale =[] + + fnames=fieldnames(dvar) + for i in range(np.size(fnames)): + nvar=nvar+np.size(vars(dvar)[fnames[i]]) + pmatrix=[pmatrix,prop_matrix(vars(dvar)[fnames[i]])] + plower =[plower ,prop_lower(vars(dvar)[fnames[i]]) ] + pupper =[pupper ,prop_upper(vars(dvar)[fnames[i]]) ] + ptarget=[ptarget,prop_target(vars(dvar)[fnames[i]])] + pstype =[pstype ,prop_stype(vars(dvar)[fnames[i]]) ] + pscale =[pscale ,prop_scale(vars(dvar)[fnames[i]]) ] + + + # write linear constraints + + print ' Writing '+str(nvar)+' '+cstring+' linear constraints.' + + if len(pmatrix) != 0: + fidi.write('\t '+cstring2+'_matrix =\n') + vector_write(fidi,'\t ',pmatrix,6,76) + + if len(plower) != 0: + fidi.write('\t '+cstring2+'_lower_bounds =\n') + vector_write(fidi,'\t ',plower ,6,76) + + if len(pupper) != 0: + fidi.write('\t '+cstring2+'_upper_bounds =\n') + vector_write(fidi,'\t ',pupper ,6,76) + + if len(ptarget) != 0: + fidi.write('\t '+cstring2+'_targets =\n') + vector_write(fidi,'\t ',ptarget,6,76) + + if len(pstype) != 0: + fidi.write('\t '+cstring2+'_scale_types =\n') + vector_write(fidi,'\t ',pstype ,6,76) + + if len(pscale) != 0: + fidi.write('\t '+cstring2+'_scales =\n') + vector_write(fidi,'\t ',pscale ,6,76) + + + Index: /issm/trunk-jpl/src/m/qmu/param_write.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/param_write.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/param_write.py (revision 23095) @@ -0,0 +1,27 @@ +#move this later: +from helpers import * + +def param_write(fidi,sbeg,pname,smid,s,params): + ''' +function to write a parameter +''' + if not isfield(params,pname): + print 'WARNING: param_write:param_not_found: Parameter '+str(pname)+' not found in structure.' + return + + params_pname = vars(params)[pname] + + if type(params_pname) == bool and (not params_pname): + return + + if type(params_pname) == bool: + fidi.write(sbeg+str(pname)+s) + + elif type(params_pname) in [str]: + fidi.write(sbeg+str(pname)+smid+params_pname+s) + + elif type(params_pname) in [int, float]: + fidi.write(sbeg+str(pname)+smid+str(params_pname)+s) + + + Index: /issm/trunk-jpl/src/m/qmu/postqmu.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/postqmu.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/postqmu.py (revision 23095) @@ -0,0 +1,60 @@ +from os import system,getpid,stat +from os.path import isfile +from subprocess import Popen + +from dakota_out_parse import * +from helpers import * + +def postqmu(md,qmufile,qmudir='qmu'+str(getpid())): + ''' +Deal with dakota output results in files. + +INPUT function + md=postqmu(md,qmufile,qmudir) + +By default: qmudir = 'qmu'+pid (eg. 'qmu2189') +''' + + # check to see if dakota returned errors in the err file + qmuerrfile=str(md.miscellaneous.name)+'.qmu.err' + + if isfile(qmuerrfile) and stat(qmuerrfile).st_size > 0: + with open(qmuerrfile,'r') as fide: + fline=fide.read() + print fline + + raise RuntimeError('Dakota returned error in '+str(qmuerrfile)+' file. '+str(qmudir)+' directory retained.') + + # parse inputs and results from dakota + qmuinfile=str(md.miscellaneous.name)+'.qmu.in' + qmuoutfile=str(md.miscellaneous.name)+'.qmu.out' + + # unused and unimplemented + #[method,dvar,dresp_in]=dakota_in_parse(qmuinfile) + #dakotaresults.method =method + #dakotaresults.dvar =dvar + #dakotaresults.dresp_in =dresp_in + + [method,dresp_out,scm,pcm,srcm,prcm]=dakota_out_parse(qmuoutfile) + dakotaresults = struct() + dakotaresults.dresp_out=dresp_out + dakotaresults.scm =scm + dakotaresults.pcm =pcm + dakotaresults.srcm =srcm + dakotaresults.prcm =prcm + + if isfile('dakota_tabular.dat'): + # only need a subset of the outputs; dakota_out_parse handles .dat seperately + [method,dresp_dat,_,_,_,_]=dakota_out_parse('dakota_tabular.dat') + dakotaresults.dresp_dat=dresp_dat + + # put dakotaresults in their right location. + md.results.dakota=dakotaresults + + # move all the individual function evalutations into zip files + if not md.qmu.isdakota: + Popen('zip -mq params.in.zip params.in.[1-9]*',shell=True) + Popen('zip -mq results.out.zip results.out.[1-9]*',shell=True) + Popen('zip -mq matlab.out.zip matlab*.out.[1-9]*',shell=True) + + return md Index: /issm/trunk-jpl/src/m/qmu/preqmu.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/preqmu.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/preqmu.py (revision 23095) @@ -0,0 +1,142 @@ +import os +from MatlabFuncs import * +from expandvariables import * +from expandresponses import * +from helpers import * +from dakota_in_data import * +from process_qmu_response_data import * + +def preqmu(md,options): + '''QMU - apply Quantification of Margins and Uncertainties techniques + to a solution sequence (like stressbalance.py, progonstic.py, etc ...), + using the Dakota software from Sandia. + + options come from the solve.py routine. They can include Dakota options: + + qmudir: any directory where to run the qmu analysis + qmufile: input file for Dakota + + (ivap, iresp, imethod, and iparams are currently unimplemented) + 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 before analysis +''' + + print 'preprocessing dakota inputs' + qmudir = options.getfieldvalue('qmudir','qmu'+str(os.getpid())) + # qmudir = ['qmu_' datestr(now,'yyyymmdd_HHMMSS')] + qmufile = options.getfieldvalue('qmufile','qmu') + # qmufile cannot be changed unless ????script.sh is also changed + overwrite = options.getfieldvalue('overwrite','n') + ivar = options.getfieldvalue('ivar',0) + iresp = options.getfieldvalue('iresp',0) + imethod = options.getfieldvalue('imethod',0) + iparams = options.getfieldvalue('iparams',0) + + # first create temporary directory in which we will work + if strncmpi(overwrite,'y',1): + os.system('rm -rf '+qmudir+'/*') + else: + # does the directory exist? if so, then error out + if os.path.isdir(qmudir): + raise RuntimeError('Existing '+str(options.qmudir)+' directory, cannot overwrite. Specify "overwrite","y" option in solve arguments.') + + + os.makedirs(qmudir) + os.chdir(qmudir) + + # when running in library mode, the in file needs to be called md.miscellaneous.name.qmu.in + qmufile=md.miscellaneous.name + + # retrieve variables and resposnes for this particular analysis. + #print type(md.qmu.variables) + #print md.qmu.variables.__dict__ + #print ivar + variables=md.qmu.variables#[ivar] + responses=md.qmu.responses#[iresp] + + # expand variables and responses + #print variables.__dict__ + #print responses.__dict__ + variables=expandvariables(md,variables) + responses=expandresponses(md,responses) + + # go through variables and responses, and check they don't have more than + # md.qmu.numberofpartitions values. Also determine numvariables and numresponses + #[[[ + numvariables=0 + variable_fieldnames=fieldnames(variables) + for i in range(len(variable_fieldnames)): + field_name=variable_fieldnames[i] + fieldvariables=vars(variables)[field_name] + for j in range(np.size(fieldvariables)): + if strncmpi(fieldvariables[j].descriptor,'\'scaled_',8) and str2int(fieldvariables[j].descriptor,'last')>md.qmu.numberofpartitions: + raise RuntimeError('preqmu error message: one of the expanded variables has more values than the number of partitions (setup in md.qmu.numberofpartitions)') + + numvariables=numvariables+np.size(vars(variables)[field_name]) + + numresponses=0 + response_fieldnames=fieldnames(responses) + for i in range(len(response_fieldnames)): + field_name=response_fieldnames[i] + fieldresponses=vars(responses)[field_name] + for j in range(np.size(fieldresponses)): + if strncmpi(fieldresponses[j].descriptor,'\'scaled_',8) and str2int(fieldresponses[j].descriptor,'last')>md.qmu.numberofpartitions: + raise RuntimeError('preqmu error message: one of the expanded responses has more values than the number of partitions (setup in md.qmu.numberofpartitions)') + + numresponses=numresponses+np.size(vars(responses)[field_name]) + + #]]] + + # create in file for dakota + #dakota_in_data(md.qmu.method[imethod],variables,responses,md.qmu.params[iparams],qmufile) + dakota_in_data(md.qmu.method,variables,responses,md.qmu.params,qmufile) + +#====================================================================================# + #REMOVED FOR DEBUGGING ONLY: + #os.system('rm -rf '+str(md.miscellaneous.name)+'.py') +#====================================================================================# + + # build a list of variables and responses descriptors. the list is not expanded. + #[[[ + variabledescriptors=[] + # variable_fieldnames=fieldnames(md.qmu.variables[ivar]) + variable_fieldnames=fieldnames(md.qmu.variables) + for i in range(len(variable_fieldnames)): + field_name=variable_fieldnames[i] + #fieldvariables=vars(md.qmu.variables[ivar])[field_name] + fieldvariables=vars(md.qmu.variables)[field_name] + if type(fieldvariables) in [list,np.ndarray]: + for j in range(np.size(fieldvariables)): + variabledescriptors.append(fieldvariables[j].descriptor) + else: + variabledescriptors.append(fieldvariables.descriptor) + + + responsedescriptors=[] + # response_fieldnames=fieldnames(md.qmu.responses[iresp]) + response_fieldnames=fieldnames(md.qmu.responses) + for i in range(len(response_fieldnames)): + field_name=response_fieldnames[i] + #fieldresponses=vars(md.qmu.responses[iresp])[field_name] + fieldresponses=vars(md.qmu.responses)[field_name] + for j in range(np.size(fieldresponses)): + responsedescriptors.append(fieldresponses[j].descriptor) + + + #]]] + + # register the fields that will be needed by the Qmu model. + md.qmu.numberofresponses=numresponses + md.qmu.variabledescriptors=variabledescriptors + md.qmu.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) + + return md Index: /issm/trunk-jpl/src/m/qmu/process_qmu_response_data.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/process_qmu_response_data.m (revision 23094) +++ /issm/trunk-jpl/src/m/qmu/process_qmu_response_data.m (revision 23095) @@ -43,8 +43,8 @@ %ok, process the domains named in qmu_mass_flux_profiles, to build a list of segments (MatArray) md.qmu.mass_flux_segments=cell(num_mass_flux,1); - + md.qmu.mass_flux_segments for i=1:num_mass_flux, md.qmu.mass_flux_segments{i}=MeshProfileIntersection(md.mesh.elements,md.mesh.x,md.mesh.y,[md.qmu.mass_flux_profile_directory '/' md.qmu.mass_flux_profiles{i}]); end - + md.qmu.mass_flux_segments end Index: /issm/trunk-jpl/src/m/qmu/process_qmu_response_data.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/process_qmu_response_data.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/process_qmu_response_data.py (revision 23095) @@ -0,0 +1,51 @@ +from MatlabFuncs import * +import numpy as np +from MeshProfileIntersection import * + +from helpers import empty_nd_list + +def process_qmu_response_data(md): + ''' +PROCESS_QMU_RESPONSE_DATA - process any data necessary for the solutions to process the data. + + Usage: md=process_qmu_response_data(md) + + See also PREQMU, PRESOLVE +''' + + # preliminary data + process_mass_flux_profiles=0 + num_mass_flux=0 + + # loop through response descriptors, and act accordingly + for i in range(np.size(md.qmu.responsedescriptors)): + + # Do we have to process mass flux profiles? + if strncmpi(md.qmu.responsedescriptors[i],'indexed_MassFlux',16): + num_mass_flux+=1 + process_mass_flux_profiles=1 + + # deal with mass flux profiles + if process_mass_flux_profiles: + # we need a profile of points on which to compute the mass_flux, is it here? + if type(md.qmu.mass_flux_profiles) == float and np.isnan(md.qmu.mass_flux_profiles): + raise RuntimeError('process_qmu_response_data error message: could not find a mass_flux exp profile!') + + if type(md.qmu.mass_flux_profiles) != list: + raise RuntimeError('process_qmu_response_data error message: qmu_mass_flux_profiles field should be a list of domain outline names') + + if np.size(md.qmu.mass_flux_profiles) == 0: + raise RuntimeError('process_qmu_response_data error message: qmu_mass_flux_profiles cannot be empty!') + + if num_mass_flux!=np.size(md.qmu.mass_flux_profiles): + raise RuntimeError('process_qmu_response_data error message: qmu_mass_flux_profiles should be of the same size as the number of MassFlux responses asked for in the Qmu analysis') + + # ok, process the domains named in qmu_mass_flux_profiles, + # to build a list of segments (MatArray) + md.qmu.mass_flux_segments = empty_nd_list((num_mass_flux,1)) + + for i in range(num_mass_flux): + md.qmu.mass_flux_segments[i]=np.array(MeshProfileIntersection(md.mesh.elements,md.mesh.x,md.mesh.y,md.qmu.mass_flux_profile_directory+'/'+md.qmu.mass_flux_profiles[i])[0]) + + return md + Index: /issm/trunk-jpl/src/m/qmu/rlev_write.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/rlev_write.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/rlev_write.py (revision 23095) @@ -0,0 +1,107 @@ +import numpy as np + +#move this later +from helpers import * + +from vector_write import * +from param_write import * +#import relevent qmu classes + +from MatlabArray import * + +def rlevi_write(fidi,ltype,levels): + ''' + function to each type of response level +''' + + fidi.write('\t num_'+str(ltype)+' =') + levels = np.array(levels) + + if len(levels) > 0 and type(levels[0]) in [list,np.ndarray]: + for i in range(len(levels)): + fidi.write(' ' + str(len(levels[i]))) + else: + fidi.write(' ' + str(len(levels))) + + fidi.write('\n') + + fidi.write('\t '+str(ltype)+' =\n') + + # check if we have a vector of vectors, or just 1 vector + if np.size(levels) > 0 and type(levels[0]) in [list,np.ndarray]: + for i in range(len(levels)): + if len(levels[i]) != 0: + vector_write(fidi,'\t ',levels[i],8,76) + else: + vector_write(fidi,'\t ',levels,8,76) + + return + +def rlev_write(fidi,dresp,cstring,params): + ''' + function to write response levels +''' + from response_function import * + + if len(dresp) == 0 or len(fieldnames(dresp[0])) == 0: + return + + if type(dresp) in [list,np.ndarray]: + if len(dresp) > 0 and type(dresp[0]) == struct: + func = eval(cstring) + else: + func = type(dresp[0]) + elif type(dresp) == struct: + # type is defined within struct's contents + func = None + dresp = [dresp] + else: + func = type(dresp) + dresp = [dresp] + + # put responses into lists for writing + + nresp = 0 + respl = [] + probl = [] + rell = [] + grell = [] + + # assume all fields in dvar[0:n] are consistent (ex. all are normal_uncertain) + # which will always be true since this is called per field + fnames=fieldnames(dresp[0]) + for j in range(len(dresp)): + for i in range(np.size(fnames)): + if func == None: + func = type(vars(dresp[j])[fnames[i]]) + + nresp+=1 + [respli,probli,relli,grelli]=func.prop_levels([vars(dresp[j])[fnames[i]]]) + respl.extend(respli) + probl.extend(probli) + rell.extend(relli) + grell.extend(grelli) + + + # write response levels + + respl=allempty(respl) + probl=allempty(probl) + rell =allempty(rell) + grell=allempty(grell) + + param_write(fidi,'\t ','distribution',' ','\n',params) + if len(respl) != 0: + rlevi_write(fidi,'response_levels',respl) + param_write(fidi,'\t ','compute',' ','\n',params) + + if len(probl) != 0: + rlevi_write(fidi,'probability_levels',probl) + + if len(rell) != 0: + rlevi_write(fidi,'reliability_levels',rell) + + if len(grell) != 0: + rlevi_write(fidi,'gen_reliability_levels',grell) + + return Index: /issm/trunk-jpl/src/m/qmu/rlist_write.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/rlist_write.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/rlist_write.py (revision 23095) @@ -0,0 +1,94 @@ +import numpy as np + +#move this later +from helpers import * + +from vector_write import * +from response_function import * + +def rlist_write(fidi,cstring,cstring2,dresp,rdesc): + ''' + function to write response list +''' + + if dresp == None: + return + + func = eval(cstring) + + if type(dresp) not in [list,np.ndarray]: + dresp = [dresp] + + # put responses into lists for writing + # (and accumulate descriptors into list for subsequent writing) + + nresp=0 + pstype =[] + pscale =[] + pweight=[] + plower =[] + pupper =[] + ptarget=[] + + # assume all fields in dvar[0:n] are consistent (ex. all are normal_uncertain) + # which will always be true since this is called per field + fnames=fieldnames(dresp[0]) + for j in range(len(dresp)): + for i in range(np.size(fnames)): + nresp=nresp+np.size(vars(dresp[j])[fnames[i]]) + pstype.extend( func.prop_stype( vars(dresp[j])[fnames[i]])) + pscale.extend( func.prop_scale( vars(dresp[j])[fnames[i]])) + pweight.extend(func.prop_weight(vars(dresp[j])[fnames[i]])) + plower.extend( func.prop_lower( vars(dresp[j])[fnames[i]])) + pupper.extend( func.prop_upper( vars(dresp[j])[fnames[i]])) + ptarget.extend(func.prop_target(vars(dresp[j])[fnames[i]])) + rdesc.extend( func.prop_desc( vars(dresp[j])[fnames[i]],fnames[i])) + + + # write responses + + print ' Writing '+str(nresp)+' '+cstring+' responses.' + + if strcmp(cstring,'calibration_terms')==1: + fidi.write('\t'+cstring+' = '+str(nresp)+'\n') + + else: + fidi.write('\tnum_'+cstring+'s = '+str(nresp)+'\n') + + if not isempty(pstype): + fidi.write('\t '+cstring2+'_scale_types =\n') + vector_write(fidi,'\t ',pstype ,6,76) + + if not isempty(pscale): + fidi.write('\t '+cstring2+'_scales =\n') + vector_write(fidi,'\t ',pscale ,6,76) + + if not isempty(pweight): + if cstring2 == 'objective_function': + fidi.write('\t multi_objective_weights =\n') + vector_write(fidi,'\t ',pweight,6,76) + elif cstring2 == 'least_squares_term': + fidi.write('\t least_squares_weights =\n') + vector_write(fidi,'\t ',pweight,6,76) + + if not isempty(plower): + fidi.write('\t '+cstring2+'_lower_bounds =\n') + vector_write(fidi,'\t ',plower ,6,76) + + if not isempty(pupper): + fidi.write('\t '+cstring2+'_upper_bounds =\n') + vector_write(fidi,'\t ',pupper ,6,76) + + if not isempty(ptarget): + fidi.write('\t '+cstring2+'_targets =\n') + vector_write(fidi,'\t ',ptarget,6,76) + + # because qmu in files need '' for strings + for i in range(len(rdesc)): + if type(rdesc[i]) in [list,np.ndarray]: + for j in range(len(rdesc[i])): + rdesc[i][j] = "'" + rdesc[i][j] + "'" + else: + rdesc[i] = "'" + rdesc[i] + "'" + + return rdesc Index: /issm/trunk-jpl/src/m/qmu/setupdesign/QmuSetupResponses.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/setupdesign/QmuSetupResponses.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/setupdesign/QmuSetupResponses.py (revision 23095) @@ -0,0 +1,23 @@ +from MatlabFuncs import * +from copy import deepcopy + +def QmuSetupResponses(md,dresp,responses): + + #get descriptor + descriptor=responses.descriptor + + #decide whether this is a distributed response, which will drive whether we expand it into npart values, + #or if we just carry it forward as is. + + #ok, key off according to type of descriptor: + if strncmp(descriptor,'scaled_',7): + #we have a scaled response, expand it over the partition. + #ok, dealing with semi-discrete distributed response. Distribute according to how many + #partitions we want + for j in range(md.qmu.numberofpartitions): + dresp.append(deepcopy(responses)) + dresp[-1].descriptor=str(responses.descriptor)+'_'+str(j+1) + else: + dresp.append(responses) + + return dresp Index: /issm/trunk-jpl/src/m/qmu/setupdesign/QmuSetupVariables.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/setupdesign/QmuSetupVariables.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/setupdesign/QmuSetupVariables.py (revision 23095) @@ -0,0 +1,79 @@ +from MatlabFuncs import * +from uniform_uncertain import* +from normal_uncertain import * +from copy import deepcopy + +def QmuSetupVariables(md,dvar,variables): + + #get descriptor + descriptor=variables.descriptor + + #decide whether this is a distributed variable, which will drive whether we expand it into npart values, + #or if we just carry it forward as is. + + #ok, key off according to type of descriptor: + if strncmp(descriptor,'scaled_',7): + #we have a scaled variable, expand it over the partition. + + if isinstance(variables,uniform_uncertain): + if ((type(variables.lower) in [list,np.ndarray] and len(variables.lower) > md.qmu.numberofpartitions) or (type(variables.upper) in [list,np.ndarray] and len(variables.upper) > md.qmu.numberofpartitions)): + raise RuntimeError('QmuSetupDesign error message: upper and lower should be either a scalar or a "npart" length vector') + + elif isinstance(variables,normal_uncertain): + if type(variables.stddev) in [list,np.ndarray] and len(variables.stddev) > md.qmu.numberofpartitions: + raise RuntimeError('QmuSetupDesign error message: stddev should be either a scalar or a "npart" length vector') + + #ok, dealing with semi-discrete distributed variable. Distribute according to how many + #partitions we want + + for j in range(md.qmu.numberofpartitions): + dvar.append(deepcopy(variables)) + # "'" is because qmu.in files need for strings to be in actual '' + # must also account for whether we are given 1 instance or an array of instances + + # handle descriptors for everything + if type(dvar[-1].descriptor) in [list,np.ndarray] and len(variables.descriptor) > 1 and len(variables.upper) != md.qmu.numberofpartitions: + if type(variables.descriptor) == np.ndarray: + dvar[-1].descriptor = np.append(dvar[-1].descriptor,"'"+str(variables.descriptor)+'_'+str(j+1)+"'") + else: + dvar[-1].descriptor.append("'"+str(variables.descriptor)+'_'+str(j+1)+"'") + else: + dvar[-1].descriptor = "'"+str(variables.descriptor)+'_'+str(j+1)+"'" + + # handle uniform_uncertain + if isinstance(variables,uniform_uncertain): + if type(variables.lower) in [list,np.ndarray] and len(variables.lower) > 1 and len(variables.upper) != md.qmu.numberofpartitions: + if type(variables.lower) == np.ndarray: + dvar[-1].lower = np.append(dvar[-1].lower, variables.lower[j]) + else: + dvar[-1].lower.append(variables.lower[j]) + else: + dvar[-1].lower = variables.lower + + if type(variables.upper) in [list,np.ndarray] and len(variables.upper) > 1 and len(variables.upper) != md.qmu.numberofpartitions: + if type(variables.upper) == np.ndarray: + dvar[-1].upper = np.append(dvar[-1].upper, variables.upper[j]) + else: + dvar[-1].upper.append(variables.upper[j]) + else: + dvar[-1].upper = variables.upper + + # handle normal_uncertain + elif isinstance(variables,normal_uncertain): + if type(variables.stddev) in [list,np.ndarray] and len(variables.stddev) > 1 and len(variables.upper) != md.qmu.numberofpartitions: + if type(variables.stddev) == np.ndarray: + dvar[-1].stddev = np.append(dvar[-1].stddev, variables.stddev[j]) + else: + dvar[-1].stddev.append(variables.stddev[j]) + else: + dvar[-1].stddev = variables.stddev + + # running with a single instance, and therefore length 1 arrays of qmu classes + else: + dvar.append(variables) + # text parsing in dakota requires literal "'identifier'" not just "identifier" + for v in dvar: + v.descriptor = "'"+str(v.descriptor)+"'" + + return dvar + Index: /issm/trunk-jpl/src/m/qmu/vector_write.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/vector_write.m (revision 23094) +++ /issm/trunk-jpl/src/m/qmu/vector_write.m (revision 23095) @@ -1,5 +1,5 @@ - -%% function to write a vector on multiple lines - +% +% function to write a vector on multiple lines +% function []=vector_write(fidi,sbeg,vec,nmax,cmax) @@ -22,5 +22,4 @@ % assemble each line, flushing when necessary - for i=1:numel(vec) if isnumeric(vec(i)) Index: /issm/trunk-jpl/src/m/qmu/vector_write.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/vector_write.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/vector_write.py (revision 23095) @@ -0,0 +1,47 @@ +import numpy as np + +def vector_write(fidi,sbeg,vec,nmax,cmax): + ''' +function to write a vector on multiple lines +''' + if nmax == None: + nmax=np.inf + + if cmax == None: + cmax=np.inf + + # set up first iteration + svec =[] + nitem=nmax + lsvec=cmax + + # transpose vector from column-wise to row-wise + vec=np.array(vec).conj().T + + # assemble each line, flushing when necessary + for i in range(np.size(vec,0)): + + # [[1],[1],[1]...] should be [1,1,1,...] + if type(vec[i]) in [list,np.ndarray] and len(vec[i]) == 1: + sitem = str(vec[i][0]) + else: + sitem = str(vec[i]) + + nitem=nitem+1 + lsvec=lsvec+1+len(sitem) + + if (nitem <= nmax) and (lsvec <= cmax): + svec=str(svec)+' '+str(sitem) + else: + if len(svec) > 0: + fidi.write(str(svec)+'\n') + + svec=str(sbeg)+str(sitem) + nitem=1 + lsvec=len(svec) + + # flush buffer at , if necessary + if len(svec) > 0: + fidi.write(str(svec)+'\n') + + return Index: /issm/trunk-jpl/src/m/qmu/vlist_write.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/vlist_write.py (revision 23095) +++ /issm/trunk-jpl/src/m/qmu/vlist_write.py (revision 23095) @@ -0,0 +1,130 @@ +import numpy as np +#move this later +from helpers import * + +from vector_write import * + +from uniform_uncertain import * +from normal_uncertain import * + +def check(a,l,p): + '''in the event that a and b are equal, return a; +in the event that a and b are not equal, return their concatenation + + This is used for when both the input dvar and the 'cstring' variables have non-1 length +''' + + if np.size(a) == l: + if p == 0: + if type(a) in [list,np.ndarray]: + return a + else: + return [a] + else: + return [] + elif np.size(a) == 1: + if type(a) in [list,np.ndarray]: + return a + else: + return [a] + elif np.size(a) == 0: + return [] + else: + raise RuntimeError('ERROR vlist_write: input field had size '+str(np.size(a))+'; must have size of 0, 1, or match size of provided dvar ('+str(l)+').') + + return + +def vlist_write(fidi,cstring,cstring2,dvar): + ''' +function to write variable list +''' + if dvar == None: + return + from uniform_uncertain import * + func = eval(cstring) + + # put variables into lists for writing + + if type(dvar) not in [list,np.ndarray]: + dvar = [dvar] + + # assume all fields in dvar[0:n] are consistent (ex. all are normal_uncertain) + # which will always be true since this vlist_write is called per field + fnames=fieldnames(dvar[0]) + + nvar=0 + pinitpt=[[] for i in range(len(fnames))] + plower =[[] for i in range(len(fnames))] + pupper =[[] for i in range(len(fnames))] + pmean =[[] for i in range(len(fnames))] + pstddev=[[] for i in range(len(fnames))] + pinitst=[[] for i in range(len(fnames))] + pstype =[[] for i in range(len(fnames))] + pscale =[[] for i in range(len(fnames))] + pdesc =[[] for i in range(len(fnames))] + + for i in range(len(fnames)): + nvar += len(dvar) + for j in dvar: + j = vars(j)[fnames[i]] + pinitpt[i].extend(check(func.prop_initpt(j),len(dvar),len(pinitpt[i]))) + plower[i].extend(check(func.prop_lower(j),len(dvar),len(plower[i]))) + pupper[i].extend(check(func.prop_upper(j),len(dvar),len(pupper[i]))) + pmean[i].extend(check(func.prop_mean(j),len(dvar),len(pmean[i]))) + pstddev[i].extend(check(func.prop_stddev(j),len(dvar),len(pstddev[i]))) + pinitst[i].extend(check(func.prop_initst(j),len(dvar),len(pinitst[i]))) + pstype[i].extend(check(func.prop_stype(j),len(dvar),len(pstype[i]))) + pscale[i].extend(check(func.prop_scale(j),len(dvar),len(pscale[i]))) + pdesc[i].extend(check(func.prop_desc(j,fnames[i]),len(dvar),len(pdesc[i]))) + + pinitpt=allempty(pinitpt) + plower =allempty(plower) + pupper =allempty(pupper) + pmean =allempty(pmean) + pstddev=allempty(pstddev) + pinitst=allempty(pinitst) + pstype =allempty(pstype) + pscale =allempty(pscale) + pdesc =allempty(pdesc) + + # write variables + print ' Writing '+str(nvar)+' '+cstring+' variables.' + + fidi.write('\t'+cstring+' = '+str(nvar)+'\n') + if not isempty(pinitpt): + fidi.write('\t '+cstring2+'_initial_point =\n') + vector_write(fidi,'\t ',pinitpt,6,76) + + if not isempty(plower): + fidi.write('\t '+cstring2+'_lower_bounds =\n') + vector_write(fidi,'\t ',plower ,6,76) + + if not isempty(pupper): + fidi.write('\t '+cstring2+'_upper_bounds =\n') + vector_write(fidi,'\t ',pupper ,6,76) + + if not isempty(pmean): + fidi.write('\t '+cstring2+'_means =\n') + vector_write(fidi,'\t ',pmean ,6,76) + + if not isempty(pstddev): + fidi.write('\t '+cstring2+'_std_deviations =\n') + vector_write(fidi,'\t ',pstddev,6,76) + + if not isempty(pinitst): + fidi.write('\t '+cstring2+'_initial_state =\n') + vector_write(fidi,'\t ',pinitst,6,76) + + if not isempty(pstype): + fidi.write('\t '+cstring2+'_scale_types =\n') + vector_write(fidi,'\t ',pstype ,6,76) + + if not isempty(pscale): + fidi.write('\t '+cstring2+'_scales =\n') + vector_write(fidi,'\t ',pscale ,6,76) + + if not isempty(pdesc): + fidi.write('\t '+cstring2+'_descriptors =\n') + vector_write(fidi,'\t ',pdesc ,6,76) + + return Index: /issm/trunk-jpl/src/m/solve/WriteData.py =================================================================== --- /issm/trunk-jpl/src/m/solve/WriteData.py (revision 23094) +++ /issm/trunk-jpl/src/m/solve/WriteData.py (revision 23095) @@ -332,5 +332,4 @@ fid.write(struct.pack('i',len(data))) - #write each matrix: for matrix in data: if isinstance(matrix,(bool,int,long,float)): @@ -342,5 +341,6 @@ s=matrix.shape - if np.ndim(data) == 1: + + if np.ndim(matrix) == 1: fid.write(struct.pack('i',s[0])) fid.write(struct.pack('i',1)) Index: /issm/trunk-jpl/src/m/solve/loadresultsfromcluster.py =================================================================== --- /issm/trunk-jpl/src/m/solve/loadresultsfromcluster.py (revision 23094) +++ /issm/trunk-jpl/src/m/solve/loadresultsfromcluster.py (revision 23095) @@ -3,4 +3,6 @@ import platform from loadresultsfromdisk import loadresultsfromdisk + +from helpers import * def loadresultsfromcluster(md,runtimename=False): @@ -23,6 +25,6 @@ filelist.append(md.miscellaneous.name+'.qmu.err') filelist.append(md.miscellaneous.name+'.qmu.out') - if 'tabular_graphics_data' in md.qmu.params: - if md.qmu.params['tabular_graphics_data']: + if 'tabular_graphics_data' in fieldnames(md.qmu.params): + if md.qmu.params.tabular_graphics_data: filelist.append('dakota_tabular.dat') else: @@ -31,12 +33,21 @@ #If we are here, no errors in the solution sequence, call loadresultsfromdisk. - if os.path.getsize(md.miscellaneous.name+'.outbin')>0: - md=loadresultsfromdisk(md,md.miscellaneous.name+'.outbin') - else: - print 'WARNING, outbin file is empty for run '+md.miscellaneous.name + if os.path.exists(md.miscellaneous.name+'.outbin'): + if os.path.getsize(md.miscellaneous.name+'.outbin')>0: + md=loadresultsfromdisk(md,md.miscellaneous.name+'.outbin') + else: + print 'WARNING, outbin file is empty for run '+md.miscellaneous.name + elif not md.qmu.isdakota: + print 'WARNING, outbin file does not exist '+md.miscellaneous.name #erase the log and output files if md.qmu.isdakota: - filename=os.path.join('qmu'+str(os.getpid()),md.miscellaneous.name) + #filename=os.path.join('qmu'+str(os.getpid()),md.miscellaneous.name) + filename = md.miscellaneous.name + + # this will not work like normal as dakota doesn't generate outbin files, + # instead calls postqmu to store results directly in the model + # at md.results.dakota via dakota_out_parse + md=loadresultsfromdisk(md,md.miscellaneous.name+'.outbin') else: filename=md.miscellaneous.name @@ -52,5 +63,6 @@ if hostname==cluster.name: if md.qmu.isdakota: - filename=os.path.join('qmu'+str(os.getpid()),md.miscellaneous.name) + #filename=os.path.join('qmu'+str(os.getpid()),md.miscellaneous.name) + filename = md.miscellaneous.name TryRem('.queue',filename) else: @@ -62,10 +74,16 @@ TryRem('.bat',filename) + # remove this for bin file debugging TryRem('.bin',filename) + + #cwd = os.getcwd().split('/')[-1] + if md.qmu.isdakota: + os.chdir('..') + #TryRem('',cwd) return md def TryRem(extension,filename): - try: + try: os.remove(filename+extension) except OSError: Index: /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.py =================================================================== --- /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.py (revision 23094) +++ /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.py (revision 23095) @@ -3,4 +3,5 @@ from parseresultsfromdisk import parseresultsfromdisk import MatlabFuncs as m +from postqmu import * def loadresultsfromdisk(md,filename): @@ -57,6 +58,5 @@ #post processes qmu results if necessary else: - md=postqmu(md) - os.chdir('..') + md=postqmu(md,filename) return md Index: /issm/trunk-jpl/src/m/solve/marshall.py =================================================================== --- /issm/trunk-jpl/src/m/solve/marshall.py (revision 23094) +++ /issm/trunk-jpl/src/m/solve/marshall.py (revision 23095) @@ -11,6 +11,6 @@ marshall(md) """ - - print "marshalling file '%s.bin'." % md.miscellaneous.name + if md.verbose.solution: + print "marshalling file '%s.bin'." % md.miscellaneous.name #open file for binary writing Index: /issm/trunk-jpl/src/m/solve/solve.py =================================================================== --- /issm/trunk-jpl/src/m/solve/solve.py (revision 23094) +++ /issm/trunk-jpl/src/m/solve/solve.py (revision 23095) @@ -7,4 +7,6 @@ from waitonlock import waitonlock from loadresultsfromcluster import loadresultsfromcluster +from preqmu import * +#from MatlabFuncs import * def solve(md,solutionstring,*args): @@ -88,5 +90,6 @@ #check model consistency if options.getfieldvalue('checkconsistency','yes')=='yes': - print "checking model consistency" + if md.verbose.solution: + print "checking model consistency" ismodelselfconsistent(md) @@ -146,10 +149,11 @@ print 'The results must be loaded manually with md=loadresultsfromcluster(md).' else: #load results - print 'loading results from cluster' + if md.verbose.solution: + print 'loading results from cluster' md=loadresultsfromcluster(md) #post processes qmu results if necessary if md.qmu.isdakota: - if not strncmpi(options['keep'],'y',1): + if not strncmpi(options.getfieldvalue('keep','y'),'y',1): shutil.rmtree('qmu'+str(os.getpid()))