Changeset 27318 for issm/trunk-jpl/src/m/classes/SMBarma.py

Timestamp:

10/21/22 11:51:59 (3 years ago)

Author:

vverjans

Message:

NEW piecewise polynomials implemented in ARMA model schemes

File:

• issm/trunk-jpl/src/m/classes/SMBarma.py (modified) (6 diffs)

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

@@ -16 +16 @@
     def __init__(self, *args):  # {{{
         self.num_basins = 0
-        self.const = np.nan
-        self.trend = np.nan
+        self.num_params = 0
+        self.num_breaks = 0
+        self.polynomialparams = np.nan
         self.ar_order = 0
-        self.arma_initialtime = 0
-        self.arma_timestep = 0
         self.arlag_coefs = np.nan
         self.malag_coefs = np.nan
@@ -42 +41 @@
         s += '{}\n'.format(fielddisplay(self, 'num_basins', 'number of different basins [unitless]'))
         s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]'))
-        s += '{}\n'.format(fielddisplay(self, 'const', 'basin-specific constant values [m ice eq./yr]'))
-        s += '{}\n'.format(fielddisplay(self, 'trend', 'basin-specific trend values [m ice eq. yr^(-2)]'))
+        s += '{}\n'.format(fielddisplay(self, 'num_breaks', 'number of different breakpoints in the piecewise-polynomial (separating num_breaks+1 periods)'))
+        s += '{}\n'.format(fielddisplay(self, 'num_params', 'number of different parameters in the piecewise-polynomial (1:intercept only, 2:with linear trend, 3:with quadratic trend, etc.)'))
+        s += '{}\n'.format(fielddisplay(self, 'polynomialparams', 'coefficients for the polynomial (const,trend,quadratic,etc.),dim1 for basins,dim2 for periods,dim3 for orders, ex: polyparams=cat(num_params,intercepts,trendlinearcoefs,trendquadraticcoefs)'))
+        s += '{}\n'.format(fielddisplay(self, 'datebreaks', 'dates at which the breakpoints in the piecewise polynomial occur (1 row per basin) [yr]'))
         s += '{}\n'.format(fielddisplay(self, 'ar_order', 'order of the autoregressive model [unitless]'))
         s += '{}\n'.format(fielddisplay(self, 'ma_order', 'order of the moving-average model [unitless]'))
-        s += '{}\n'.format(fielddisplay(self, 'arma_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]'))
         s += '{}\n'.format(fielddisplay(self, 'arma_timestep', 'time resolution of the ARMA model [yr]'))
         s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of AR lag coefficients [unitless]'))
@@ -83 +83 @@
             self.arlag_coefs = np.zeros((self.num_basins, self.ar_order)) # Autoregression coefficients all set to 0
             print('      smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0')
-        if self.arma_initialtime == 0:
-            self.arma_initialtime = md.timestepping.start_time # ARMA model has no prescribed initial time
-            print('      smb.arma_initialtime (initial time in the ARMA model parameterization) not specified: set to md.timestepping.start_time')
         if self.arma_timestep == 0:
             self.arma_timestep = md.timestepping.time_step # ARMA model has no prescribed time step
@@ -100 +97 @@
     def checkconsistency(self, md, solution, analyses):  # {{{
         if 'MasstransportAnalysis' in analyses:
+            nbas = md.smb.num_basins;
+            nprm = md.smb.num_params;
+            nbrk = md.smb.num_breaks;
             md = checkfield(md, 'fieldname', 'smb.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0)
+            md = checkfield(md, 'fieldname', 'smb.num_params', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0)
+            md = checkfield(md, 'fieldname', 'smb.num_breaks', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0)
             md = checkfield(md, 'fieldname', 'smb.basin_id', 'Inf', 1, '>=', 0, '<=', md.smb.num_basins, 'size', [md.mesh.numberofelements])
-            if len(np.shape(self.const)) == 1:
-                self.const = np.array([self.const])
-                self.trend = np.array([self.trend])
-            md = checkfield(md, 'fieldname', 'smb.const', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) # Scheme fails if passed as column vector
-            md = checkfield(md, 'fieldname', 'smb.trend', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) # Scheme fails if passed as column vector; NOTE: As opposed to MATLAB implementation, pass list
+            if len(np.shape(self.polynomialparams)) == 1:
+                self.polynomialparams = np.array([[self.polynomialparams]])
+            if(nbas>1 and nbrk>=1 and nprm>1):
+                md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm)
+            elif(nbas==1):
+                md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nprm,nbrk+1],'numel',nbas*(nbrk+1)*nprm)
+            elif(nbrk==0):
+                md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nprm],'numel',nbas*(nbrk+1)*nprm)
+            elif(nprm==1):
+                md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk],'numel',nbas*(nbrk+1)*nprm)
             md = checkfield(md, 'fieldname', 'smb.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0)
-            md = checkfield(md, 'fieldname', 'smb.arma_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1)
             md = checkfield(md, 'fieldname', 'smb.arma_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep
             md = checkfield(md, 'fieldname', 'smb.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ar_order])
             md = checkfield(md, 'fieldname', 'smb.malag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ma_order])
-            
+            if(nbrk>0):
+                md = checkfield(md, 'fieldname', 'smb.datebreaks', 'NaN', 1, 'Inf', 1, 'size', [nbas,nbrk])
+            elif(np.size(md.smb.datebreaks)==0 or np.all(np.isnan(md.smb.datebreaks))):
+                pass
+            else:
+                raise RuntimeError('md.smb.num_breaks is 0 but md.smb.datebreaks is not empty')
+
             if(np.any(np.isnan(self.refelevation) is False) or np.size(self.refelevation) > 1):
                 if len(np.shape(self.refelevation)) == 1:
@@ -149 +161 @@
     def marshall(self, prefix, md, fid):  # {{{
         yts = md.constants.yts
-
+        nbas = md.smb.num_basins;
+        nprm = md.smb.num_params;
+        nper = md.smb.num_breaks+1;
         templapserates    = np.copy(md.smb.lapserates)
         tempelevationbins = np.copy(md.smb.elevationbins)
         temprefelevation  = np.copy(md.smb.refelevation)
+        # Scale the parameters #
+        polyparamsScaled   = np.copy(md.smb.polynomialparams)
+        polyparams2dScaled = np.zeros((nbas,nper*nprm))
+        if(nprm>1):
+            # Case 3D #
+            if(nbas>1 and nper>1):
+                for ii in range(nprm):
+                    polyparamsScaled[:,:,ii] = polyparamsScaled[:,:,ii]*(1/yts)**(ii+1)
+                # Fit in 2D array #
+                for ii in range(nprm):
+                    polyparams2dScaled[:,ii*nper:(ii+1)*nper] = 1*polyparamsScaled[:,:,ii]
+            # Case 2D and higher-order params at increasing row index #
+            elif(nbas==1):
+                for ii in range(nprm):
+                    polyparamsScaled[ii,:] = polyparamsScaled[ii,:]*(1/yts)**(ii+1)
+                # Fit in row array #
+                for ii in range(nprm):
+                    polyparams2dScaled[0,ii*nper:(ii+1)*nper] = 1*polyparamsScaled[ii,:]
+            # Case 2D and higher-order params at incrasing column index #
+            elif(nper==1):
+                for ii in range(nprm):
+                    polyparamsScaled[:,ii] = polyparamsScaled[:,ii]*(1/yts)**(ii+1)
+                # 2D array is already in correct format #
+                polyparams2dScaled = np.copy(polyparamsScaled)
+        else:
+            polyparamsScaled   = polyparamsScaled*(1/yts)
+            # 2D array is already in correct format #
+            polyparams2dScaled = np.copy(polyparamsScaled)
+        
+        if(nper==1):
+            dbreaks = np.zeros((nbas,1))
+        else:
+            dbreaks = np.copy(md.smb.datebreaks)
+
         if(np.any(np.isnan(md.smb.lapserates))):
             templapserates = np.zeros((md.smb.num_basins,2))
@@ -172 +220 @@
         WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 13, 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'num_basins', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'num_breaks', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'num_params', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ar_order', 'format', 'Integer')
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ma_order', 'format', 'Integer')
-        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'arma_initialtime', 'format', 'Double', 'scale', yts)
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'arma_timestep', 'format', 'Double', 'scale', yts)
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'basin_id', 'data', self.basin_id - 1, 'name', 'md.smb.basin_id', 'format', 'IntMat', 'mattype', 2)  # 0-indexed
-        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'const', 'format', 'DoubleMat', 'name', 'md.smb.const', 'scale', 1 / yts, 'yts', yts)
-        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'trend', 'format', 'DoubleMat', 'name', 'md.smb.trend', 'scale', 1 / (yts ** 2), 'yts', yts)
+        WriteData(fid, prefix, 'data', polyparams2dScaled, 'name', 'md.smb.polynomialparams', 'format', 'DoubleMat')
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.smb.arlag_coefs', 'yts', yts)
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'malag_coefs', 'format', 'DoubleMat', 'name', 'md.smb.malag_coefs', 'yts', yts)
+        WriteData(fid, prefix, 'data', dbreaks, 'name', 'md.smb.datebreaks', 'format', 'DoubleMat','scale',yts)
         WriteData(fid, prefix, 'data', templapserates, 'name', 'md.smb.lapserates', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts)
         WriteData(fid, prefix, 'data', tempelevationbins, 'name', 'md.smb.elevationbins', 'format', 'DoubleMat')