Changeset 27251


Ignore:
Timestamp:
08/31/22 06:32:32 (3 years ago)
Author:
vverjans
Message:

CHG Changing autoregression schemes towards ARMA schemes: work is ongoing (forgot two files)

Location:
issm/trunk-jpl/src/m/classes
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r27165 r27251  
    5252                        for field=self.fields
    5353                                %Checking agreement of classes
    54                                 if(contains(field,'SMBautoregression'))
     54                                if(contains(field,'SMBarma'))
    5555                                        mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field))));
    5656                                        if~(isequal(class(md.smb),char(mdname)))
    … …  
    8888                                        end
    8989                                end
    90                                 if(contains(field,'BasalforcingsDeepwaterMeltingRateAutoregression'))
     90                                if(contains(field,'BasalforcingsDeepwaterMeltingRatearma'))
    9191                                        mdname = structstoch.mdnames(find(strcmp(structstoch.fields,char(field))));
    9292                                        if~(isequal(class(md.basalforcings),char(mdname)))
    … …  
    111111            end
    112112                                %Checking for specific dimensions
    113                                 if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression') || strcmp(field,'BasalforcingsDeepwaterMeltingRateAutoregression'))
     113                                if ~(strcmp(field,'SMBarma') || strcmp(field,'FrontalForcingsRignotarma') || strcmp(field,'BasalforcingsDeepwaterMeltingRatearma'))
    114114                                        checkdefaults = true; %field with non-specific dimensionality
    115115                                end
    … …  
    118118                        %Retrieve all the field dimensionalities
    119119                        dimensions = self.defaultdimension*ones(1,num_fields);
    120                         indSMBar   = -1; %about to check for index of SMBautoregression
    121                         indTFar   = -1; %about to check for index of FrontalForcingsRignotAutoregression
    122                         indBDWar          = -1; %about to check for index of BasalforcingsDeepwaterMeltingRateAutoregression
    123                         if any(contains(self.fields,'SMBautoregression'))
    124                                 indSMBar = find(contains(self.fields,'SMBautoregression')); %index of SMBar, now check for consistency with other ar timesteps
    125                                 dimensions(indSMBar) = md.smb.num_basins;
    126                                 if(md.smb.ar_timestep<self.stochastictimestep)
    127                                         error('SMBautoregression cannot have a timestep shorter than stochastictimestep');
    128                                 end
    129                         end
    130                         if any(contains(self.fields,'FrontalForcingsRignotAutoregression'))
    131                                 indTFar = find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with other ar timesteps
    132                                 dimensions(indTFar) = md.frontalforcings.num_basins;
    133                                 if(md.frontalforcings.ar_timestep<self.stochastictimestep)
    134                                         error('FrontalForcingsRignotAutoregression cannot have a timestep shorter than stochastictimestep');
    135                                 end
    136                         end
    137                         if any(contains(self.fields,'BasalforcingsDeepwaterMeltingRateAutoregression'))
    138                                 indBDWar        = find(contains(self.fields,'BasalforcingsDeepwaterMeltingRateAutoregression')); %index of BDWar, now check for consistency with other ar timesteps
    139                                 dimensions(indBDWar) = md.basalforcings.num_basins;
    140                                 if(md.basalforcings.ar_timestep<self.stochastictimestep)
    141                                         error('BasalforcingsDeepwaterMeltingRateAutoregression cannot have a timestep shorter than stochastictimestep');
     120                        indSMBarma = -1; %about to check for index of SMBarma
     121                        indTFarma  = -1; %about to check for index of FrontalForcingsRignotarma
     122                        indBDWarma  = -1; %about to check for index of BasalforcingsDeepwaterMeltingRatearma
     123                        if any(contains(self.fields,'SMBarma'))
     124                                indSMBarma = find(contains(self.fields,'SMBarma')); %index of SMBarma, now check for consistency with other arma timesteps
     125                                dimensions(indSMBarma) = md.smb.num_basins;
     126                                if(md.smb.arma_timestep<self.stochastictimestep)
     127                                        error('SMBarma cannot have a timestep shorter than stochastictimestep');
     128                                end
     129                        end
     130                        if any(contains(self.fields,'FrontalForcingsRignotarma'))
     131                                indTFarma       = find(contains(self.fields,'FrontalForcingsRignotarma')); %index of TFarma, now check for consistency with other arma timesteps
     132                                dimensions(indTFarma) = md.frontalforcings.num_basins;
     133                                if(md.frontalforcings.arma_timestep<self.stochastictimestep)
     134                                        error('FrontalForcingsRignotarma cannot have a timestep shorter than stochastictimestep');
     135                                end
     136                        end
     137                        if any(contains(self.fields,'BasalforcingsDeepwaterMeltingRatearma'))
     138                                indBDWarma      = find(contains(self.fields,'BasalforcingsDeepwaterMeltingRatearma')); %index of BDWarma, now check for consistency with other arma timesteps
     139                                dimensions(indBDWarma) = md.basalforcings.num_basins;
     140                                if(md.basalforcings.arma_timestep<self.stochastictimestep)
     141                                        error('BasalforcingsDeepwaterMeltingRatearma cannot have a timestep shorter than stochastictimestep');
    142142                                end
    143143                        end
    144144                        size_tot = sum(dimensions);
    145145
    146                         if(indSMBar~=-1 && indTFar~=-1) %both autoregressive models are used: check autoregressive time step consistency
    147                                 if(md.smb.ar_timestep~=md.frontalforcings.ar_timestep)
    148                                         crossentries = reshape(self.covariance(1+sum(dimensions(1:indSMBar-1)):sum(dimensions(1:indSMBar)),1+sum(dimensions(1:indTFar-1)):sum(dimensions(1:indTFar))),1,[]);
     146                        if(indSMBarma~=-1 && indTFarma~=-1) %both ARMA models are used: check ARMA time step consistency
     147                                if(md.smb.arma_timestep~=md.frontalforcings.arma_timestep)
     148                                        crossentries = reshape(self.covariance(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma))),1,[]);
    149149                                        if any(crossentries~=0)
    150                                                 error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance');
    151                                         end
    152                                 end
    153                         end
    154                         if(indSMBar~=-1 && indBDWar~=-1) %both autoregressive models are used: check autoregressive time step consistency
    155                                 if(md.smb.ar_timestep~=md.basalforcings.ar_timestep)
    156                                         crossentries = reshape(self.covariance(1+sum(dimensions(1:indSMBar-1)):sum(dimensions(1:indSMBar)),1+sum(dimensions(1:indBDWar-1)):sum(dimensions(1:indBDWar))),1,[]);
     150                                                error('SMBarma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance');
     151                                        end
     152                                end
     153                        end
     154                        if(indSMBarma~=-1 && indBDWarma~=-1) %both ARMA models are used: check ARMA time step consistency
     155                                if(md.smb.arma_timestep~=md.basalforcings.arma_timestep)
     156                                        crossentries = reshape(self.covariance(1+sum(dimensions(1:indSMBarma-1)):sum(dimensions(1:indSMBarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]);
    157157                                        if any(crossentries~=0)
    158                                                 error('SMBautoregression and BasalforcingsDeepwaterMeltingRateAutoregression have different ar_timestep and non-zero covariance');
    159                                         end
    160                                 end
    161                         end
    162                         if(indTFar~=-1 && indBDWar~=-1) %both autoregressive models are used: check autoregressive time step consistency
    163                                 if(md.frontalforcings.ar_timestep~=md.basalforcings.ar_timestep)
    164                                         crossentries = reshape(self.covariance(1+sum(dimensions(1:indTFar-1)):sum(dimensions(1:indTFar)),1+sum(dimensions(1:indBDWar-1)):sum(dimensions(1:indBDWar))),1,[]);
     158                                                error('SMBarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance');
     159                                        end
     160                                end
     161                        end
     162                        if(indTFarma~=-1 && indBDWarma~=-1) %both ARMA models are used: check ARMA time step consistency
     163                                if(md.frontalforcings.arma_timestep~=md.basalforcings.arma_timestep)
     164                                        crossentries = reshape(self.covariance(1+sum(dimensions(1:indTFarma-1)):sum(dimensions(1:indTFarma)),1+sum(dimensions(1:indBDWarma-1)):sum(dimensions(1:indBDWarma))),1,[]);
    165165                                        if any(crossentries~=0)
    166                                                 error('FrontalForcingsRignotAutoregression and BasalforcingsDeepwaterMeltingRateAutoregression have different ar_timestep and non-zero covariance');
     166                                                error('FrontalForcingsRignotarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance');
    167167                                        end
    168168                                end
    … …  
    181181                        disp(sprintf('   stochasticforcing parameters:'));
    182182                        fielddisplay(self,'isstochasticforcing','is stochasticity activated?');
    183                         fielddisplay(self,'fields','fields with stochasticity applied, ex: [{''SMBautoregression''}], or [{''SMBforcing''},{''DefaultCalving''}]');
     183                        fielddisplay(self,'fields','fields with stochasticity applied, ex: [{''SMBarma''}], or [{''SMBforcing''},{''DefaultCalving''}]');
    184184                        fielddisplay(self,'defaultdimension','dimensionality of the noise terms (does not apply to fields with their specific dimension)');
    185185                        fielddisplay(self,'default_id','id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)');
    … …  
    212212                                for field=self.fields
    213213                                        %Checking for specific dimensions
    214                                         if(strcmp(field,'SMBautoregression'))
     214                                        if(strcmp(field,'SMBarma'))
    215215                                                dimensions(ind) = md.smb.num_basins;
    216216                                        end
    217                                         if(strcmp(field,'FrontalForcingsRignotAutoregression'))
     217                                        if(strcmp(field,'FrontalForcingsRignotarma'))
    218218                                                dimensions(ind) = md.frontalforcings.num_basins;
    219219                                        end
    220                                         if(strcmp(field,'BasalforcingsDeepwaterMeltingRateAutoregression'))
     220                                        if(strcmp(field,'BasalforcingsDeepwaterMeltingRatearma'))
    221221                                                dimensions(ind) = md.basalforcings.num_basins;
    222222                                        end
    … …  
    225225
    226226                                %Scaling covariance matrix (scale column-by-column and row-by-row)
    227                                 scaledfields = {'BasalforcingsDeepwaterMeltingRateAutoregression','BasalforcingsSpatialDeepwaterMeltingRate','DefaultCalving','FloatingMeltRate','SMBautoregression','SMBforcing'}; %list of fields that need scaling *1/yts
     227                                scaledfields = {'BasalforcingsDeepwaterMeltingRatearma','BasalforcingsSpatialDeepwaterMeltingRate','DefaultCalving','FloatingMeltRate','SMBarma','SMBforcing'}; %list of fields that need scaling *1/yts
    228228                                tempcovariance = self.covariance; %copy of covariance to avoid writing back in member variable
    229229                                for i=1:num_fields
    … …  
    266266        % supported and corresponding md names
    267267        structure.fields = {...
    268                 'BasalforcingsDeepwaterMeltingRateAutoregression',...
     268                'BasalforcingsDeepwaterMeltingRatearma',...
    269269                'BasalforcingsSpatialDeepwaterMeltingRate',...
    270270                'DefaultCalving',...
    … …  
    273273                'FrictionCoulombWaterPressure',...
    274274                'FrictionSchoofWaterPressure',...
    275                 'FrontalForcingsRignotAutoregression',...
    276                 'SMBautoregression',...
     275                'FrontalForcingsRignotarma',...
     276                'SMBarma',...
    277277                'SMBforcing'
    278278                };
    279279        structure.mdnames = {...
    280                 'autoregressionlinearbasalforcings',...
     280                'linearbasalforcingsarma',...
    281281                'spatiallinearbasalforcings',...
    282282                'calving',...
    … …  
    285285                'frictioncoulomb',...
    286286                'frictionschoof',...
    287                 'frontalforcingsrignotautoregression',...
    288                 'SMBautoregression',...
     287                'frontalforcingsrignotarma',...
     288                'SMBarma',...
    289289                'SMBforcing'
    290290        };

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

    r27165 r27251  
    4040        s += '   FloatingMeltRate\n'
    4141        s += '   FrictionWaterPressure\n'
    42         s += '   FrontalForcingsRignotAutoregression (thermal forcing)\n'
    43         s += '   SMBautoregression\n'
     42        s += '   FrontalForcingsRignotarma (thermal forcing)\n'
     43        s += '   SMBarma\n'
    4444        s += '   SMBforcing\n'
    4545        return s
    … …  
    7575        for field in self.fields:
    7676            # Checking agreement of classes
    77             if 'SMBautoregression' in field:
     77            if 'SMBarma' in field:
    7878                mdname = structstoch[field]
    7979                if (type(md.smb).__name__ != mdname):
    … …  
    9999                if (type(md.basalforcings).__name__ != mdname):
    100100                    raise TypeError('md.basalforcings does not agree with stochasticforcing field {}'.format(field))
    101             if 'BasalforcingsDeepwaterMeltingRateAutoregression' in field:
     101            if 'BasalforcingsDeepwaterMeltingRatearma' in field:
    102102                mdname = structstoch[field]
    103103                if (type(md.basalforcings).__name__ != mdname):
    … …  
    115115
    116116            # Checking for specific dimensions
    117             if field not in['SMBautoregression', 'FrontalForcingsRignotAutoregression','BasalforcingsDeepwaterMeltingRateAutoregression']:
     117            if field not in['SMBarma', 'FrontalForcingsRignotarma','BasalforcingsDeepwaterMeltingRatearma']:
    118118                checkdefaults = True  # field with non-specific dimensionality
    119119
    120120        # Retrieve sum of all the field dimensionalities
    121121        dimensions = self.defaultdimension*np.ones((num_fields))
    122         indSMBar   = -1  # About to check for index of SMBautoregression
    123         indTFar    = -1  # About to check for index of FrontalForcingsRignotAutoregression
    124         indBDWar   = -1  # About to check for index of BasalforcingsDeepwaterMeltingRateAutoregression
    125         if ('SMBautoregression' in self.fields):
    126             indSMBar = self.fields.index('SMBautoregression')  # Index of SMBar, now check for consistency with other timesteps
    127             dimensions[indSMBar] = md.smb.num_basins
    128             if(md.smb.ar_timestep<self.stochastictimestep):
    129                 raise TypeError('SMBautoregression cannot have a timestep shorter than stochastictimestep')
    130         if ('FrontalForcingsRignotAutoregression' in self.fields):
    131             indTFar = self.fields.index('FrontalForcingsRignotAutoregression')  # Index of TFar, now check for consistency with other timesteps
    132             dimensions[indTFar] = md.frontalforcings.num_basins
    133             if(md.frontalforcings.ar_timestep<self.stochastictimestep):
    134                 raise TypeError('FrontalForcingsRignotAutoregression cannot have a timestep shorter than stochastictimestep')
    135         if ('BasalforcingsDeepwaterMeltingRateAutoregression' in self.fields):
    136             indBDWar = self.fields.index('BasalforcingsDeepwaterMeltingRateAutoregression')  # Index of BDWar, now check for consistency with other timesteps
    137             dimensions[indTFar] = md.basalforcings.num_basins
    138             if(md.basalforcings.ar_timestep<self.stochastictimestep):
    139                 raise TypeError('BasalforcingsDeepwaterMeltingRateAutoregression cannot have a timestep shorter than stochastictimestep')
     122        indSMBarma   = -1  # About to check for index of SMBarma
     123        indTFarma    = -1  # About to check for index of FrontalForcingsRignotarma
     124        indBDWarma   = -1  # About to check for index of BasalforcingsDeepwaterMeltingRatearma
     125        if ('SMBarma' in self.fields):
     126            indSMBarma = self.fields.index('SMBarma')  # Index of SMBarma, now check for consistency with other timesteps
     127            dimensions[indSMBarma] = md.smb.num_basins
     128            if(md.smb.arma_timestep<self.stochastictimestep):
     129                raise TypeError('SMBarma cannot have a timestep shorter than stochastictimestep')
     130        if ('FrontalForcingsRignotarma' in self.fields):
     131            indTFarma = self.fields.index('FrontalForcingsRignotarma')  # Index of TFarma, now check for consistency with other timesteps
     132            dimensions[indTFarma] = md.frontalforcings.num_basins
     133            if(md.frontalforcings.arma_timestep<self.stochastictimestep):
     134                raise TypeError('FrontalForcingsRignotarma cannot have a timestep shorter than stochastictimestep')
     135        if ('BasalforcingsDeepwaterMeltingRatearma' in self.fields):
     136            indBDWarma = self.fields.index('BasalforcingsDeepwaterMeltingRatearma')  # Index of BDWarma, now check for consistency with other timesteps
     137            dimensions[indTFarma] = md.basalforcings.num_basins
     138            if(md.basalforcings.arma_timestep<self.stochastictimestep):
     139                raise TypeError('BasalforcingsDeepwaterMeltingRatearma cannot have a timestep shorter than stochastictimestep')
    140140        size_tot = np.sum(dimensions)
    141141
    142         if (indSMBar != -1 and indTFar != -1):  # Both autoregressive models are used: check autoregressive time step consistency
    143             covsum = self.covariance[np.sum(dimensions[0:indSMBar]).astype(int):np.sum(dimensions[0:indSMBar + 1]).astype(int), np.sum(dimensions[0:indTFar]).astype(int):np.sum(dimensions[0:indTFar + 1]).astype(int)]
    144             if((md.smb.ar_timestep != md.frontalforcings.ar_timestep) and np.any(covsum != 0)):
    145                 raise IOError('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance')
    146         if (indSMBar != -1 and indBDWar != -1):  # Both autoregressive models are used: check autoregressive time step consistency
    147             covsum = self.covariance[np.sum(dimensions[0:indSMBar]).astype(int):np.sum(dimensions[0:indSMBar + 1]).astype(int), np.sum(dimensions[0:indBDWar]).astype(int):np.sum(dimensions[0:indBDWar + 1]).astype(int)]
    148             if((md.smb.ar_timestep != md.basalforcings.ar_timestep) and np.any(covsum != 0)):
    149                 raise IOError('SMBautoregression and BasalforcingsDeepwaterMeltingRateAutoregression have different ar_timestep and non-zero covariance')
    150         if (indTFar != -1 and indBDWar != -1):  # Both autoregressive models are used: check autoregressive time step consistency
    151             covsum = self.covariance[np.sum(dimensions[0:indTFar]).astype(int):np.sum(dimensions[0:indTFar + 1]).astype(int), np.sum(dimensions[0:indBDWar]).astype(int):np.sum(dimensions[0:indBDWar + 1]).astype(int)]
    152             if((md.frontalforcings.ar_timestep != md.basalforcings.ar_timestep) and np.any(covsum != 0)):
    153                 raise IOError('FrontalForcingsRignotAutoregression and BasalforcingsDeepwaterMeltingRateAutoregression have different ar_timestep and non-zero covariance')
     142        if (indSMBarma != -1 and indTFarma != -1):  # Both ARMA models are used: check ARMA time step consistency
     143            covsum = self.covariance[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int)]
     144            if((md.smb.arma_timestep != md.frontalforcings.arma_timestep) and np.any(covsum != 0)):
     145                raise IOError('SMBarma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance')
     146        if (indSMBarma != -1 and indBDWarma != -1):  # Both ARMA models are used: check ARMA time step consistency
     147            covsum = self.covariance[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)]
     148            if((md.smb.arma_timestep != md.basalforcings.arma_timestep) and np.any(covsum != 0)):
     149                raise IOError('SMBarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance')
     150        if (indTFarma != -1 and indBDWarma != -1):  # Both ARMA models are used: check ARMA time step consistency
     151            covsum = self.covariance[np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)]
     152            if((md.frontalforcings.arma_timestep != md.basalforcings.arma_timestep) and np.any(covsum != 0)):
     153                raise IOError('FrontalForcingsRignotarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance')
    154154
    155155        md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1])
    … …  
    184184            for ind, field in enumerate(self.fields):
    185185                # Checking for specific dimensions
    186                 if (field == 'SMBautoregression'):
     186                if (field == 'SMBarma'):
    187187                    dimensions[ind] = md.smb.num_basins
    188                 if (field == 'FrontalForcingsRignotAutoregression'):
     188                if (field == 'FrontalForcingsRignotarma'):
    189189                    dimensions[ind] = md.frontalforcings.num_basins
    190                 if (field == 'BasalforcingsDeepwaterMeltingRateAutoregression'):
     190                if (field == 'BasalforcingsDeepwaterMeltingRatearma'):
    191191                    dimensions[ind] = md.basalforcings.num_basins
    192192
    193193            # Scaling covariance matrix (scale column-by-column and row-by-row)
    194             scaledfields = ['BasalforcingsDeepwaterMeltingRateAutoregression','BasalforcingsSpatialDeepwaterMeltingRate','DefaultCalving', 'FloatingMeltRate', 'SMBautoregression', 'SMBforcing']  # list of fields that need scaling * 1/yts
     194            scaledfields = ['BasalforcingsDeepwaterMeltingRatearma','BasalforcingsSpatialDeepwaterMeltingRate','DefaultCalving', 'FloatingMeltRate', 'SMBarma', 'SMBforcing']  # list of fields that need scaling * 1/yts
    195195            tempcovariance = np.copy(self.covariance)
    196196            for i in range(num_fields):
    … …  
    229229           supported and corresponding md names
    230230        """
    231         structure = {'BasalforcingsDeepwaterMeltingRateAutoregression': 'autoregressionlinearbasalforcings',
     231        structure = {'BasalforcingsDeepwaterMeltingRatearma': 'linearbasalforcingsarma',
    232232                     'BasalforcingsSpatialDeepwaterMeltingRate': 'spatiallinearbasalforcings',
    233233                     'DefaultCalving': 'calving',
    … …  
    236236                     'FrictionCoulombWaterPressure': 'frictioncoulomb',
    237237                     'FrictionSchoofWaterPressure': 'frictionschoof',
    238                      'FrontalForcingsRignotAutoregression': 'frontalforcingsrignotautoregression',
    239                      'SMBautoregression': 'SMBautoregression',
     238                     'FrontalForcingsRignotarma': 'frontalforcingsrignotarma',
     239                     'SMBarma': 'SMBarma',
    240240                     'SMBforcing': 'SMBforcing'}
    241241        return structure
Note: See TracChangeset for help on using the changeset viewer.