Coverage for calorine/tools/analysis.py: 100%

1import numpy as np 

2import pandas as pd 

3from ase import units 

4from scipy import stats 

5 

6 

7def analyze_data(data: np.ndarray, max_lag: int = None) -> dict: 

8 """ Carries out an extensive analysis of the data series. 

9 

10 Parameters 

11 ---------- 

12 data 

13 data series to compute autocorrelation function for 

14 max_lag 

15 maximum lag between two data points, used for computing autocorrelation 

16 

17 Returns 

18 ------- 

19 dict 

20 calculated properties of the data including, mean, standard deviation, 

21 correlation length and a 95% error estimate. 

22 """ 

23 summary = dict(mean=data.mean(), 

24 std=data.std()) 

25 acf = get_autocorrelation_function(data, max_lag) 

26 correlation_length = _estimate_correlation_length_from_acf(acf) 

27 if correlation_length is not None: 

28 error_estimate = _estimate_error(data, correlation_length, confidence=0.95) 

29 summary['correlation_length'] = correlation_length 

30 summary['error_estimate'] = error_estimate 

31 else: 

32 summary['correlation_length'] = np.nan 

33 summary['error_estimate'] = np.nan 

34 return summary 

35 

36 

37def get_autocorrelation_function(data: np.ndarray, max_lag: int = None) -> np.ndarray: 

38 """ Returns autocorrelation function. 

39 

40 The autocorrelation function is computed using `pandas.Series.autocorr 

41 <https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Series.autocorr.html>`_. 

42 

43 Parameters 

44 ---------- 

45 data 

46 data series to compute autocorrelation function for 

47 max_lag 

48 maximum lag between two data points 

49 

50 Returns 

51 ------- 

52 calculated autocorrelation function 

53 """ 

54 if max_lag is None: 

55 max_lag = len(data) - 1 

56 if max_lag < 1 or max_lag >= len(data): 

57 raise ValueError('max_lag should be between 1 and len(data)-1.') 

58 series = pd.Series(data) 

59 acf = [series.autocorr(lag) for lag in range(0, max_lag)] 

60 return np.array(acf) 

61 

62 

63def get_correlation_length(data: np.ndarray) -> int | None: 

64 """ Returns estimate of the correlation length of data. 

65 

66 The correlation length is taken as the first point where the 

67 autocorrelation functions is less than :math:`\\exp(-2)`. If the 

68 correlation function never drops below :math:`\\exp(-2)` ``np.nan`` is 

69 returned. 

70 

71 If the correlation length cannot be computed since the auto-correlation 

72 function is unconverged the function returns `None`. 

73 

74 Parameters 

75 ---------- 

76 data 

77 data series for which to the compute autocorrelation function 

78 

79 Returns 

80 ------- 

81 correlation length 

82 """ 

83 

84 acf = get_autocorrelation_function(data) 

85 correlation_length = _estimate_correlation_length_from_acf(acf) 

86 if correlation_length is None: 

87 return None 

88 return correlation_length 

89 

90 

91def get_error_estimate(data: np.ndarray, confidence: float = 0.95) -> float | None: 

92 """Returns estimate of standard error :math:`\\mathrm{error}` 

93 with confidence interval. 

94 

95 .. math:: 

96 

97 \\mathrm{error} = t_\\mathrm{factor} * \\mathrm{std}(\\mathrm{data}) / \\sqrt{N_s} 

98 

99 where :math:`t_{factor}` is the factor corresponding to the confidence 

100 interval and :math:`N_s` is the number of independent measurements 

101 (with correlation taken into account). 

102 

103 If the correlation length cannot be computed since the auto-correlation 

104 function is unconverged the function returns `None`. 

105 

106 Parameters 

107 ---------- 

108 data 

109 data series for which to estimate the error 

110 confidence 

111 confidence level of the interval 

112 

113 Returns 

114 ------- 

115 error estimate 

116 

117 """ 

118 correlation_length = get_correlation_length(data) 

119 if correlation_length is None: 

120 return None 

121 error_estimate = _estimate_error(data, correlation_length, confidence) 

122 return error_estimate 

123 

124 

125def _estimate_correlation_length_from_acf(acf: np.ndarray) -> int | None: 

126 """Estimates correlation length from acf. Returns None if the auto-correlation 

127 function is uncoverged.""" 

128 for i, a in enumerate(acf): 

129 if a < np.exp(-2): 

130 return i 

131 return None # np.nan 

132 

133 

134def _estimate_error(data: np.ndarray, 

135 correlation_length: int, 

136 confidence: float) -> float: 

137 """ Estimates error using correlation length. """ 

138 t_factor = stats.t.ppf((1 + confidence) / 2, len(data)-1) # type: float 

139 error = t_factor * np.std(data) / np.sqrt(len(data) / correlation_length) # type: float 

140 return error 

141 

142 

143def get_rtc_from_hac(hac: np.ndarray, V: float, T: float, dt: float) -> np.ndarray: 

144 """Returns the running thermal conductivity (RTC) in W/m/K using a 

145 heat auto-correlation (HAC) as input. 

146 

147 Parameters 

148 ---------- 

149 hac 

150 The HAC :math:`\\langle j(t)j(0)\\rangle` in units of eV\\ :sup:`3`/amu as given by GPUMD 

151 V 

152 Volume of cell in Å\\ :sup:`3` 

153 T 

154 Temperature in Kelvin 

155 dt 

156 Time step in the HAC in ps 

157 running thermal conductivity in W/m/K 

158 """ 

159 

160 hac = np.asarray(hac) 

161 

162 if hac.ndim > 1: 

163 raise ValueError('hac must be 1D') 

164 

165 kB = units.kB # eV/K 

166 

167 J_per_eV = 1 / units.J 

168 kg_per_u = 1 / units.kg 

169 s_per_ps = 1e-12 

170 m_per_Å = 1 / units.m 

171 

172 # The integration and normalization 

173 rtc = np.cumsum(hac) * dt / (kB * T**2 * V) 

174 

175 # From eV/ps/Å/K to J/s/m/K 

176 rtc *= (J_per_eV**3 / kg_per_u) * s_per_ps / (J_per_eV * m_per_Å**3) 

177 

178 return rtc