Changeset 26635

Timestamp:

11/18/21 01:08:12 (3 years ago)

Author:

bdef

Message:

BUG:fixing Stochasticforcing translation

File:

• issm/trunk-jpl/src/m/classes/stochasticforcing.py (modified) (5 diffs)

issm/trunk-jpl/src/m/classes/stochasticforcing.py

@@ -1 +1 @@
 import numpy as np
-
 from checkfield import checkfield
 from fielddisplay import fielddisplay
-from MatlabFuncs import *
 from WriteData import WriteData
+
 
 class stochasticforcing(object):
@@ -15 +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 +42 @@
     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 +53 @@
             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 +111 @@
         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():