diff --git a/lib/eofs/cdms.py b/lib/eofs/cdms.py index 9b7453b..24a8c76 100644 --- a/lib/eofs/cdms.py +++ b/lib/eofs/cdms.py @@ -659,3 +659,348 @@ def getWeights(self): """ return self._solver.getWeights() + + +class EEof(object): + """EEOF analysis (meta-data enabled `cdms2` interface)""" + + def __init__(self, dataset, lag, weights=None, center=True, ddof=1): + """Create an EEof object. + + The EEOF solution is computed at initialization time. Method + calls are used to retrieve computed quantities. + + **Argument:** + + *dataset* + A `cdms2` variable containing the data to be analysed. Time + must be the first dimension. Missing values are allowed + provided that they are constant with time (e.g., values of + an oceanographic field over land). + + *lag* + A no of lag window is used to reconstruct input matrix as + extended input matrix to analysis Eof on it. + window = lag + 1, so lag must be greater or equal to zero. + zero lag will produced same result as Eof, instead of EEof. + + **Optional arguments:** + + *weights* + Sets the weighting method. The following pre-defined + weighting methods are available: + + * *'area'* : Square-root of grid cell area normalized by + total grid area. Requires a latitude-longitude grid to be + present in the `cdms2` variable *dataset*. This is a + fairly standard weighting strategy. If you are unsure + which method to use and you have gridded data then this + should be your first choice. + + * *'coslat'* : Square-root of cosine of latitude. Requires a + latitude dimension to be present in the `cdms2` variable + *dataset*. + + * *None* : Equal weights for all grid points (*'none'* is + also accepted). + + Alternatively an array of weights whose shape is compatible + with the `cdms2` variable *dataset* may be supplied instead + of specifying a weighting method. + + *center* + If *True*, the mean along the first axis of *dataset* (the + time-mean) will be removed prior to analysis. If *False*, + the mean along the first axis will not be removed. Defaults + to *True* (mean is removed). + + The covariance interpretation relies on the input data being + anomalies with a time-mean of 0. Therefore this option + should usually be set to *True*. Setting this option to + *True* has the useful side effect of propagating missing + values along the time dimension, ensuring that a solution + can be found even if missing values occur in different + locations at different times. + + *ddof* + 'Delta degrees of freedom'. The divisor used to normalize + the covariance matrix is *N - ddof* where *N* is the + number of samples. Defaults to *1*. + + **Returns:** + + *solver* + An `EEof` instance. + + **Examples:** + + EEOF analysis with grid-cell-area weighting for the input field:: + + from eofs.cdms import EEof + solver = EEof(dataset, weights='area') + + Authors : Arulalan.T, Dileep.K + + """ + # Check that dataset is recognised by cdms2 as a variable. + if not cdms2.isVariable(dataset): + raise TypeError('the input data must be a cdms2 variable') + # Store the time axis as an instance variable. + self._timeax = dataset.getTime() + + # Verify that a time axis was found, getTime returns None when a + # time axis is not found. + if self._timeax is None: + raise ValueError('time axis not found') + + self.window = lag + 1 + if (lag < 0): + raise ValueError('lag window should not be less than 0') + elif (lag == 0): + print "lag is 0. So EEof result will be same as Eof" + self._lagtimeax = self._timeax + else: + # genearate lag time axis + # Remove last window length from original time axis + lag_time_axis = self._timeax[:-self.window+1] + lag_time_axis = cdms2.createAxis(lag_time_axis) + lag_time_axis.units = self._timeax.units + lag_time_axis.id = self._timeax.id + lag_time_axis.designateTime() + self._lagtimeax = lag_time_axis + + # Check the dimension order of the input, time must be the first + # dimension. + order = dataset.getOrder() + if order[0] != 't': + raise ValueError('time must be the first dimension, ' + 'consider using the reorder() method') + # Verify the presence of at least one spatial dimension. The + # instance variable channels will also be used as a partial axis + # list when constructing meta-data. It contains the spatial + # dimensions. + self._channels = dataset.getAxisList() + self._channels.remove(self._timeax) + if len(self._channels) < 1: + raise ValueError('one or more spatial dimensions are required') + # Store the missing value attribute of the data set in an + # instance variable so that it is recoverable later. + self._missing_value = dataset.getMissing() + # Generate an appropriate set of weights for the input dataset. There + # are several weighting schemes. The 'area' weighting scheme requires + # a latitude-longitude grid to be present, the 'cos_lat' scheme only + # requires a latitude dimension. + if weights in ('none', None): + # No weights requested, set the weight array to None. + wtarray = None + else: + try: + # Generate a weights array of the appropriate kind, with a + # shape compatible with the data set. + scheme = weights.lower() + wtarray = weights_array(dataset, scheme=scheme) + except AttributeError: + # Weights is not a string, assume it is an array. + wtarray = weights + except ValueError, err: + # Weights is not recognized, raise an error. + raise ValueError(err) + # Cast the wtarray to the same type as the dataset. This prevents the + # promotion of 32-bit input to 64-bit on multiplication with the + # weight array when not required. This will fail with a AttributeError + # exception if the weights array is None, which it may be if no + # weighting was requested. + try: + wtarray = wtarray.astype(dataset.dtype) + except AttributeError: + pass + # Create an EEofSolver object using appropriate arguments for this + # data set. The object will be used for the decomposition and + # for returning the results. + self._solver = standard.EEof(dataset.asma(), + lag=lag, + weights=wtarray, + center=center, + ddof=ddof) + #: Number of EEOFs in the solution. + self.neeofs = self._solver.neeofs + # name for the dataset. + self._dataset_name = cdms2_name(dataset).replace(' ', '_') + self._dataset_id = dataset.id + + def pcs(self, pcscaling=0, npcs=None): + """Principal component time series (PCs) of EEOF. + + **Optional arguments:** + + *pcscaling* + Set the scaling of the retrieved PCs. The following + values are accepted: + + * *0* : Un-scaled principal components (default). + * *1* : Principal components are scaled to unit variance + (divided by the square-root of their eigenvalue). + * *2* : Principal components are multiplied by the + square-root of their eigenvalue. + + *npcs* + Number of PCs to retrieve. Defaults to all the PCs. If the + number of requested PCs is more than the number that are + available, then all available PCs will be returned. + + **Returns:** + + *pcs* + A `cdms2` variable containing the ordered PCs. The PCs are + numbered from 0 to *npcs* - 1. + + **Examples:** + + All un-scaled PCs:: + + pcs = solver.pcs() + + First 3 PCs scaled to unit variance:: + + pcs = solver.pcs(npcs=3, pcscaling=1) + + """ + pcs = self._solver.pcs(pcscaling, npcs) + pcsax = cdms2.createAxis(range(pcs.shape[1]), id='pc') + pcsax.long_name = 'pc_number' + axlist = [self._lagtimeax, pcsax] + pcs = cdms2.createVariable(pcs, id='pcs', axes=axlist) + pcs.long_name = 'extended_eof_principal_components' + return pcs + + def eeofs(self, eeofscaling=0, neeofs=None): + """Extended Emipirical orthogonal functions (EEOFs). + + **Optional arguments:** + + *eofscaling* + Sets the scaling of the EEOFs. The following values are + accepted: + + * *0* : Un-scaled EOFs (default). + * *1* : EEOFs are divided by the square-root of their + eigenvalues. + * *2* : EEOFs are multiplied by the square-root of their + eigenvalues. + + *neeofs* + Number of EEOFs to return. Defaults to all EEOFs. If the + number of EEOFs requested is more than the number that are + available, then all available EEOFs will be returned. + + **Returns:** + + *eeofs* + A `cdms2` variable containing the ordered EEOFs. The EEOFs are + numbered from 0 to *neeofs* - 1. + + **Examples:** + + All EEOFs with no scaling:: + + eeofs = solver.eeofs() + + First 3 EEOFs with scaling applied:: + + eeofs = solver.eeofs(neeofs=3, eeofscaling=1) + + """ + eeofs = self._solver.eeofs(eofscaling, neofs) + eeofs.fill_value = self._missing_value + eeofax = cdms2.createAxis(range(len(eeofs)), id='eeof') + eeofax.long_name = 'eeof_number' + + lagax = cdms2.createAxis(range(self.window), id='lag') + + axlist = [eeofax, lagax] + self._channels + eeofs = cdms2.createVariable(eeofs, + id='eeofs', + axes=axlist, + fill_value=self._missing_value) + eeofs.long_name = 'extended_empirical_orthogonal_functions' + return eeofs + + def eigenvalues(self, neigs=None): + """Eigenvalues (decreasing variances) associated with each EEOF. + + **Optional argument:** + + *neigs* + Number of eigenvalues to return. Defaults to all + eigenvalues.If the number of eigenvalues requested is more + than the number that are available, then all available + eigenvalues will be returned. + + **Returns:** + + *eigenvalues* + A `cdms2` variable containing the eigenvalues arranged + largest to smallest. + + **Examples:** + + All eigenvalues:: + + eigenvalues = solver.eigenvalues() + + The first eigenvalue:: + + eigenvalue1 = solver.eigenvalues(neigs=1) + + """ + lambdas = self._solver.eigenvalues(neigs=neigs) + eofax = cdms2.createAxis(range(len(lambdas)), id='eigenvalue') + eofax.long_name = 'eigenvalue_number' + axlist = [eofax] + lambdas = cdms2.createVariable(lambdas, id='eigenvalues', axes=axlist) + lambdas.long_name = 'eigenvalues' + return lambdas + + def varianceFraction(self, neigs=None): + """Fractional EEOF variances. + + The fraction of the total variance explained by each EOF mode, + values between 0 and 1 inclusive. + + **Optional argument:** + + *neigs* + Number of eigenvalues to return the fractional variance for. + Defaults to all eigenvalues. If the number of eigenvalues + requested is more than the number that are available, then + fractional variances for all available eigenvalues will be + returned. + + **Returns:** + + *variance_fractions* + A `cdms2` variable containing the fractional variances for + each eigenvalue. The eigenvalues are numbered from 0 to + *neigs* - 1. + + **Examples:** + + The fractional variance represented by each eigenvalue:: + + variance_fractions = solver.varianceFraction() + + The fractional variance represented by the first 3 eigenvalues:: + + variance_fractions = solver.VarianceFraction(neigs=3) + + """ + vfrac = self._solver.varianceFraction(neigs=neigs) + eofax = cdms2.createAxis(range(len(vfrac)), id='eigenvalue') + eofax.long_name = 'eigenvalue_number' + axlist = [eofax] + vfrac = cdms2.createVariable(vfrac, id='variance_fractions', + axes=axlist) + vfrac.long_name = 'variance_fractions' + return vfrac + + diff --git a/lib/eofs/iris.py b/lib/eofs/iris.py index f43188d..969ce27 100644 --- a/lib/eofs/iris.py +++ b/lib/eofs/iris.py @@ -656,3 +656,147 @@ def getWeights(self): """ return self._solver.getWeights() + +class EEof(object): + """EEOF analysis (meta-data enabled `iris` interface)""" + + def __init__(self, cube, lag, weights=None, center=True, ddof=1): + """Create an Eof object. + + The EEOF solution is computed at initialization time. Method + calls are used to retrieve computed quantities. + + **Argument:** + + *dataset* + A `~iris.cube.Cube` instance containing the data to be + analysed. Time must be the first dimension. Missing values + are allowed provided that they are constant with time (e.g., + values of an oceanographic field over land). + + *lag* + A no of lag window is used to reconstruct input matrix as + extended input matrix to analysis Eof on it. + window = lag + 1, so lag must be greater or equal to zero. + zero lag will produced same result as Eof, instead of EEof. + + **Optional arguments:** + + *weights* + Sets the weighting method. The following pre-defined + weighting methods are available: + + * *'area'* : Square-root of grid cell area normalized by + total grid area. Requires a latitude-longitude grid to be + present in the `~iris.cube.Cube` *dataset*. This is a + fairly standard weighting strategy. If you are unsure + which method to use and you have gridded data then this + should be your first choice. + + * *'coslat'* : Square-root of cosine of latitude. Requires a + latitude dimension to be present in the `~iris.cube.Cube` + *dataset*. + + * *None* : Equal weights for all grid points (*'none'* is + also accepted). + + Alternatively an array of weights whose shape is compatible + with the `~iris.cube.Cube` *dataset* may be supplied instead + of specifying a weighting method. + + *center* + If *True*, the mean along the first axis of *dataset* (the + time-mean) will be removed prior to analysis. If *False*, + the mean along the first axis will not be removed. Defaults + to *True* (mean is removed). + + The covariance interpretation relies on the input data being + anomalies with a time-mean of 0. Therefore this option + should usually be set to *True*. Setting this option to + *True* has the useful side effect of propagating missing + values along the time dimension, ensuring that a solution + can be found even if missing values occur in different + locations at different times. + + *ddof* + 'Delta degrees of freedom'. The divisor used to normalize + the covariance matrix is *N - ddof* where *N* is the + number of samples. Defaults to *1*. + + **Returns:** + + *solver* + An `EEof` instance. + + **Examples:** + + EOF analysis with grid-cell-area weighting for the input field:: + + from eofs.iris import EEof + solver = EEof(cube, weights='area') + + """ + # Check that the input is an Iris cube. + if type(cube) is not Cube: + raise TypeError('the input must be an iris cube') + # Check for a time coordinate, raise an error if there isn't one. + # The coord_and_dim function will raise a ValuerError with a + # useful message so no need to handle it explicitly here. + self._time, self._time_dim = coord_and_dim(cube, 'time') + if self._time_dim != 0: + raise ValueError('time must be the first dimension, ' + 'consider using the transpose() method') + + self.window = lag + 1 + if (lag < 0): + raise ValueError('lag window should not be less than 0') + elif (lag == 0): + print "lag is 0. So EEof result will be same as Eof" + self._lagtimeax = self._time + else: + # genearate lag time axis + # Remove last window length from original time axis + lag_time_axis = self._time[:-self.window+1] + + # TODO : Kindly write proper syntax to create time axis in iris + # Above lag_time_axis variable must be iris support time axis with + # proper id, units, designateTime and all. + + self._lagtimeax = lag_time_axis + + + # Get the cube coordinates and remove time, leaving just the other + # dimensions. + self._coords = list(copy(cube.dim_coords)) + self._coords.remove(self._time) + if len(self._coords) < 1: + raise ValueError('one or more non-time dimensions are required') + # Define the weights array for the cube. + if weights is None: + wtarray = None + else: + try: + scheme = weights.lower() + wtarray = weights_array(cube, scheme=scheme) + except AttributeError: + wtarray = weights + try: + # Ensure weights are the same type as the cube data. + wtarray = wtarray.astype(cube.data.dtype) + except AttributeError: + pass + # Initialize a solver. + # Create an EEofSolver object using appropriate arguments for this + # data set. The object will be used for the decomposition and + # for returning the results. + self._solver = standard.EEof(cube.data, + lag=lag, + weights=wtarray, + center=center, + ddof=ddof) + #: Number of EOFs in the solution. + self.neeofs = self._solver.neeofs + # Get the name of the cube to refer to later. + self._cube_name = cube.name(default='dataset').replace(' ', '_') + self._cube_var_name = cube.var_name + diff --git a/lib/eofs/standard.py b/lib/eofs/standard.py index 58162b3..6324f21 100644 --- a/lib/eofs/standard.py +++ b/lib/eofs/standard.py @@ -22,6 +22,7 @@ import numpy as np import numpy.ma as ma +from numpy.lib.stride_tricks import as_strided from .tools.standard import correlation_map, covariance_map @@ -748,3 +749,400 @@ def getWeights(self): """ return self._weights + + +class EEof(Eof, object): + """Extended EOF analysis (EEOF) (`numpy` interface)""" + + def __init__(self, dataset, lag, weights=None, center=True, ddof=1): + """Create an Eof object. + + The EEOF solution is computed at initialization time. Method + calls are used to retrieve computed quantities. + + **Arguments:** + + *dataset* + A `numpy.ndarray` or `numpy.ma.MaskedArray` with two or more + dimensions containing the data to be analysed. The first + dimension is assumed to represent time. Missing values are + permitted, either in the form of a masked array, or + `numpy.nan` values. Missing values must be constant with time + (e.g., values of an oceanographic field over land). + + *lag* + A no of lag window is used to reconstruct input matrix as + extended input matrix to analysis Eof on it. + window = lag + 1, so lag must be greater or equal to zero. + zero lag will produced same result as Eof, instead of EEof. + + **Optional arguments:** + + *weights* + An array of weights whose shape is compatible with those of + the input array *dataset*. The weights can have the same + shape as *dataset* or a shape compatible with an array + broadcast (i.e., the shape of the weights can can match the + rightmost parts of the shape of the input array *dataset*). + If the input array *dataset* does not require weighting then + the value *None* may be used. Defaults to *None* (no + weighting). + + *center* + If *True*, the mean along the first axis of *dataset* (the + time-mean) will be removed prior to analysis. If *False*, + the mean along the first axis will not be removed. Defaults + to *True* (mean is removed). + + The covariance interpretation relies on the input data being + anomaly data with a time-mean of 0. Therefore this option + should usually be set to *True*. Setting this option to + *True* has the useful side effect of propagating missing + values along the time dimension, ensuring that a solution + can be found even if missing values occur in different + locations at different times. + + *ddof* + 'Delta degrees of freedom'. The divisor used to normalize + the covariance matrix is *N - ddof* where *N* is the + number of samples. Defaults to *1*. + + **Returns:** + + *solver* + An `EEof` instance. + + **Examples:** + + EEOF analysis with no weighting:: + + from eofs.standard import EEof + solver = EEof(data) + + EEOF analysis of a data array with spatial dimensions that + represent latitude and longitude with weighting. In this example + the data array is dimensioned (ntime, nlat, nlon), and in order + for the latitude weights to be broadcastable to this shape, an + extra length-1 dimension is added to the end:: + + from eofs.standard import EEof + import numpy as np + latitude = np.linspace(-90, 90, 73) + weights_array = np.cos(np.deg2rad(latitude))[:, np.newaxis] + solver = EEof(data, weights=weight_array) + + Authors : Arulalan.T, Dileep.K + + """ + # Store the input data in an instance variable. + if dataset.ndim < 2: + raise ValueError('the input data set must be ' + 'at least two dimensional') + self._data = dataset.copy() + # Check if the input is a masked array. If so fill it with NaN. + try: + self._data = self._data.filled(fill_value=np.nan) + self._filled = True + except AttributeError: + self._filled = False + # Store information about the shape/size of the input data. + self._records = self._data.shape[0] + self._originalshape = self._data.shape[1:] + self._lag = lag + self._window = lag + 1 + channels = np.product(self._originalshape) + # Weight the data set according to weighting argument. + if weights is not None: + try: + # The broadcast_arrays call returns a list, so the second index + # is retained, but also we want to remove the time dimension + # from the weights so the the first index from the broadcast + # array is taken. + self._weights = np.broadcast_arrays( + self._data[0:1], weights)[1][0] + self._data = self._data * self._weights + except ValueError: + raise ValueError('weight array dimensions are incompatible') + except TypeError: + raise TypeError('weights are not a valid type') + else: + self._weights = None + # Remove the time mean of the input data unless explicitly told + # not to by the "center" argument. + self._centered = center + if center: + self._data = self._center(self._data) + # Reshape to two dimensions (time, space) creating the design matrix. + self._data = self._data.reshape([self._records, channels]) + # Find the indices of values that are not missing in one row. All the + # rows will have missing values in the same places provided the + # array was centered. If it wasn't then it is possible that some + # missing values will be missed and the singular value decomposition + # will produce not a number for everything. + if not self._valid_nan(self._data): + raise ValueError('missing values detected in different ' + 'locations at different times') + + # Get covariance matrix of eeof by passing reverse time axis input data + cov_matrix_eeof = self._embed_dimension(self._data[::-1], self._window) + nonMissingIndex = np.where(np.isnan(cov_matrix_eeof[0]) == False)[0] + # Remove missing values from the design matrix. + dataNoMissing = cov_matrix_eeof[:, nonMissingIndex] + # make memory free + del cov_matrix_eeof + # new channels dimension + new_channels = self._window * channels + + # Compute the singular value decomposition of the design matrix. + try: + A, Lh, E = np.linalg.svd(dataNoMissing, full_matrices=False) + except (np.linalg.LinAlgError, ValueError): + raise ValueError('error encountered in SVD, check that missing ' + 'values are in the same places at each time and ' + 'that all the values are not missing') + # Singular values are the square-root of the eigenvalues of the + # covariance matrix. Construct the eigenvalues appropriately and + # normalize by N-ddof where N is the number of observations. This + # corresponds to the eigenvalues of the normalized covariance matrix. + self._ddof = ddof + normfactor = float(self._records - self._lag - self._ddof) + self._L = Lh * Lh / normfactor + # Store the number of eigenvalues (and hence EEOFs) that were actually + # computed. + self.neeofs = len(self._L) + # Re-introduce missing values into the eigenvectors in the same places + # as they exist in the input maps. Create an array of not-a-numbers + # and then introduce data values where required. We have to use the + # astype method to ensure the eigenvectors are the same type as the + # input dataset since multiplication by np.NaN will promote to 64-bit. + self._flatE = np.ones([self.neeofs, new_channels], + dtype=self._data.dtype) * np.NaN + self._flatE = self._flatE.astype(self._data.dtype) + self._flatE[:, nonMissingIndex] = E + # Remove the scaling on the principal component time-series that is + # implicitily introduced by using SVD instead of eigen-decomposition. + # The PCs may be re-scaled later if required. + self._P = A * Lh + + def _embed_dimension(self, array, window): + """ + Embed a given length window from the leading dimension of an array. + + **Arguments:** + + *array* + A 2-dimensional (nxm) `numpy.ndarray` or `numpy.ma.MaskedArray`. + + *window* + An integer specifying the length of the embedding window. + + **Returns:** + + A 2-dimenensional ((n-window+1) x (m*window)) `numpy.ndarray` or + `numpy.ma.MaskedArray` which is a view on the input *array*. + + **Example:** + + data = np.arange(4*3).reshape(4, 3) + >>> data + array([[ 0, 1, 2], + [ 3, 4, 5], + [ 6, 7, 8], + [ 9, 10, 11]]) + + >>> _embed_dimension(data, window=2) + array([[ 0, 1, 2, 3, 4, 5], + [ 3, 4, 5, 6, 7, 8], + [ 6, 7, 8, 9, 10, 11]]) + + >>> _embed_dimension(data, window=3) + array([[ 0, 1, 2, 3, 4, 5, 6, 7, 8], + [ 3, 4, 5, 6, 7, 8, 9, 10, 11]]) + + For general case, input array shape of (n, m), this + _embed_dimension returns shape of (n - window + 1, m * window), + where n is no of rows and m is no of columns. + + Note1 : For window is 1, then there would be no changes in + input data shape. i.e. For lag=0, EEof is equal to + normal Eof. + + >>> _embed_dimension(data, window=1) + array([[ 0, 1, 2], + [ 3, 4, 5], + [ 6, 7, 8], + [ 9, 10, 11]]) + + Note 2 : window = lag + 1 + + References : A.Hannachi, 2004, "A Primer for EOF Analysis of Climate + Data", Department of Meteorology, University of Reading Reading + RG6 6BB, U.K. (page numbers 15-28) + Link : http://eros.eas.gatech.edu/eas-8803/lectures/EOFs/eofprimer.pdf + + Author: Dr Andrew Dawson + Date: 18-11-2013 + + """ + if array.ndim != 2: + raise ValueError('array must have exactly 2 dimensions') + if window >= array.shape[0]: + raise ValueError('embedding window must be shorter than the ' + 'first dimension of the array') + n, m = array.shape + nwin = n - window + 1 + shape = (nwin, window) + array.shape[1:] + + strides = (array.strides[0], array.strides[0]) + array.strides[1:] + windowed = as_strided(array, shape=shape, strides=strides) + if ma.isMaskedArray(array): + if array.mask is ma.nomask: + windowed_mask = array.mask + else: + strides = ((array.mask.strides[0], array.mask.strides[0]) + + array.mask.strides[1:]) + windowed_mask = as_strided(array.mask, shape=shape, + strides=strides) + # end of if array.mask is ma.nomask: + windowed = ma.array(windowed, mask=windowed_mask) + # end of if ma.isMaskedArray(array): + out_shape = (nwin, window * array.shape[1]) + + return windowed.reshape(out_shape) + + def eeofs(self, eofscaling=0, neeofs=None): + """Extended Empirical orthogonal functions (EEOFs). + + **Optional arguments:** + + *eofscaling* + Sets the scaling of the EOFs. The following values are + accepted: + + * *0* : Un-scaled EEOFs (default). + * *1* : EEOFs are divided by the square-root of their + eigenvalues. + * *2* : EEOFs are multiplied by the square-root of their + eigenvalues. + + *neeofs* + Number of EEOFs to return. Defaults to all EEOFs. If the + number of EEOFs requested is more than the number that are + available, then all available EEOFs will be returned. + + **Returns:** + + *eofs* + An array with the ordered EEOFs along the first dimension. + + **Examples:** + + All EEOFs with no scaling:: + + eeofs = solver.eeofs() + + The leading EEOF with scaling applied:: + + eeof1 = solver.eeofs(neeofs=1, eofscaling=1) + + """ + if neeofs is None or neeofs > self.neeofs: + neeofs = self.neeofs + slicer = slice(0, neeofs) + neeofs = neeofs or self.neeofs + flat_eofs = self._flatE[slicer].copy() + if eofscaling == 0: + # No modification. A copy needs to be returned in case it is + # modified. If no copy is made the internally stored eigenvectors + # could be modified unintentionally. + #rval = self._flatE[slicer].copy() + rval = flat_eofs + elif eofscaling == 1: + # Divide by the square-root of the eigenvalues. + rval = flat_eofs / np.sqrt(self._L[slicer])[:, np.newaxis] + elif eofscaling == 2: + # Multiply by the square-root of the eigenvalues. + rval = flat_eofs * np.sqrt(self._L[slicer])[:, np.newaxis] + else: + raise ValueError('invalid eof scaling option: ' + '{!s}'.format(eofscaling)) + if self._filled: + rval = ma.array(rval, mask=np.where(np.isnan(rval), True, False)) + return rval.reshape(((neeofs, self._window,) + self._originalshape)) + + def eigenvalues(self, neigs=None): + """Eigenvalues (decreasing variances) associated with each EEOF. + + **Optional argument:** + + *neigs* + Number of eigenvalues to return. Defaults to all + eigenvalues. If the number of eigenvalues requested is more + than the number that are available, then all available + eigenvalues will be returned. + + **Returns:** + + *eigenvalues* + An array containing the eigenvalues arranged largest to + smallest. + + **Examples:** + + All eigenvalues:: + + eigenvalues = solver.eigenvalues() + + The first eigenvalue:: + + eigenvalue1 = solver.eigenvalues(neigs=1) + + """ + # Create a slicer and use it on the eigenvalue array. A copy must be + # returned in case the slicer takes all elements, in which case a + # reference to the eigenvalue array is returned. If this is modified + # then the internal eigenvalues array would then be modified as well. + if neigs is None or neigs > self.neeofs: + neigs = self.neeofs + slicer = slice(0, neigs) + return self._L[slicer].copy() + + def varianceFraction(self, neigs=None): + """Fractional EEOF mode variances. + + The fraction of the total variance explained by each EOF mode, + values between 0 and 1 inclusive. + + **Optional argument:** + + *neigs* + Number of eigenvalues to return the fractional variance for. + Defaults to all eigenvalues. If the number of eigenvalues + requested is more than the number that are available, then + fractional variances for all available eigenvalues will be + returned. + + **Returns:** + + *variance_fractions* + An array containing the fractional variances. + + **Examples:** + + The fractional variance represented by each EEOF mode:: + + variance_fractions = solver.varianceFraction() + + The fractional variance represented by the first EEOF mode:: + + variance_fraction_mode_1 = solver.VarianceFraction(neigs=1) + + """ + # Return the array of eigenvalues divided by the sum of the + # eigenvalues. + if neigs is None or neigs > self.neeofs: + neigs = self.neeofs + slicer = slice(0, neigs) + return sum(self._L[slicer]) / self._L.sum() + +