<commit_before># coding=utf-8
"""
Main Linear Inverse Model classes and methods.

Author: <NAME>
"""

import numpy as np
from numpy.linalg import pinv, cond, eigvals
from math import ceil
import tables as tb
import cPickle as cpk
import logging
import time

from Stats import calc_eofs
import DataTools as Dt


logger = logging.getLogger(__name__)


def _calc_m(x0, xt, tau=1):
    """Calculate either L or G for forecasting (using nomenclature
    from Newman 2013"""
    
    # These represent the C(tau) and C(0) covariance matrices
    #    Note: x is an anomaly vector, no division by N-1 because it's undone
    #    in the inversion anyways
    
    x0x0 = np.dot(x0, x0.T)
    x0xt = np.dot(xt, x0.T)

    # Calculate the mapping term G_tau
    G = np.dot(x0xt, pinv(x0x0))

    # Calculate the forcing matrix to check that all modes are damped
    L = (1/tau) * np.log(G)
    Leigs = eigvals(L)

    if np.any(Leigs.real >= 0):
        logger.debug('L eigenvalues: \n' + str(Leigs))
        raise ValueError('Positive eigenvalues detected in forecast matrix L.')

    if cond(x0x0) > 20:
        logger.warn(('C_0 condition number is large ({:2.2f}). Too many '
                     'features (or EOFs) are likely provided and further '
                     'truncation is suggested.'))

    return G


def _create_h5_fcast_grps(h5f, dest, atom, shape, fcast_times):
    """
    Helper method for creating forecast groups in the hdf5 file.
    """

    out_fcast = []

    # if not h5f.__contains__(dest):
    #     grp, node = path.split(dest)
    #     h5f.create_group(grp, node)

    for lead in fcast_times:
        out_fcast.append(
            Dt.empty_hdf5_carray(h5f,
                                 dest,
                                 'f{:d}'.format(lead),
                                 atom,
                                 shape,
                                 createparents=True,
                                 title='{:d} Year Forecast'.format(lead)))

    return out_fcast


class LIM(object):
    """Linear inverse forecast model.
    
    This class uses a calibration dataset to make simple linear forecasts.
    
    Notes
    -----
    Based on the LIM described by <NAME> (2013) [1].  Right now it
    assumes the use of monthly data (i.e. each timestep should represent a
    single month).
    
    References
    ----------
    .. [1] <NAME>. (2013), An Empirical Benchmark for Decadal Forecasts of 
       Global Surface Temperature Anomalies, <NAME>., 26(14), 5260–5269, 
       doi:10.1175/JCLI-D-12-00590.1.
    ....
    """

    def __init__(self, calib_data, nelem_in_tau1, h5file=None):
        """
        Parameters
        ----------
        calib_data: ndarray
            Data for calibrating the LIM.  Expects
            a 2D MxN matrix where M (rows) represent the sampling dimension and
            N(columns) represents the feature dimension (e.g. spatial grid
            points).
        nelem_in_tau1: int
            Number of time samples that span tau=1.  E.g. for monthly data when
            a forecast tau is equivalent to 1 year, nelem_in_tau should be 12.
        h5file: HDF5_Object, Optional
            File object to store LIM output.  It will create a series of
            directories under the given group
        """
        logger.info('Initializing LIM forecasting object...')

        if calib_data.ndim != 2:
            logger.error(('LIM calibration data is not 2D '
                          '(Contained ndim={:d}').format(calib_data.ndim))
            raise ValueError('Input LIM calibration data is not 2D')

        self._h5file = h5file
        self._from_precalib = False
        self._nelem_in_tau1 = nelem_in_tau1
        self._eof_var_stats = {}

        x0 = calib_data[:, 0:-nelem_in_tau1]
        x1 = calib_data[:, nelem_in_tau1:]<commit_msg>calc_m fix<commit_after># coding=utf-8
"""
Main Linear Inverse Model classes and methods.

Author: <NAME>
"""

import numpy as np
from numpy.linalg import pinv, cond, eigvals
from math import ceil
import tables as tb
import cPickle as cpk
import logging
import time

from Stats import calc_eofs
import DataTools as Dt


logger = logging.getLogger(__name__)


def _calc_m(x0, xt, tau=1):
    """Calculate either L or G for forecasting (using nomenclature
    from Newman 2013
    
    Parameters
    ----------
    x0: ndarray
        State at time=0.  MxN where M is number of samples, and N is the 
        number of features.
    xt: ndarray
        State at time=tau.  MxN where M is number of samples, and N is the 
        number of fatures.
    tau: float
        lag time (in units of tau) that we are calculating G for.  This is 
        used to check that all modes of L are damped.
        
        
    """
    
    # These represent the C(tau) and C(0) covariance matrices
    #    Note: x is an anomaly vector, no division by N-1 because it's undone
    #    in the inversion anyways
    
    x0x0 = np.dot(x0.T, x0)
    x0xt = np.dot(xt.T, x0)

    # Calculate the mapping term G_tau
    G = np.dot(x0xt, pinv(x0x0))

    # Calculate the forcing matrix to check that all modes are damped
    L = (1./tau) * np.log(G)
    Leigs = eigvals(L)

    if np.any(Leigs.real >= 0):
        logger.debug('L eigenvalues: \n' + str(Leigs))
        raise ValueError('Positive eigenvalues detected in forecast matrix L.')

    if cond(x0x0) > 20:
        logger.warn(('C_0 condition number is large ({:2.2f}). Too many '
                     'features (or EOFs) are likely provided and further '
                     'truncation is suggested.'))

    return G


def _create_h5_fcast_grps(h5f, dest, atom, shape, fcast_times):
    """
    Helper method for creating forecast groups in the hdf5 file.
    """

    out_fcast = []

    # if not h5f.__contains__(dest):
    #     grp, node = path.split(dest)
    #     h5f.create_group(grp, node)

    for lead in fcast_times:
        out_fcast.append(
            Dt.empty_hdf5_carray(h5f,
                                 dest,
                                 'f{:d}'.format(lead),
                                 atom,
                                 shape,
                                 createparents=True,
                                 title='{:d} Year Forecast'.format(lead)))

    return out_fcast


class LIM(object):
    """Linear inverse forecast model.
    
    This class uses a calibration dataset to make simple linear forecasts.
    
    Notes
    -----
    Based on the LIM described by <NAME> (2013) [1].  Right now it
    assumes the use of monthly data (i.e. each timestep should represent a
    single month).
    
    References
    ----------
    .. [1] <NAME>. (2013), An Empirical Benchmark for Decadal Forecasts of 
       Global Surface Temperature Anomalies, J. Clim., 26(14), 5260–5269, 
       doi:10.1175/JCLI-D-12-00590.1.
    ....
    """

    def __init__(self, calib_data, nelem_in_tau1, h5file=None):
        """
        Parameters
        ----------
        calib_data: ndarray
            Data for calibrating the LIM.  Expects
            a 2D MxN matrix where M (rows) represent the sampling dimension and
            N(columns) represents the feature dimension (e.g. spatial grid
            points).
        nelem_in_tau1: int
            Number of time samples that span tau=1.  E.g. for monthly data when
            a forecast tau is equivalent to 1 year, nelem_in_tau should be 12.
        h5file: HDF5_Object, Optional
            File object to store LIM output.  It will create a series of
            directories under the given group
        """
        logger.info('Initializing LIM forecasting object...')

        if calib_data.ndim != 2:
            logger.error(('LIM calibration data is not 2D '
                          '(Contained ndim={:d}').format(calib_data.ndim))
            raise ValueError('Input LIM calibration data is not 2D')

        self._h5file = h5file
        self._from_precalib = False
        self._nelem_in_tau1 = nelem_in_tau1
        self._eof_var_stats = {}

        x0 = calib_data[0:-nelem_in_tau1, :]
        x1 = calib_data[nelem_in_tau1:, :]
{'fraction_non_alphanumeric': 0.09193846153846154, 'fraction_numerical': 0.02289230769230769, 'mean_word_length': 1.9156799425905993, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1245506', 'n_tokens_mistral': 2644, 'n_tokens_neox': 2440, 'n_words': 820}