Index: /issm/trunk-jpl/src/m/classes/stochasticforcing.py =================================================================== --- /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 26634) +++ /issm/trunk-jpl/src/m/classes/stochasticforcing.py (revision 26635) @@ -1,8 +1,7 @@ import numpy as np - from checkfield import checkfield from fielddisplay import fielddisplay -from MatlabFuncs import * from WriteData import WriteData + class stochasticforcing(object): @@ -15,14 +14,14 @@ def __init__(self, *args): # {{{ self.isstochasticforcing = 0 - self.fields = np.nan - self.defaultdimension = 0 - self.default_id = np.nan - self.covariance = np.nan - self.randomflag = 1 + self.fields = np.nan + self.defaultdimension = 0 + self.default_id = np.nan + self.covariance = np.nan + self.randomflag = 1 if len(args) == 0: self.setdefaultparameters() else: - error('constructor not supported') + raise RuntimeError('constructor not supported for stochasticforcing') def __repr__(self): # {{{ @@ -43,7 +42,7 @@ def setdefaultparameters(self): # {{{ # Type of stabilization used - self.isstochasticforcing = 0 # stochasticforcing is turned off by default - self.fields = [] # Need to initialize to list to avoid "RuntimeError: object of type 'float' has no len()" on import of class - self.randomflag = 1 # true randomness is implemented by default + self.isstochasticforcing = 0 # stochasticforcing is turned off by default + self.fields = [] # Need to initialize to list to avoid "RuntimeError: object of type 'float' has no len()" on import of class + self.randomflag = 1 # true randomness is implemented by default return self #}}} @@ -54,47 +53,45 @@ return md - num_fields = numel(self.fields) + num_fields = len(self.fields) - #Check that covariance matrix is positive definite - try: - np.linalg.cholesky(self.covariance); - except: - error('md.stochasticforcing.covariance is not positive definite'); + #Check that covariance matrix is positive definite this is done internaly by linalg + np.linalg.cholesky(self.covariance) # Check that all fields agree with the corresponding md class checkdefaults = False for field in self.fields: - if (contains(field, 'SMB')): - if not (type(md.smb) == field): - error('md.smb does not agree with stochasticforcing field {}'.format(field)) - if (contains(field, 'frontalforcings')): - if not (type(md.frontalforcings) == field): - error('md.frontalforcings does not agree with stochasticforcing field {}'.format(field)) + if 'SMB' in field: + if type(md.smb).__name__ != field: + raise TypeError('md.smb does not agree with stochasticforcing field {}'.format(field)) + if 'frontalforcings' in field: + if type(md.frontalforcings).__name__ != field: + raise TypeError('md.frontalforcings does not agree with stochasticforcing field {}'.format(field)) #Checking for specific dimensions - if (field!='SMBautoregression' or field!='FrontalForcingsRignotAutoregression'): - checkdefaults = True #field with non-specific dimensionality + if not (field == 'SMBautoregression' or field == 'FrontalForcingsRignotAutoregression'): + checkdefaults = True #field with non-specific dimensionality #Retrieve sum of all the field dimensionalities - size_tot = self.defaultdimension*num_fields; - indSMBar = -1; #about to check for index of SMBautoregression - indTFar = -1; #about to check for index of FrontalForcingsRignotAutoregression + size_tot = self.defaultdimension * num_fields + indSMBar = -1 #about to check for index of SMBautoregression + indTFar = -1 #about to check for index of FrontalForcingsRignotAutoregression if ('SMBautoregression' in self.fields): - size_tot = size_tot-self.defaultdimension+md.smb.num_basins - indSMBar = np.where(self.fields=='SMBautoregression')[0][0] + size_tot = size_tot - self.defaultdimension + md.smb.num_basins + indSMBar = self.fields.index('SMBautoregression') if ('FrontalForcingsRignotAutoregression' in self.fields): - size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins - indSMBar = np.where(self.fields=='FrontalForcingsRignotAutoregression')[0][0] - if (indSMBar!=-1 and indTFar!=-1): - if((md.smb.ar_timestep!=md.frontalforcings.ar_timestep) and np.any(self.covariance[np.sum(self.dimensions[0:indSMBar]).astype(int):np.sum(self.dimensions[0:indSMBar+1]).astype(int),np.sum(self.dimensions[0:indTFar]).astype(int):np.sum(self.dimensions[0:indTFar+1]).astype(int)]!=0)): - error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance') + size_tot = size_tot - self.defaultdimension + md.frontalforcings.num_basins + indSMBar = self.fields.index('FrontalForcingsRignotAutoregression') + if (indSMBar != -1 and indTFar != -1): + covsum = self.covariance[np.sum(self.defaultdimensions[0:indSMBar]).astype(int):np.sum(self.defaultdimensions[0:indSMBar + 1]).astype(int), np.sum(self.defaultdimensions[0:indTFar]).astype(int):np.sum(self.defaultdimensions[0:indTFar + 1]).astype(int)] + if((md.smb.ar_timestep != md.frontalforcings.ar_timestep) and np.any(covsum != 0)): + raise IOError('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance') md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1]) md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', supportedstochforcings()) #md = checkfield(md, 'fieldname', 'stochasticforcing.dimensions', 'NaN', 1, 'Inf', 1, '>', 0, 'size', [num_fields]) # specific dimension for each field; NOTE: As opposed to MATLAB implementation, pass list - md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot]) # global covariance matrix + md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot]) # global covariance matrix md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1]) if (checkdefaults): md = checkfield(md, 'fieldname', 'stochasticforcing.defaultdimension', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) - md = checkfield(md, 'fieldname', 'stochasticforcing.default_id', 'Inf', 1, '>=', 0, '<=', self.defaultdimension, 'size', [md.mesh.numberofelements,1]) + md = checkfield(md, 'fieldname', 'stochasticforcing.default_id', 'Inf', 1, '>=', 0, '<=', self.defaultdimension, 'size', [md.mesh.numberofelements, 1]) return md # }}} @@ -114,33 +111,34 @@ else: #Retrieve dimensionality of each field - dimensions = self.defaultdimension*np.ones(num_fields) - for ind,field in enumerate(self.fields): + dimensions = self.defaultdimension * np.ones(num_fields) + for ind, field in enumerate(self.fields): #Checking for specific dimensions - if (field=='SMBautoregression'): + if (field == 'SMBautoregression'): dimensions[ind] = md.smb.num_basins - if (field=='FrontalForcingsRignotAutoregression'): + if (field == 'FrontalForcingsRignotAutoregression'): dimensions[ind] = md.frontalforcings.num_basins # Scaling covariance matrix (scale column-by-column and row-by-row) - scaledfields = ['DefaultCalving','SMBautoregression'] # list of fields that need scaling * 1/yts + scaledfields = ['DefaultCalving', 'SMBautoregression'] # list of fields that need scaling * 1/yts tempcovariance = np.copy(self.covariance) for i in range(num_fields): if self.fields[i] in scaledfields: - inds = range(int(np.sum(self.dimensions[0:i])), int(np.sum(self.dimensions[0:i + 1]))) - for row in inds: # scale rows corresponding to scaled field + inds = range(int(np.sum(dimensions[0:i])), int(np.sum(dimensions[0:i + 1]))) + for row in inds: # scale rows corresponding to scaled field tempcovariance[row, :] = 1 / yts * tempcovariance[row, :] - for col in inds: # scale columns corresponding to scaled field + for col in inds: # scale columns corresponding to scaled field tempcovariance[:, col] = 1 / yts * tempcovariance[:, col] #Set dummy default_id vector if defaults not used if np.isnan(self.default_id): - self.default_id = np.zeros(md.mesh.numberofelements) + self.default_id = np.zeros(md.mesh.numberofelements) WriteData(fid, prefix, 'data', num_fields, 'name', 'md.stochasticforcing.num_fields', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'fields', 'format', 'StringArray') - WriteData(fid, prefix, 'data', 'dimensions', 'name', 'md.stochasticforcing.dimensions', 'format', 'IntMat') - WriteData(fid, prefix, 'object', self, 'fieldname', 'default_id', 'format', 'IntMat') #12Nov2021 make sure this is zero-indexed! - WriteData(fid, prefix, 'object', self, 'fieldname', 'defaultdimension', 'format', 'Integer') + WriteData(fid, prefix, 'data', dimensions, 'name', 'md.stochasticforcing.dimensions', 'format', 'IntMat') + WriteData(fid, prefix, 'object', self, 'fieldname', 'default_id', 'format', 'IntMat') #12Nov2021 make sure this is zero-indexed! + WriteData(fid, prefix, 'object', self, 'fieldname', 'defaultdimension', 'format', 'Integer') WriteData(fid, prefix, 'data', tempcovariance, 'name', 'md.stochasticforcing.covariance', 'format', 'DoubleMat') WriteData(fid, prefix, 'object', self, 'fieldname', 'randomflag', 'format', 'Boolean') # }}} + def supportedstochforcings():