import os
import numpy as np
import pandas as pd
from logging import info
from pycif.utils.path import init_dir
from .....utils.check.errclass import CifError, CifIOError
[docs]
def build_tcorrelations(
period,
subperiod,
dates,
sigma_t,
sigma_type,
evalmin=0.5,
dump=False,
dir_dump="",
crop_chi=False,
corr_plg=None,
tracer=None,
**kwargs
):
"""Build temporal correlation matrix based on timedelta between periods.
For period i and j, the corresponding correlation is:
c(i,j) = exp(-timedelta(i, j) / sigma)
Args:
period (int): period duration
subperiod (int): sub-period duration
dates (np.array): dates sub-dividing the control vector periods
sigma_t (float): decay distance for correlation between periods
(in days)
evalmin (float): flag out all eigenvalues below this value.
Default is 0.5
dump (bool): dumps computed correlations if True
dir_dump (str): directory where correlation matrices are stored
Return:
tuple with:
- square roots of eigenvalues
- eigenvectors
"""
# Check ambiguous units
if sigma_t[-1].lower() == "m":
raise CifError(
"WARNING! Your `sigma_t` unit end by 'm' which is an ambiguous unit. \n"
"if you mean a temporal correlation with monthly unit, please convert it "
"into days unit. If you want minute units, please state the full `min` unit."
)
# Try reading existing file
try:
evalues, evectors = read_tcorr(period, subperiod, dates,
sigma_t, sigma_type, dir_dump)
# Else build correlations
except IOError:
info("Computing temporal correlations")
temp_distance = np.abs(
pd.DatetimeIndex(dates).values[:, np.newaxis]
- pd.DatetimeIndex(dates).values[np.newaxis, :]
)
corr = np.eye(dates.size)
if sigma_type == "isotrope":
# Compute matrix of distance
dt = temp_distance / pd.to_timedelta(sigma_t)
# Compute the correlation matrix itself
corr = np.exp(-dt)
elif sigma_type == "frequency":
freq = pd.to_timedelta(corr_plg.freq).to_numpy()
mask = np.mod(temp_distance, freq) == np.timedelta64(0)
corr[mask] = np.exp(
-np.abs(temp_distance[mask] / pd.to_timedelta(sigma_t)))
elif sigma_type == "category":
scale = corr_plg.scale
if scale == "hourofday":
category = pd.DatetimeIndex(dates).hour.values
ntimes = 24
dist_scale = np.timedelta64(1, "h")
elif scale == "dayofweek":
category = pd.DatetimeIndex(dates).dayofweek.values
ntimes = 7
dist_scale = np.timedelta64(1, "D")
elif scale == "monthofyear":
category = pd.DatetimeIndex(dates).month.values
ntimes = 12
dist_scale = np.timedelta64(30, "D")
else:
raise CifError(
f"Scale for 'category' temporal correlations unknown:{scale}")
# Compute and loop distance
dist = category[:, np.newaxis] - category[np.newaxis, :]
dist = (dist + ntimes / 2) % ntimes - ntimes / 2
dist = dist.astype(int) * dist_scale
corr = np.exp(-np.abs(dist / pd.to_timedelta(sigma_t)))
else:
raise CifError(
f"Type for temporal correlations is not recognized:{sigma_type}")
# Component analysis
evalues, evectors = np.linalg.eigh(corr)
# Re-ordering values
# (not necessary in principle in recent numpy versions)
index = np.argsort(evalues)[::-1]
evalues = evalues[index]
evectors = evectors[:, index]
# Dumping to a txt file
if dump:
dump_tcorr(period, subperiod, dates, sigma_t, sigma_type,
evalues, evectors,
f"{tracer.workdir}/controlvect/correlations/")
except Exception as e:
raise e
# Truncating values < evalmin
mask = evalues >= evalmin
if crop_chi:
return evalues[mask] ** 0.5, evectors[:, mask]
else:
evalues[~mask] = 0
return evalues ** 0.5, evectors
[docs]
def dump_tcorr(period, subperiod, dates, sigma_t, sigma_type, evalues, evectors,
dir_dump, overwrite=False):
"""Dumps eigenvalues and vectors to a bin file. The default file format
is:
f"{dir_dump}/tempcor_{datei.strftime('%Y%m%d%H%M')}_{datef.strftime('%Y%m%d%H%M')}_per{period}-{subperiod}_ct{sigma_t}_{sigma_type}.bin"
Args:
period (int): period duration
subperiod (int): subperiod duration
dates (np.array): dates sub-dividing the control vector periods
sigma_t (float): decay distance for correlation between periods
(in days)
"""
datei = dates[0]
datef = dates[-1]
file_dump = f"{dir_dump}/tempcor_{datei.strftime('%Y%m%d%H%M')}_{datef.strftime('%Y%m%d%H%M')}_per{period}-{subperiod}_ct{sigma_t}_{sigma_type}.bin"
if os.path.isfile(file_dump) and overwrite:
raise CifIOError(
f"Warning: {file_dump} already exists. I don't want to overwrite it"
)
# Initiate path if needed
if not os.path.isdir(dir_dump):
init_dir(dir_dump)
datasave = np.concatenate((evalues[np.newaxis, :], evectors), axis=0)
datasave.tofile(file_dump)
[docs]
def read_tcorr(period, subperiod, dates, sigma_t, sigma_type, dir_dump):
"""Reads temporal correlations from existing bin file
Args:
period (int): period duration
subperiod (int): subperiod duration
dates (np.array): dates sub-dividing the control vector periods
sigma_t (float): decay distance for correlation between periods
(in days)
dir_dump (str): where the horizontal correlations have been stored
"""
datei = dates[0]
datef = dates[-1]
file_dump = f"{dir_dump}/tempcor_{datei.strftime('%Y%m%d%H%M')}_{datef.strftime('%Y%m%d%H%M')}_per{period}-{subperiod}_ct{sigma_t}_{sigma_type}.bin"
if not os.path.isfile(file_dump):
raise CifIOError(
f"{file_dump} does not exist. Please compute correlations from scratch"
)
data = np.fromfile(file_dump).reshape((-1, len(dates)))
evalues = data[0]
evectors = data[1:]
evalues[evalues < 0] = 0.0
return evalues, evectors
