From d6238fd337ffa019e1bc5d843830e5e589ba0c0c Mon Sep 17 00:00:00 2001 From: Colladito Date: Tue, 23 Sep 2025 15:04:24 +0200 Subject: [PATCH 01/18] [VCF] Changing axis limit in pca plot --- bluemath_tk/teslakit/plotting/pcs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bluemath_tk/teslakit/plotting/pcs.py b/bluemath_tk/teslakit/plotting/pcs.py index d0288fc..67eaaa5 100755 --- a/bluemath_tk/teslakit/plotting/pcs.py +++ b/bluemath_tk/teslakit/plotting/pcs.py @@ -99,7 +99,7 @@ def axplot_PCs_3D_WTs(ax, d_PCs, wt_colors, ttl='PCs'): ax.set_ylabel('PC2', {'fontsize':10}) ax.set_zlabel('PC3', {'fontsize':10}) - lim = np.nanmax([np.nanmax(PC1), np.nanmax(PC2), np.nanmax(PC3)]) + lim = np.nanmax([np.nanmax(np.abs(PC1)), np.nanmax(np.abs(PC2)), np.nanmax(np.abs(PC3))]) # ax.set_xlim([-3,3]) # ax.set_ylim([-3,3]) From a3ebaf0fb05692be6334673895265694cc2924ae Mon Sep 17 00:00:00 2001 From: Colladito Date: Thu, 2 Oct 2025 20:03:22 +0200 Subject: [PATCH 02/18] [VCF] Update NonStatGEV plots --- bluemath_tk/distributions/nonstat_gev.py | 435 ++++++++++++----------- 1 file changed, 226 insertions(+), 209 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index 7be7fe1..912dc1b 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -4797,7 +4797,7 @@ def _quantile(self, prob=None, harm=False) -> np.ndarray: return Q - def plot(self, return_plot: bool = True, save: bool = False): + def plot(self, return_plot: bool = True, save: bool = False, init_year=2000): """ Plot the location, scale and shape parameters, also the PP plot and QQ plot @@ -4946,6 +4946,8 @@ def plot(self, return_plot: bool = True, save: bool = False): # Location and Scale parameter plotting t_anual = np.mod(self.t, 1) quan95 = self._quantile() + rp10 = self._quantile(1-1/10) + rp100 = self._quantile(1-1/100) if ( self.betaT is None @@ -4965,7 +4967,7 @@ def plot(self, return_plot: bool = True, save: bool = False): label="Data", ) ax2 = ax1.twinx() - l1 = ax1.plot( + ax1.plot( t_anual[t_ord], mut[t_ord], label="Location", @@ -4980,7 +4982,7 @@ def plot(self, return_plot: bool = True, save: bool = False): color=self.colors[0], alpha=0.3, ) - l2 = ax2.plot( + ax2.plot( t_anual[t_ord], psit[t_ord], label="Scale", @@ -4995,7 +4997,7 @@ def plot(self, return_plot: bool = True, save: bool = False): color=self.colors[1], alpha=0.3, ) - l3 = ax1.plot( + ax1.plot( t_anual[t_ord], quan95[t_ord], linestyle="dashed", @@ -5004,8 +5006,8 @@ def plot(self, return_plot: bool = True, save: bool = False): ) else: fig, ax1 = plt.subplots(figsize=(20, 6)) - l0 = ax1.plot( - self.t, + ax1.plot( + self.t + init_year, self.xt, marker="+", linestyle="None", @@ -5013,9 +5015,8 @@ def plot(self, return_plot: bool = True, save: bool = False): markersize=5, label="Data", ) - ax2 = ax1.twinx() - l1 = ax1.plot( - self.t, + ax1.plot( + self.t+init_year, mut, label="Location", linewidth=2, @@ -5023,31 +5024,44 @@ def plot(self, return_plot: bool = True, save: bool = False): color="tab:blue", alpha=1, ) - # ax1.fill_between( - # self.t, ci_low_mut, ci_up_mut, color=self.colors[0], alpha=0.3 - # ) + ax1.fill_between( + self.t+init_year, + mut - psit, + mut + psit, + label=r"Location $\pm$ Scale", + color="tab:blue", + alpha=0.4 + ) - l2 = ax2.plot( - self.t, - psit, - label="Scale", - linewidth=2, - # color=self.colors[1], + ax1.plot( + self.t+init_year, + rp10, + label="10 years", + linewidth=1, + color="tab:red", + linestyle="--", + alpha=.9, + ) + + ax1.plot( + self.t+init_year, + rp100, + label="100 years", + linewidth=1, color="tab:green", - alpha=1, + linestyle="--", + alpha=.9, ) - # ax2.fill_between( - # self.t, ci_low_psit, ci_up_psit, color=self.colors[1], alpha=0.3 - # ) - l3 = ax1.plot( - self.t, - quan95, - linestyle="dashed", - # color=self.colors[2], - color="tab:orange", + + ax1.plot( + self.t+init_year, + epst, + label="Shape", linewidth=1, - label="Quantile 95%", + color="tab:orange", + alpha=.9, ) + # TODO: Add aggregated return period lines # rt_10 = np.zeros(40) # for year in range(40): @@ -5075,58 +5089,21 @@ def plot(self, return_plot: bool = True, save: bool = False): # ) ax1.set_xlabel("Time (years)") - ax1.set_ylabel(r"$\mu_t(m), Hs(m)$") - ax2.set_ylabel(r"$\psi_t(m)$") + ax1.set_ylabel(rf"$\mu_t(m)$, {self.var_name}") # ax1.set_title(f"Evolution of location and scale parameters ({self.var_name})") - ax1.set_title(f"Evolution of adjustment ({self.var_name})") + ax1.set_title(f"Evolution of parameters ({self.var_name})") ax1.grid(True) - handles = [ - # art for art in l0 + l1 + l2 + l3 + l4 + l5 + l6 + l7 if not art.get_label().startswith("_") - art - for art in l0 + l1 + l2 + l3 - if not art.get_label().startswith("_") - ] - ax1.legend(handles=handles, loc="best") + ax1.legend(loc="best") # ax2.set_ylim(0,1.5) # ax1.set_xlim(-0.07, 40) ax1.margins(x=0.01) + plt.tight_layout() if save: os.makedirs("Figures", exist_ok=True) # plt.savefig(f"Figures/Location_Scale_Parameters_{self.var_name}.png", dpi=300) plt.savefig(f"Figures/Adjustment_Evolution_{self.var_name}.png", dpi=300) plt.show() - # mu, mu-phi, mu+phi, points - fig, ax1 = plt.subplots(figsize=(10, 6)) - l0 = ax1.plot( - self.t, - self.xt, - marker="+", - linestyle="None", - color="black", - markersize=5, - label=r"$H_s^{max}$", - ) - l1 = ax1.plot( - self.t, mut, label=r"$\mu_t$", linewidth=2, color=self.colors[0], alpha=1 - ) - ax1.fill_between( - self.t, mut - psit, mut + psit, color=self.colors[1], alpha=0.3 - ) - l3 = ax1.plot( - self.t, quan95, linestyle="dashed", color=self.colors[2], linewidth=2 - ) - ax1.set_xlabel("Time (yearly scale)") - ax1.set_ylabel(r"$\mu_t$") - ax1.set_title(f"Location parameter ({self.var_name})") - ax1.grid(True) - handles = [art for art in l0 + l1 + l3 if not art.get_label().startswith("_")] - ax1.legend(handles=handles, loc="best") - ax1.margins(x=0.01) - if save: - plt.savefig(f"Figures/Location_Parameter_{self.var_name}.png", dpi=300) - plt.show() - month_initials = ["J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"] month_positions = [(i + 0.5) / 12 for i in range(12)] @@ -5140,7 +5117,7 @@ def plot(self, return_plot: bool = True, save: bool = False): linestyle="None", color="black", markersize=5, - label=r"$H_s^{max}$", + label=f"{self.var_name}", ) # ax2 = ax1.twinx() l1 = ax1.plot( @@ -5198,7 +5175,7 @@ def plot(self, return_plot: bool = True, save: bool = False): linestyle="None", color="black", markersize=5, - label=r"$H_s^{max}$", + label=f"{self.var_name}", ) # ax2 = ax1.twinx() l1 = ax1.plot( @@ -5221,6 +5198,7 @@ def plot(self, return_plot: bool = True, save: bool = False): uppermaxs, color=self.colors[0], alpha=0.3, + label="Location +- scale" ) # l2 = ax2.plot(self.t[mask_month], psit[mask_month], label=r'$\psi_t$', linewidth=2, color=self.colors[1], alpha=1) # ax2.fill_between(self.t[mask_month], ci_low_psit[mask_month], ci_up_psit[mask_month], color=self.colors[1], @@ -5246,148 +5224,187 @@ def plot(self, return_plot: bool = True, save: bool = False): ax1.legend(handles=handles, loc="best") ax1.margins(x=0.01) if save: - plt.savefig(f"Figures/Shape_Parameter_{self.var_name}.png", dpi=300) + plt.savefig(f"Figures/Evolution_Location_FirstYear{self.var_name}.png", dpi=300) plt.show() ### Shape parameter plot - if self.ngamma > 0: - t_ord = np.argsort(t_anual) - - # Confidence interval for epst - ci_up = ( - epst + norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdepst - ) - ci_low = ( - epst - norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdepst - ) - - plt.figure(figsize=(10, 6)) - plt.plot(t_anual[t_ord], epst[t_ord], color=self.colors[0]) - plt.fill_between( - t_anual[t_ord], - ci_low[t_ord], - ci_up[t_ord], - color=self.colors[0], - alpha=0.3, - label=r"$\xi_t$ Confidence Interval", - ) - plt.title(f"Shape parameter ({self.var_name})") - plt.xlabel("Time (yearly scale)") - plt.ylabel(r"$\xi_t$") - plt.xticks(month_positions, month_initials) - plt.grid(True) - plt.show() - - ### Harmonic Location parameter plot - if self.nmu > 0: - t_ord = np.argsort(t_anual) - quan95_2 = self._quantile(harm=True) - - mut2 = self._parametro(self.beta0, self.beta) - # Confidence interval for mut - ci_up = mut2 + norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdmut - ci_low = ( - mut2 - norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdmut - ) - - plt.figure(figsize=(10, 6)) - plt.plot( - t_anual[t_ord], - self.xt[t_ord], - marker="+", - linestyle="None", - color="black", - markersize=5, - label=r"$H_s^{max}$", - ) - plt.plot( - t_anual[t_ord], - mut2[t_ord], - label=r"$\mu_t$", - linewidth=2, - color=self.colors[0], - ) - plt.fill_between( - t_anual[t_ord], - ci_low[t_ord], - ci_up[t_ord], - color=self.colors[0], - alpha=0.3, - ) - # Confidence interval for the quantile - ci_up = ( - quan95_2 + norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdDq - ) - ci_low = ( - quan95_2 - norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdDq - ) - plt.plot( - t_anual[t_ord], - quan95_2[t_ord], - linestyle="dashed", - color=self.colors[1], - markersize=5, - label=rf"$q_{self.quanval}$", - ) - plt.fill_between( - t_anual[t_ord], - ci_low[t_ord], - ci_up[t_ord], - color=self.colors[1], - alpha=0.3, - ) - plt.title(f"Harmonic part of Location parameter ({self.var_name})") - plt.xlabel("Time (yearly scale)") - plt.ylabel(r"$\mu_t$") - plt.xticks(month_positions, month_initials) - plt.legend(loc="best") - plt.grid(True) - if save: - plt.savefig( - f"Figures/Harmonic_Location_Parameter_{self.var_name}.png", dpi=300 - ) - plt.show() - - ### Scale parameter plot - if self.npsi > 0: - t_ord = np.argsort(t_anual) + # Confidence interval for epst + # ci_up = ( + # epst + norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdepst + # ) + # ci_low = ( + # epst - norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdepst + # ) + + # LOCATION + plt.figure(figsize=(20, 6)) + plt.plot(self.t + init_year, mut, color="tab:blue") + # plt.fill_between( + # t_anual[t_ord], + # ci_low[t_ord], + # ci_up[t_ord], + # color=self.colors[0], + # alpha=0.3, + # label=r"$\xi_t$ Confidence Interval", + # ) + plt.title(f"Location parameter ({self.var_name})") + plt.xlabel("Time (yearly scale)") + plt.ylabel(r"$\mu_t$") + plt.grid(True) + if save: + plt.savefig(f"Figures/Evolution_Location{self.var_name}.png", dpi=300) + plt.show() - psit2 = np.exp(self._parametro(self.alpha0, self.alpha)) - # Confidence interval for psit - ci_up = ( - psit2 + norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdpsit - ) - ci_low = ( - psit2 - norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdpsit - ) + # SCALE + plt.figure(figsize=(20, 6)) + plt.plot(self.t + init_year, psit, color="tab:green") + # plt.fill_between( + # t_anual[t_ord], + # ci_low[t_ord], + # ci_up[t_ord], + # color=self.colors[0], + # alpha=0.3, + # label=r"$\xi_t$ Confidence Interval", + # ) + plt.title(f"Scale parameter ({self.var_name})") + plt.xlabel("Time (yearly scale)") + plt.ylabel(r"$\psi_t$") + plt.grid(True) + if save: + plt.savefig(f"Figures/Evolution_Scale{self.var_name}.png", dpi=300) + plt.show() - plt.figure(figsize=(10, 6)) - plt.plot( - t_anual[t_ord], - psit2[t_ord], - label=r"$\psi_t$", - linewidth=2, - color=self.colors[0], - ) - plt.fill_between( - t_anual[t_ord], - ci_low[t_ord], - ci_up[t_ord], - color=self.colors[0], - alpha=0.3, - label=r"$\psi_t$ Confidence Interval", - ) - # plt.plot(t_anual[t_ord], quan95[t_ord], linestyle='dashed', color=self.colors[2], markersize=5, label=fr"$q_{self.quanval}$") - plt.title(f"Harmonic part of Scale parameter ({self.var_name})") - plt.xlabel("Time (yearly scale)") - plt.xticks(month_positions, month_initials) - plt.ylabel(r"$\psi_t$") - plt.grid(True) - if save: - plt.savefig( - f"Figures/Harmonic_Scale_Parameter_{self.var_name}.png", dpi=300 - ) - plt.show() + # SHAPE + plt.figure(figsize=(20, 6)) + plt.plot(self.t + init_year, epst, color="tab:orange") + # plt.fill_between( + # t_anual[t_ord], + # ci_low[t_ord], + # ci_up[t_ord], + # color=self.colors[0], + # alpha=0.3, + # label=r"$\xi_t$ Confidence Interval", + # ) + plt.title(f"Shape parameter ({self.var_name})") + plt.xlabel("Time (yearly scale)") + plt.ylabel(r"$\xi_t$") + plt.grid(True) + if save: + plt.savefig(f"Figures/Evolution_Shape{self.var_name}.png", dpi=300) + plt.show() + + + + # ### Harmonic Location parameter plot + # if self.nmu > 0: + # t_ord = np.argsort(t_anual) + # quan95_2 = self._quantile(harm=True) + + # mut2 = self._parametro(self.beta0, self.beta) + # # Confidence interval for mut + # ci_up = mut2 + norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdmut + # ci_low = ( + # mut2 - norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdmut + # ) + + # plt.figure(figsize=(10, 6)) + # plt.plot( + # t_anual[t_ord], + # self.xt[t_ord], + # marker="+", + # linestyle="None", + # color="black", + # markersize=5, + # label=f"{self.var_name}", + # ) + # plt.plot( + # t_anual[t_ord], + # mut2[t_ord], + # label=r"$\mu_t$", + # linewidth=2, + # color=self.colors[0], + # ) + # plt.fill_between( + # t_anual[t_ord], + # ci_low[t_ord], + # ci_up[t_ord], + # color=self.colors[0], + # alpha=0.3, + # ) + # # Confidence interval for the quantile + # ci_up = ( + # quan95_2 + norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdDq + # ) + # ci_low = ( + # quan95_2 - norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdDq + # ) + # plt.plot( + # t_anual[t_ord], + # quan95_2[t_ord], + # linestyle="dashed", + # color=self.colors[1], + # markersize=5, + # label=rf"$q_{self.quanval}$", + # ) + # plt.fill_between( + # t_anual[t_ord], + # ci_low[t_ord], + # ci_up[t_ord], + # color=self.colors[1], + # alpha=0.3, + # ) + # plt.title(f"Harmonic part of Location parameter ({self.var_name})") + # plt.xlabel("Time (yearly scale)") + # plt.ylabel(r"$\mu_t$") + # plt.xticks(month_positions, month_initials) + # plt.legend(loc="best") + # plt.grid(True) + # if save: + # plt.savefig( + # f"Figures/Harmonic_Location_Parameter_{self.var_name}.png", dpi=300 + # ) + # plt.show() + + # ### Scale parameter plot + # if self.npsi > 0: + # t_ord = np.argsort(t_anual) + + # psit2 = np.exp(self._parametro(self.alpha0, self.alpha)) + # # Confidence interval for psit + # ci_up = ( + # psit2 + norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdpsit + # ) + # ci_low = ( + # psit2 - norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdpsit + # ) + + # plt.figure(figsize=(10, 6)) + # plt.plot( + # t_anual[t_ord], + # psit2[t_ord], + # label=r"$\psi_t$", + # linewidth=2, + # color=self.colors[0], + # ) + # plt.fill_between( + # t_anual[t_ord], + # ci_low[t_ord], + # ci_up[t_ord], + # color=self.colors[0], + # alpha=0.3, + # label=r"$\psi_t$ Confidence Interval", + # ) + # # plt.plot(t_anual[t_ord], quan95[t_ord], linestyle='dashed', color=self.colors[2], markersize=5, label=fr"$q_{self.quanval}$") + # plt.title(f"Harmonic part of Scale parameter ({self.var_name})") + # plt.xlabel("Time (yearly scale)") + # plt.xticks(month_positions, month_initials) + # plt.ylabel(r"$\psi_t$") + # plt.grid(True) + # if save: + # plt.savefig( + # f"Figures/Harmonic_Scale_Parameter_{self.var_name}.png", dpi=300 + # ) + # plt.show() #### PP Plot self.PPplot(save=save) From 9c03395b633858e232caf363b4874120a6cfdf48 Mon Sep 17 00:00:00 2001 From: Colladito Date: Mon, 6 Oct 2025 12:55:21 +0200 Subject: [PATCH 03/18] [VCF] Add summary and paramplot method --- bluemath_tk/distributions/nonstat_gev.py | 179 ++++++++++++++++++++++- 1 file changed, 177 insertions(+), 2 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index 912dc1b..06b79af 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -1157,8 +1157,6 @@ def auto_adjust(self, max_iter: int = 1000, plot: bool = False) -> dict: + self.ntrend_sh + self.nind_sh ) - fit_result["n_params"] = n_params - fit_result["AIC"] = self._AIC(-fit_result["negloglikelihood"], n_params) # Compute the final loglikelihood and the information matrix f, Jx, Hxx = self._loglikelihood( @@ -1505,6 +1503,8 @@ def _fit( fit_result["x"] = result.x # Optimal parameters vector fit_result["negloglikelihood"] = result.fun # Optimal loglikelihood + fit_result["AIC"] = self._AIC(-fit_result["negloglikelihood"], n_params) + fit_result["n_params"] = n_params fit_result["success"] = result.success fit_result["message"] = result.message fit_result["grad"] = result.grad @@ -1632,6 +1632,8 @@ def _fit( else: fit_result["gamma_cov"] = np.empty(0) + self.fit_result = fit_result + return fit_result def _auxmin_loglikelihood( @@ -5412,6 +5414,9 @@ def plot(self, return_plot: bool = True, save: bool = False, init_year=2000): #### QQ plot self.QQplot(save=save) + #### Parameters Heatmap plot + # self.paramplot(save=save) + #### Return periods if ( self.ntrend_loc == 0 @@ -5422,6 +5427,93 @@ def plot(self, return_plot: bool = True, save: bool = False, init_year=2000): ) and return_plot: self.ReturnPeriodPlot() + def paramplot(self, save: bool=False): + """ + Create a heatmap of parameter sensitivities for location, scale and shape. + """ + # Get parameters + loc_params = [] + param_names = [] + if self.fit_result["beta0"] is not None: + loc_params.append(self.fit_result["beta0"]) + param_names.append("Intercept") + if self.fit_result["beta"].size > 0: + for i, b in enumerate(self.fit_result["beta"]): + loc_params.append(b) + param_names.append(f"PC{i}") + if self.fit_result["beta_cov"].size > 0: + for i, b in enumerate(self.fit_result["beta_cov"]): + loc_params.append(b) + param_names.append(f"{self.covariates.columns[i]}") + if self.fit_result["betaT"].size > 0: + loc_params.append(self.fit_result["betaT"]) + param_names.append("Trend") + + # Scale parameters + scale_params = [] + if self.fit_result["alpha0"] is not None: + scale_params.append(self.fit_result["alpha0"]) + if self.fit_result["alpha"].size > 0: + scale_params.extend(self.fit_result["alpha"]) + if self.fit_result["alpha_cov"].size > 0: + scale_params.extend(self.fit_result["alpha_cov"]) + if self.fit_result["alphaT"].size > 0: + scale_params.append(self.fit_result["alphaT"]) + + # Shape parameters + shape_params = [] + if self.fit_result["gamma0"] is not None: + shape_params.append(self.fit_result["gamma0"]) + if self.fit_result["gamma"].size > 0: + shape_params.extend(self.fit_result["gamma"]) + if self.fit_result["gamma_cov"].size > 0: + shape_params.extend(self.fit_result["gamma_cov"]) + if self.fit_result["gammaT"].size > 0: + shape_params.append(self.fit_result["gammaT"]) + + # Create data matrix + max_len = max(len(param_names), len(scale_params), len(shape_params)) + data = np.zeros((3, len(param_names))) + data[0,:len(loc_params)] = loc_params + data[1,:len(scale_params)] = scale_params + data[2,:len(shape_params)] = shape_params + + # Create figure + fig, ax = plt.subplots(figsize=(20,4)) + + # Create heatmap + im = ax.imshow(data, cmap='RdBu', aspect='auto', vmin=-2, vmax=2) + + # Add colorbar + cbar = plt.colorbar(im) + cbar.set_label('Coefficient value') + + # Add text annotations + for i in range(data.shape[0]): + for j in range(data.shape[1]): + text = ax.text(j, i, f'{data[i, j]:.2f}', + ha="center", va="center", color="black") + + # Configure axes + ax.set_xticks(np.arange(len(param_names))) + ax.set_yticks(np.arange(3)) + ax.set_xticklabels(param_names) + ax.set_yticklabels(['Location', 'Scale', 'Shape']) + + # Rotate x labels for better readability + plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor") + + # Add title and adjust layout + ax.set_title('Time-dependent GEV Parameters') + ax.set_xlabel('Components') + + plt.tight_layout() + + if save: + plt.savefig(f"Figures/Parameters_Heatmap_{self.var_name}.png", dpi=300, bbox_inches='tight') + + plt.show() + def QQplot(self, save: bool = False): """ QQ plot @@ -6950,3 +7042,86 @@ def _ConfidInterQuanAggregate(self, q, t0, t1) -> np.ndarray: stdQuan = np.sqrt(jacob.T @ self.invI0 @ jacob) return stdQuan + + def summary(self): + """ + Print a summary of the fitted model, including parameter estimates, standard errors and fit statistics. + """ + + std_params = np.sqrt(np.diag(self.invI0)) + param_idx = 0 + + print(f"\nFitted Time-Dependent GEV model for {self.var_name}") + print("="*50) + print("\nLocation Parameters") + print("-"*20) + print(f"Beta0: {self.beta0:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + + # Harmonic terms for location + for i in range(self.nmu): + print(f"Beta{i+1} (sin): {self.beta[2*i]:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + print(f"Beta{i+1} (cos): {self.beta[2*i+1]:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + + # Trend in location + if self.ntrend_loc > 0: + print(f"BetaT: {self.betaT:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + + # Covariates in location + for i in range(self.nind_loc): + print(f"Beta_cov{i+1}: {self.beta_cov[i]:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + + print("\nScale Parameters") + print("-"*20) + print(f"Alpha0: {self.alpha0:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + + # Harmonic terms for scale + for i in range(self.npsi): + print(f"Alpha{i+1} (sin): {self.alpha[2*i]:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + print(f"Alpha{i+1} (cos): {self.alpha[2*i+1]:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + + # Trend in scale + if self.ntrend_sc > 0: + print(f"AlphaT: {self.alphaT:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + + # Covariates in scale + for i in range(self.nind_sc): + print(f"Alpha_cov{i+1}: {self.alpha_cov[i]:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + + print("\nShape Parameters") + print("-"*20) + if self.ngamma0 > 0: + print(f"Gamma0: {self.gamma0:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + + # Harmonic terms for shape + for i in range(self.ngamma): + print(f"Gamma{i+1} (sin): {self.gamma[2*i]:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + print(f"Gamma{i+1} (cos): {self.gamma[2*i+1]:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + + # Trend in shape + if self.ntrend_sh > 0: + print(f"GammaT: {self.gammaT:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + + # Covariates in shape + for i in range(self.nind_sh): + print(f"Gamma_cov{i+1}: {self.gamma_cov[i]:.4f} ({std_params[param_idx]:.4f})") + param_idx += 1 + + print("\nFit Statistics") + print("-"*20) + print(f"Log-likelihood: {-self.fit_result['negloglikelihood']:.4f}") + print(f"AIC: {self.fit_result['AIC']:.4f}") + print(f"Number of parameters: {self.fit_result['n_params']}") \ No newline at end of file From 1030bf4e7d26a3efefcaf9b099f3f3037bb533c7 Mon Sep 17 00:00:00 2001 From: Colladito Date: Mon, 6 Oct 2025 13:45:25 +0200 Subject: [PATCH 04/18] [VCF] Improve summary in NonStatGEV --- bluemath_tk/distributions/nonstat_gev.py | 86 ++++++++++++++---------- 1 file changed, 51 insertions(+), 35 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index 06b79af..4274d34 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -7047,81 +7047,97 @@ def summary(self): """ Print a summary of the fitted model, including parameter estimates, standard errors and fit statistics. """ - std_params = np.sqrt(np.diag(self.invI0)) param_idx = 0 + z_norm = norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) print(f"\nFitted Time-Dependent GEV model for {self.var_name}") - print("="*50) + print("="*70) + + # Header format + param_header = "Parameter".ljust(15) + estimate_header = "Estimate".rjust(12) + se_header = "Std. Error".rjust(15) + ci_header = f"{self.quanval*100:.0f}% CI".center(25) + header = f"{param_header} {estimate_header} {se_header} {ci_header}" + print("\nLocation Parameters") - print("-"*20) - print(f"Beta0: {self.beta0:.4f} ({std_params[param_idx]:.4f})") - param_idx += 1 + print("-"*70) + print(header) + print("-"*70) + + # Parameter line format + def format_line(name, value, std_err): + ci_low = value - std_err*z_norm + ci_up = value + std_err*z_norm + return f"{name:<15} {value:>12.4f} {std_err:>15.4f} {f'[{ci_low:>8.4f},{ci_up:>8.4f}]':>25}" - # Harmonic terms for location + # Location parameters + print(format_line("Beta0", self.beta0, std_params[param_idx])) + param_idx += 1 + for i in range(self.nmu): - print(f"Beta{i+1} (sin): {self.beta[2*i]:.4f} ({std_params[param_idx]:.4f})") + print(format_line(f"Beta{i+1} (sin)", self.beta[2*i], std_params[param_idx])) param_idx += 1 - print(f"Beta{i+1} (cos): {self.beta[2*i+1]:.4f} ({std_params[param_idx]:.4f})") + print(format_line(f"Beta{i+1} (cos)", self.beta[2*i+1], std_params[param_idx])) param_idx += 1 - # Trend in location if self.ntrend_loc > 0: - print(f"BetaT: {self.betaT:.4f} ({std_params[param_idx]:.4f})") + print(format_line("BetaT", self.betaT, std_params[param_idx])) param_idx += 1 - # Covariates in location for i in range(self.nind_loc): - print(f"Beta_cov{i+1}: {self.beta_cov[i]:.4f} ({std_params[param_idx]:.4f})") + print(format_line(f"Beta_cov{i+1}", self.beta_cov[i], std_params[param_idx])) param_idx += 1 print("\nScale Parameters") - print("-"*20) - print(f"Alpha0: {self.alpha0:.4f} ({std_params[param_idx]:.4f})") + print("-"*70) + print(header) + print("-"*70) + + print(format_line("Alpha0", self.alpha0, std_params[param_idx])) param_idx += 1 - # Harmonic terms for scale for i in range(self.npsi): - print(f"Alpha{i+1} (sin): {self.alpha[2*i]:.4f} ({std_params[param_idx]:.4f})") + print(format_line(f"Alpha{i+1} (sin)", self.alpha[2*i], std_params[param_idx])) param_idx += 1 - print(f"Alpha{i+1} (cos): {self.alpha[2*i+1]:.4f} ({std_params[param_idx]:.4f})") + print(format_line(f"Alpha{i+1} (cos)", self.alpha[2*i+1], std_params[param_idx])) param_idx += 1 - # Trend in scale if self.ntrend_sc > 0: - print(f"AlphaT: {self.alphaT:.4f} ({std_params[param_idx]:.4f})") + print(format_line("AlphaT", self.alphaT, std_params[param_idx])) param_idx += 1 - # Covariates in scale for i in range(self.nind_sc): - print(f"Alpha_cov{i+1}: {self.alpha_cov[i]:.4f} ({std_params[param_idx]:.4f})") + print(format_line(f"Alpha_cov{i+1}", self.alpha_cov[i], std_params[param_idx])) param_idx += 1 print("\nShape Parameters") - print("-"*20) + print("-"*70) + print(header) + print("-"*70) + if self.ngamma0 > 0: - print(f"Gamma0: {self.gamma0:.4f} ({std_params[param_idx]:.4f})") + print(format_line("Gamma0", self.gamma0, std_params[param_idx])) param_idx += 1 - # Harmonic terms for shape for i in range(self.ngamma): - print(f"Gamma{i+1} (sin): {self.gamma[2*i]:.4f} ({std_params[param_idx]:.4f})") + print(format_line(f"Gamma{i+1} (sin)", self.gamma[2*i], std_params[param_idx])) param_idx += 1 - print(f"Gamma{i+1} (cos): {self.gamma[2*i+1]:.4f} ({std_params[param_idx]:.4f})") + print(format_line(f"Gamma{i+1} (cos)", self.gamma[2*i+1], std_params[param_idx])) param_idx += 1 - # Trend in shape if self.ntrend_sh > 0: - print(f"GammaT: {self.gammaT:.4f} ({std_params[param_idx]:.4f})") + print(format_line("GammaT", self.gammaT, std_params[param_idx])) param_idx += 1 - # Covariates in shape for i in range(self.nind_sh): - print(f"Gamma_cov{i+1}: {self.gamma_cov[i]:.4f} ({std_params[param_idx]:.4f})") + print(format_line(f"Gamma_cov{i+1}", self.gamma_cov[i], std_params[param_idx])) param_idx += 1 - + print("\nFit Statistics") - print("-"*20) - print(f"Log-likelihood: {-self.fit_result['negloglikelihood']:.4f}") - print(f"AIC: {self.fit_result['AIC']:.4f}") - print(f"Number of parameters: {self.fit_result['n_params']}") \ No newline at end of file + print("-"*70) + stats_width = 30 + print(f"{'Log-likelihood:':<{stats_width}} {-self.fit_result['negloglikelihood']:>.4f}") + print(f"{'AIC:':<{stats_width}} {self.fit_result['AIC']:>.4f}") + print(f"{'Number of parameters:':<{stats_width}} {self.fit_result['n_params']}") \ No newline at end of file From 7bd823794ba17d6575d52ffaf781a0cf01919f1d Mon Sep 17 00:00:00 2001 From: Colladito Date: Mon, 6 Oct 2025 17:01:26 +0200 Subject: [PATCH 05/18] [VCF] NonStatGEV return periods with bsimp --- bluemath_tk/distributions/nonstat_gev.py | 161 +++++++++++++++++++---- 1 file changed, 132 insertions(+), 29 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index 4274d34..1f51f15 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -4594,7 +4594,7 @@ def _parametro( y = y + beta_cov[i] * indicesint[i] else: for i in range(nind): - indicesintaux = self._search(times, covariates[:, i]) + indicesintaux = self._search(times, covariates[:, i], x) y = y + beta_cov[i] * indicesintaux else: for i in range(nind): @@ -4602,7 +4602,7 @@ def _parametro( return y - def _search(self, times: np.ndarray, values: np.ndarray) -> np.ndarray: + def _search(self, times: np.ndarray, values: np.ndarray, xs) -> np.ndarray: """ Function to search the nearest value of certain time to use in self._parametro function @@ -4614,12 +4614,12 @@ def _search(self, times: np.ndarray, values: np.ndarray) -> np.ndarray: Values of the covariates at those times """ n = times.shape[0] - yin = np.zeros_like(self.t) + yin = np.zeros_like(xs) pos = 0 - for j in range(len(self.t)): + for j in range(xs.size): found = 0 while found == 0 and pos < n: - if self.t[j] < times[pos]: + if xs[j] < times[pos]: yin[j] = values[pos] found = 1 else: @@ -6482,7 +6482,7 @@ def _aggquantile( gamma_cov=None, ) -> np.ndarray: """ - Function to compute the aggregated quantile for certain parameters + Function to compute the aggregated quantile between two time stamps Parameters ---------- @@ -6566,7 +6566,7 @@ def _aggquantile( cov_locint = np.zeros(len(beta_cov)) cov_scint = np.zeros(len(alpha_cov)) cov_shint = np.zeros(len(gamma_cov)) - if len(pos) > 0: + if pos.size: for i in range(len(beta_cov)): cov_locint[i] = np.mean( self.covariates.iloc[pos, self.list_loc[i]].values @@ -6580,14 +6580,14 @@ def _aggquantile( self.covariates.iloc[pos, self.list_sh[i]].values ) else: - cov_locint = None - cov_scint = None - cov_shint = None + cov_locint = np.zeros(len(beta_cov)) + cov_scint = np.zeros(len(alpha_cov)) + cov_shint = np.zeros(len(gamma_cov)) # Require quantile zqout = np.zeros(m) - media = quad( + media = bsimp( lambda x: self._parametro( beta0, beta, @@ -6600,7 +6600,7 @@ def _aggquantile( ), 0, 1, - )[0] + ) # std = quad( # lambda x: np.exp( # self._parametro( @@ -6626,7 +6626,7 @@ def _aggquantile( integ: float = 0 while err > 1e-4 and iter1 < 1000: zqold = zq - integ = quad( + integ = bsimp( lambda x: self._fzeroquanint( x, zqold, @@ -6646,12 +6646,14 @@ def _aggquantile( beta_cov, alpha_cov, gamma_cov, + self.t, + self.kt ), float(t0[il]), float(t1[il]), - )[0] + ) integ += np.log(q[il]) / 12 - dint = quad( + dint = bsimp( lambda x: self._fzeroderiquanint( x, zqold, @@ -6671,10 +6673,12 @@ def _aggquantile( beta_cov, alpha_cov, gamma_cov, + self.t, + self.kt ), float(t0[il]), float(t1[il]), - )[0] + ) zq += -integ / dint if np.abs(zq) > 1e-5: err = np.abs((zq - zqold) / zqold) @@ -6711,6 +6715,8 @@ def _fzeroquanint( beta_cov, alpha_cov, gamma_cov, + times=None, + ktold=None ) -> np.ndarray: """ Auxiliar function to solve the quantile @@ -6719,6 +6725,8 @@ def _fzeroquanint( ------ zn : np.ndarray """ + if ktold is None: + ktold=self.kt # Evaluate the location parameter at each time t as a function of the actual values of the parameters given by p mut1 = self._parametro( @@ -6763,16 +6771,17 @@ def _fzeroquanint( # The corresponding GUMBEl values are set to 1 to avoid numerical problems, note that for those cases the GUMBEL expressions are used epst[posG] = 1 - # TODO: AÑadir - #### if times is not None: - #### kt2 = spline + if times is not None: + kt2 = np.interp(t, np.asarray(self.t, float), np.asarray(ktold, float)) + else: + kt2 = np.ones_like(mut1) mut = mut1 psit = psit1 - mut[pos] = mut1[pos] + psit1[pos] * (self.kt[pos] ** epst[pos] - 1) / epst[pos] - psit[pos] = psit1[pos] * self.kt[pos] ** epst[pos] + mut[pos] = mut1[pos] + psit1[pos] * (kt2[pos] ** epst[pos] - 1) / epst[pos] + psit[pos] = psit1[pos] * kt2[pos] ** epst[pos] # Modify the parameters to include Gumbel - mut[posG] += psit[posG] * np.log(self.kt[posG]) + mut[posG] += psit[posG] * np.log(kt2[posG]) # Evaluate the auxiliary variable xn = (zq - mut) / psit @@ -6805,6 +6814,8 @@ def _fzeroderiquanint( beta_cov, alpha_cov, gamma_cov, + times=None, + ktold=None, ): """ Auxiliar Function to solve the quantile derivative @@ -6813,6 +6824,9 @@ def _fzeroderiquanint( ------ zn : np.ndarray """ + if ktold is None: + ktold = self.kt + # Evaluate the location parameter at each time t as a function of the actual values of the parameters given by p mut1 = self._parametro( beta0, @@ -6856,16 +6870,17 @@ def _fzeroderiquanint( # The corresponding GUMBEl values are set to 1 to avoid numerical problems, note that for those cases the GUMBEL expressions are used epst[posG] = 1 - # TODO: AÑadir - #### if times is not None: - #### kt2 = spline + if times is not None: + kt2 = np.interp(t, np.asarray(self.t, float), np.asarray(ktold, float)) + else: + kt2 = np.ones_like(mut1) mut = mut1 psit = psit1 - mut[pos] = mut1[pos] + psit1[pos] * (self.kt[pos] ** epst[pos] - 1) / epst[pos] - psit[pos] = psit1[pos] * self.kt[pos] ** epst[pos] + mut[pos] = mut1[pos] + psit1[pos] * (kt2[pos] ** epst[pos] - 1) / epst[pos] + psit[pos] = psit1[pos] * kt2[pos] ** epst[pos] # Modify the parameters to include Gumbel - mut[posG] += psit[posG] * np.log(self.kt[posG]) + mut[posG] += psit[posG] * np.log(kt2[posG]) # Evaluate the auxiliary variable xn = (zq - mut) / psit @@ -7140,4 +7155,92 @@ def format_line(name, value, std_err): stats_width = 30 print(f"{'Log-likelihood:':<{stats_width}} {-self.fit_result['negloglikelihood']:>.4f}") print(f"{'AIC:':<{stats_width}} {self.fit_result['AIC']:>.4f}") - print(f"{'Number of parameters:':<{stats_width}} {self.fit_result['n_params']}") \ No newline at end of file + print(f"{'Number of parameters:':<{stats_width}} {self.fit_result['n_params']}") + + + + +def bsimp(fun, a, b, n=None, epsilon=1e-8, trace=0): + """ + BSIMP Numerically evaluate integral, low order method. + I = BSIMP('F',A,B) approximates the integral of F(X) from A to B + within a relative error of 1e-3 using an iterative + Simpson's rule. 'F' is a string containing the name of the + function. Function F must return a vector of output values if given + a vector of input values.% + I = BSIMP('F',A,B,EPSILON) integrates to a total error of EPSILON. % + I = BSIMP('F',A,B,N,EPSILON,TRACE,TRACETOL) integrates to a + relative error of EPSILON, + beginning with n subdivisions of the interval [A,B],for non-zero + TRACE traces the function + evaluations with a point plot. + [I,cnt] = BSIMP(F,a,b,epsilon) also returns a function evaluation count.% + Roberto Minguez Solana% Copyright (c) 2001 by Universidad de Cantabria + """ + if n is None: + n = int(365*8*(b-a)) + + # The number of initial subintervals must be pair + if n % 2 != 0: + n = n + 1 + + # Step 1: + h = (b-a)/n + + # Step 3: + x = np.linspace(a, b, n+1) + y = fun(x) + + # print("y size: ", y.size) + # print("n: ", n) + # TODO: + # if trace: + + auxI = 0 + ainit = y[0] + y[n] + auxI1 = 0 + auxI2 = 0 + + auxI1 = auxI1 + sum(y[1:n:2]) + auxI2 = auxI2 + sum(y[2:n:2]) + + cnt = n + + # Step 4 + integral = (ainit + 4 * auxI1 + 2*auxI2) * h/3 + auxtot = auxI1 + auxI2 + + # Step 5 + integral1 = integral + epsilon*2 + j = 2 + + # Step 6 + error = 1 + while error > epsilon: + cnt = cnt + int(j*n//2) + # Step 7 + integral1 = integral + # Step 8 + x = np.linspace(a+h/j, b-h/j, int(j*n//2)) + y = fun(x) + + auxI = auxI + sum(y[0:int(j*n//2)]) + + # Step 9 + integral = (ainit + 4 * auxI + 2*auxtot) * h/(3*j) + + # Step 10 + j = j*2 + + # Step 11 + auxtot = auxtot + auxI + auxI = 0 + # Error + if np.abs(integral1) > 1e-5: + error = np.abs((integral-integral1)/integral1) + else: + error = np.abs(integral-integral1) + + + return integral + \ No newline at end of file From a657d448f1d30f014c5b30f797808ea1371dd5c3 Mon Sep 17 00:00:00 2001 From: Colladito Date: Mon, 6 Oct 2025 18:25:08 +0200 Subject: [PATCH 06/18] [VCF] Improve ReturnPeriodPlot and make it faster --- bluemath_tk/distributions/nonstat_gev.py | 265 ++++++++++++++--------- 1 file changed, 157 insertions(+), 108 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index 1f51f15..31c0a42 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -5,7 +5,7 @@ import numpy as np import pandas as pd from scipy.integrate import quad -from scipy.optimize import minimize +from scipy.optimize import minimize, root_scalar from scipy.stats import norm from ..core.models import BlueMathModel @@ -4594,7 +4594,7 @@ def _parametro( y = y + beta_cov[i] * indicesint[i] else: for i in range(nind): - indicesintaux = self._search(times, covariates[:, i], x) + indicesintaux = self._search(times, covariates[:, i], np.asarray(x).flatten()) y = y + beta_cov[i] * indicesintaux else: for i in range(nind): @@ -6347,7 +6347,7 @@ def _CDFGEVt(self): return F - def ReturnPeriodPlot(self, annualplot=True): + def ReturnPeriodPlot(self, annualplot=True, save=False): """ Funtion to plot the Aggregated Return period plot for each month and the annual Return period @@ -6357,10 +6357,10 @@ def ReturnPeriodPlot(self, annualplot=True): Whether to plot the annual return period plot """ - # Ts = np.array([2, 5, 10, 20, 25, 50, 75, 100, 200, 300, 400, 500]) - Ts = np.concatenate( - (np.arange(2, 10, 1), np.arange(10, 100, 10), np.arange(100, 501, 100)) - ) + Ts = np.array([1.1, 1.5, 2, 3, 4, 5, 7.5, 10, 15, 20, 30, 40, 50, 75, 100, 150, 200, 500, 1000]) + # Ts = np.concatenate( + # (np.arange(2, 10, 1), np.arange(10, 100, 10), np.arange(100, 501, 100)) + # ) nts = len(Ts) quanaggrA = np.zeros(nts) @@ -6371,9 +6371,10 @@ def ReturnPeriodPlot(self, annualplot=True): quanaggr[i, j] = self._aggquantile(1 - 1 / Ts[j], i / 12, (i + 1) / 12)[ 0 ] - stdQuan = self._ConfidInterQuanAggregate( - 1 - 1 / Ts[j], i / 12, (i + 1) / 12 - ) + # stdQuan = self._ConfidInterQuanAggregate( + # 1 - 1 / Ts[j], i / 12, (i + 1) / 12 + # ) + stdQuan = 0.1 stdDqX[i, j] = stdQuan * norm.ppf( 1 - (1 - self.quanval) / 2, loc=0, scale=1 ) @@ -6386,8 +6387,8 @@ def ReturnPeriodPlot(self, annualplot=True): stdup = np.zeros(nts) stdlo = np.zeros(nts) for i in range(nts): - stdQuan = self._ConfidInterQuanAggregate(1 - 1 / Ts[i], 0, 1) - # stdQuan = 0.1 + # stdQuan = self._ConfidInterQuanAggregate(1 - 1 / Ts[i], 0, 1) + stdQuan = 0.1 stdup[i] = quanaggrA[i] + stdQuan * norm.ppf( 1 - (1 - self.quanval) / 2, loc=0, scale=1 ) @@ -6425,9 +6426,10 @@ def ReturnPeriodPlot(self, annualplot=True): "#FF3333", "#33FF33", ] - plt.figure(figsize=(10, 6)) + fig = plt.figure(figsize=(10, 6)) + ax = fig.add_subplot() for i in range(12): - plt.semilogx( + ax.semilogx( Ts, quanaggr[i, :], color=colors[i], @@ -6438,7 +6440,7 @@ def ReturnPeriodPlot(self, annualplot=True): # Anual return periods if annualplot: - plt.semilogx(Ts, quanaggrA, color="black", linewidth=2, label="Annual") + ax.semilogx(Ts, quanaggrA, color="black", linewidth=2, label="Annual") ny = int(np.ceil(self.t[-1])) hmax1 = np.zeros(ny) for j in range(ny): @@ -6449,18 +6451,22 @@ def ReturnPeriodPlot(self, annualplot=True): ProHsmaxsort = np.arange(1, len(hmaxsort) + 1) / (len(hmaxsort) + 1) Tapprox = 1 / (1 - ProHsmaxsort) idx = np.where(Tapprox >= 2)[0] - plt.semilogx(Tapprox[idx], hmaxsort[idx], "ok", markersize=1.6) - plt.semilogx(Ts, stdlo, "--k", linewidth=1.1) - plt.semilogx(Ts, stdup, "--k", linewidth=1.1) - plt.xlabel("Return Period (years)") - plt.ylabel(r"$x$") - plt.legend(loc="center left", bbox_to_anchor=(1, 0.5)) - plt.xticks([1, 2, 5, 10, 20, 50, 100, 250, 500]) - plt.xlim(left=1.8, right=Ts[-1] + 50) - plt.ylim(bottom=0) - plt.title(f"Aggregate Quantiles ({self.var_name})") - plt.grid(True) - plt.margins(x=0.1) + ax.semilogx(Tapprox[idx], hmaxsort[idx], "ok", markersize=1.6) + ax.semilogx(Ts, stdlo, "--k", linewidth=1.1) + ax.semilogx(Ts, stdup, "--k", linewidth=1.1) + + ax.set_xlabel("Return Period (Years)") + ax.set_ylabel(f"{self.var_name}") + ax.legend(loc="center left", bbox_to_anchor=(1, 0.5)) + ax.set_xticks([1, 2, 5, 10, 20, 50, 100, 250, 500]) + ax.set_xticklabels([1, 2, 5, 10, 20, 50, 100, 250, 500]) + ax.set_xlim(left=1.8, right=Ts[-1] + 50) + ax.set_ylim(bottom=0) + ax.set_title(f"Aggregate Quantiles ({self.var_name})") + ax.grid(True) + plt.tight_layout() + if save: + plt.savefig(f"Figures/ReturnPeriod_{self.var_name}.png", dpi=300) plt.show() def _aggquantile( @@ -6545,6 +6551,12 @@ def _aggquantile( if gamma_cov is None: gamma_cov = self.gamma_cov + + # Deal with scalars and vectors + beta_cov = np.atleast_1d(beta_cov) + alpha_cov = np.atleast_1d(alpha_cov) + gamma_cov = np.atleast_1d(gamma_cov) + q = np.array([q]) t0 = np.array([t0]) t1 = np.array([t1]) @@ -6563,31 +6575,31 @@ def _aggquantile( # For the required period the mean value of the corresponding covariates is calculated and considered constant for the rest of the study if len(self.t) > 0: pos = np.where((self.t >= t0) & (self.t <= t1))[0] - cov_locint = np.zeros(len(beta_cov)) - cov_scint = np.zeros(len(alpha_cov)) - cov_shint = np.zeros(len(gamma_cov)) + cov_locint = np.zeros_like(beta_cov) + cov_scint = np.zeros_like(alpha_cov) + cov_shint = np.zeros_like(gamma_cov) if pos.size: - for i in range(len(beta_cov)): + for i in range(beta_cov.size): cov_locint[i] = np.mean( self.covariates.iloc[pos, self.list_loc[i]].values ) - for i in range(len(alpha_cov)): + for i in range(alpha_cov.size): cov_scint[i] = np.mean( self.covariates.iloc[pos, self.list_sc[i]].values ) - for i in range(len(gamma_cov)): + for i in range(gamma_cov.size): cov_shint[i] = np.mean( self.covariates.iloc[pos, self.list_sh[i]].values ) else: - cov_locint = np.zeros(len(beta_cov)) - cov_scint = np.zeros(len(alpha_cov)) - cov_shint = np.zeros(len(gamma_cov)) + cov_locint = np.zeros(beta_cov.size) + cov_scint = np.zeros(alpha_cov.size) + cov_shint = np.zeros(gamma_cov.size) # Require quantile zqout = np.zeros(m) - media = bsimp( + media,_ = quad( lambda x: self._parametro( beta0, beta, @@ -6618,80 +6630,117 @@ def _aggquantile( # 1, # )[0] + # CHANGE THIS + # for il in range(m): + # # for jl in range(n) + # zq = media + # err: float = 1 + # iter1: int = 1 + # integ: float = 0 + # while err > 1e-4 and iter1 < 1000: + # zqold = zq + # integ,_ = quad( + # lambda x: self._fzeroquanint( + # x, + # zqold, + # q[il], + # cov_locint, + # cov_scint, + # cov_shint, + # beta0, + # beta, + # alpha0, + # alpha, + # gamma0, + # gamma, + # betaT, + # alphaT, + # gammaT, + # beta_cov, + # alpha_cov, + # gamma_cov, + # self.t, + # self.kt + # ), + # float(t0[il]), + # float(t1[il]), + # ) + # integ = integ + np.log(q[il]) / 12 + # dint,_ = quad( + # lambda x: self._fzeroderiquanint( + # x, + # zqold, + # q[il], + # cov_locint, + # cov_scint, + # cov_shint, + # beta0, + # beta, + # alpha0, + # alpha, + # gamma0, + # gamma, + # betaT, + # alphaT, + # gammaT, + # beta_cov, + # alpha_cov, + # gamma_cov, + # self.t, + # self.kt + # ), + # float(t0[il]), + # float(t1[il]), + # ) + # zq = zq -integ / dint + # if np.abs(zq) > 1e-5: + # err = np.abs((zq - zqold) / zqold) + # else: + # err = np.abs(zq - zqold) + # iter1 += 1 + # if iter1 == 1000: + # zq = np.nan + # Warning("Maximum number of Newton iterations") + # if integ > 1e-2: + # zq = np.nan + # Warning("False zero, check it") + # zqout[il] = zq + + # return zqout + for il in range(m): - # for jl in range(n) - zq = media - err: float = 1 - iter1: int = 1 - integ: float = 0 - while err > 1e-4 and iter1 < 1000: - zqold = zq - integ = bsimp( + # function of z whose root we want + def F(z): + integ, _ = quad( lambda x: self._fzeroquanint( - x, - zqold, - q[il], - cov_locint, - cov_scint, - cov_shint, - beta0, - beta, - alpha0, - alpha, - gamma0, - gamma, - betaT, - alphaT, - gammaT, - beta_cov, - alpha_cov, - gamma_cov, - self.t, - self.kt - ), - float(t0[il]), - float(t1[il]), + x, z, q[il], + cov_locint, cov_scint, cov_shint, + beta0, beta, alpha0, alpha, + gamma0, gamma, betaT, alphaT, gammaT, + beta_cov, alpha_cov, gamma_cov, + self.t, self.kt), + float(t0[il]), float(t1[il]) ) - integ += np.log(q[il]) / 12 - dint = bsimp( - lambda x: self._fzeroderiquanint( - x, - zqold, - q[il], - cov_locint, - cov_scint, - cov_shint, - beta0, - beta, - alpha0, - alpha, - gamma0, - gamma, - betaT, - alphaT, - gammaT, - beta_cov, - alpha_cov, - gamma_cov, - self.t, - self.kt - ), - float(t0[il]), - float(t1[il]), + return integ + np.log(q[il]) / 12.0 + + try: + # root finding: start near `media`, bracket if possible + sol = root_scalar( + F, + x0=media, x1=media + 1.0, # secant starting points + method="secant", + xtol=1e-6, rtol=1e-6, + maxiter=200 ) - zq += -integ / dint - if np.abs(zq) > 1e-5: - err = np.abs((zq - zqold) / zqold) + if sol.converged: + if abs(F(sol.root)) < 1e-2: + zqout[il] = sol.root + else: + zqout[il] = np.nan # "False zero" check else: - err = np.abs(zq - zqold) - iter1 += 1 - if iter1 == 1000: - zq = np.nan - Warning("Maximum number of Newton iterations") - if integ > 1e-2: - zq = np.nan - Warning("False zero, check it") - zqout[il] = zq + zqout[il] = np.nan + except Exception: + zqout[il] = np.nan return zqout @@ -6772,7 +6821,7 @@ def _fzeroquanint( epst[posG] = 1 if times is not None: - kt2 = np.interp(t, np.asarray(self.t, float), np.asarray(ktold, float)) + kt2 = np.interp(t, np.asarray(self.t, float), np.asarray(ktold, float)).flatten() else: kt2 = np.ones_like(mut1) @@ -6871,7 +6920,7 @@ def _fzeroderiquanint( epst[posG] = 1 if times is not None: - kt2 = np.interp(t, np.asarray(self.t, float), np.asarray(ktold, float)) + kt2 = np.interp(t, np.asarray(self.t, float), np.asarray(ktold, float)).flatten() else: kt2 = np.ones_like(mut1) From 4aceeb6dbb808d9c908ae19826100ab0fa262616 Mon Sep 17 00:00:00 2001 From: Colladito Date: Mon, 6 Oct 2025 22:08:03 +0200 Subject: [PATCH 07/18] [VCF] Update aggregated quantiles in NonStatGEV --- bluemath_tk/distributions/nonstat_gev.py | 87 ++---------------------- 1 file changed, 5 insertions(+), 82 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index 31c0a42..5c8f856 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -6450,8 +6450,9 @@ def ReturnPeriodPlot(self, annualplot=True, save=False): hmaxsort = np.sort(hmax1) ProHsmaxsort = np.arange(1, len(hmaxsort) + 1) / (len(hmaxsort) + 1) Tapprox = 1 / (1 - ProHsmaxsort) - idx = np.where(Tapprox >= 2)[0] - ax.semilogx(Tapprox[idx], hmaxsort[idx], "ok", markersize=1.6) + # idx = np.where(Tapprox >= 2)[0] + # ax.semilogx(Tapprox[idx], hmaxsort[idx], "ok", markersize=1.6) + ax.semilogx(Tapprox, hmaxsort, "ok", markersize=2) ax.semilogx(Ts, stdlo, "--k", linewidth=1.1) ax.semilogx(Ts, stdup, "--k", linewidth=1.1) @@ -6460,9 +6461,9 @@ def ReturnPeriodPlot(self, annualplot=True, save=False): ax.legend(loc="center left", bbox_to_anchor=(1, 0.5)) ax.set_xticks([1, 2, 5, 10, 20, 50, 100, 250, 500]) ax.set_xticklabels([1, 2, 5, 10, 20, 50, 100, 250, 500]) - ax.set_xlim(left=1.8, right=Ts[-1] + 50) + ax.set_xlim(left=0.9, right=Ts[-1] + 50) ax.set_ylim(bottom=0) - ax.set_title(f"Aggregate Quantiles ({self.var_name})") + # ax.set_title(f"Aggregate Quantiles ({self.var_name})") ax.grid(True) plt.tight_layout() if save: @@ -6630,84 +6631,6 @@ def _aggquantile( # 1, # )[0] - # CHANGE THIS - # for il in range(m): - # # for jl in range(n) - # zq = media - # err: float = 1 - # iter1: int = 1 - # integ: float = 0 - # while err > 1e-4 and iter1 < 1000: - # zqold = zq - # integ,_ = quad( - # lambda x: self._fzeroquanint( - # x, - # zqold, - # q[il], - # cov_locint, - # cov_scint, - # cov_shint, - # beta0, - # beta, - # alpha0, - # alpha, - # gamma0, - # gamma, - # betaT, - # alphaT, - # gammaT, - # beta_cov, - # alpha_cov, - # gamma_cov, - # self.t, - # self.kt - # ), - # float(t0[il]), - # float(t1[il]), - # ) - # integ = integ + np.log(q[il]) / 12 - # dint,_ = quad( - # lambda x: self._fzeroderiquanint( - # x, - # zqold, - # q[il], - # cov_locint, - # cov_scint, - # cov_shint, - # beta0, - # beta, - # alpha0, - # alpha, - # gamma0, - # gamma, - # betaT, - # alphaT, - # gammaT, - # beta_cov, - # alpha_cov, - # gamma_cov, - # self.t, - # self.kt - # ), - # float(t0[il]), - # float(t1[il]), - # ) - # zq = zq -integ / dint - # if np.abs(zq) > 1e-5: - # err = np.abs((zq - zqold) / zqold) - # else: - # err = np.abs(zq - zqold) - # iter1 += 1 - # if iter1 == 1000: - # zq = np.nan - # Warning("Maximum number of Newton iterations") - # if integ > 1e-2: - # zq = np.nan - # Warning("False zero, check it") - # zqout[il] = zq - - # return zqout - for il in range(m): # function of z whose root we want def F(z): From d9f160613cd49d9860d9a72cbae5cb6bd103a67e Mon Sep 17 00:00:00 2001 From: Colladito Date: Tue, 7 Oct 2025 10:48:14 +0200 Subject: [PATCH 08/18] [VCF] Improve stationary GEV fitted plot result --- bluemath_tk/distributions/_base_distributions.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/bluemath_tk/distributions/_base_distributions.py b/bluemath_tk/distributions/_base_distributions.py index d1c049a..469c2e9 100644 --- a/bluemath_tk/distributions/_base_distributions.py +++ b/bluemath_tk/distributions/_base_distributions.py @@ -555,13 +555,15 @@ def return_period(self, ax: plt.axes = None) -> Tuple[plt.figure, plt.axes]: else: fig = None + return_years = np.asarray([1.001, 1.1, 1.2, 1.3, 1.4, 1.5, 1.75, 2, 3, 4, 5, 7.5, 10, 15, 20, 25, 50, 100, 250, 500, 1000]) + ecdf_fitted = 1 - 1/return_years sorted_data = np.sort(self.data) exceedance_prob = 1 - self.ecdf return_period = 1 / exceedance_prob ax.plot( - return_period, - self.dist.qf(self.ecdf, *self.params), + return_years, + self.dist.qf(ecdf_fitted, *self.params), color="tab:red", label="Fitted Distribution", ) @@ -577,7 +579,7 @@ def return_period(self, ax: plt.axes = None) -> Tuple[plt.figure, plt.axes]: ax.set_xscale("log") ax.set_xticks([1, 2, 5, 10, 25, 50, 100, 250, 1000, 10000]) ax.set_xticklabels([1, 2, 5, 10, 25, 50, 100, 500, 1000, 10000]) - ax.set_xlim(right=np.max(return_period) * 1.2) + # ax.set_xlim(right=np.max(return_period) * 1.2) ax.set_xlabel("Return Period") ax.set_ylabel("Data Values") ax.set_title("Return Period Plot") From 8a9de9bbd4179f1ed84caf939cd0385309051389 Mon Sep 17 00:00:00 2001 From: Colladito Date: Tue, 7 Oct 2025 11:02:54 +0200 Subject: [PATCH 09/18] [VCF] Improve return period plot --- bluemath_tk/distributions/nonstat_gev.py | 50 ++++++++++++------------ 1 file changed, 26 insertions(+), 24 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index 5c8f856..3781983 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -4799,7 +4799,7 @@ def _quantile(self, prob=None, harm=False) -> np.ndarray: return Q - def plot(self, return_plot: bool = True, save: bool = False, init_year=2000): + def plot(self, return_plot: bool = False, save: bool = False, init_year=2000): """ Plot the location, scale and shape parameters, also the PP plot and QQ plot @@ -5425,7 +5425,7 @@ def plot(self, return_plot: bool = True, save: bool = False, init_year=2000): and self.nind_sc == 0 and self.nind_sh == 0 ) and return_plot: - self.ReturnPeriodPlot() + self.returnperiod_plot() def paramplot(self, save: bool=False): """ @@ -6347,7 +6347,7 @@ def _CDFGEVt(self): return F - def ReturnPeriodPlot(self, annualplot=True, save=False): + def returnperiod_plot(self, annualplot=True, monthly_plot=False, save=False): """ Funtion to plot the Aggregated Return period plot for each month and the annual Return period @@ -6366,18 +6366,19 @@ def ReturnPeriodPlot(self, annualplot=True, save=False): quanaggrA = np.zeros(nts) quanaggr = np.zeros((12, nts)) stdDqX = np.zeros((12, nts)) - for i in range(12): - for j in range(nts): - quanaggr[i, j] = self._aggquantile(1 - 1 / Ts[j], i / 12, (i + 1) / 12)[ - 0 - ] - # stdQuan = self._ConfidInterQuanAggregate( - # 1 - 1 / Ts[j], i / 12, (i + 1) / 12 - # ) - stdQuan = 0.1 - stdDqX[i, j] = stdQuan * norm.ppf( - 1 - (1 - self.quanval) / 2, loc=0, scale=1 - ) + if monthly_plot: + for i in range(12): + for j in range(nts): + quanaggr[i, j] = self._aggquantile(1 - 1 / Ts[j], i / 12, (i + 1) / 12)[ + 0 + ] + # stdQuan = self._ConfidInterQuanAggregate( + # 1 - 1 / Ts[j], i / 12, (i + 1) / 12 + # ) + stdQuan = 0.1 + stdDqX[i, j] = stdQuan * norm.ppf( + 1 - (1 - self.quanval) / 2, loc=0, scale=1 + ) # If annual data has to be plotted if annualplot: @@ -6428,15 +6429,16 @@ def ReturnPeriodPlot(self, annualplot=True, save=False): ] fig = plt.figure(figsize=(10, 6)) ax = fig.add_subplot() - for i in range(12): - ax.semilogx( - Ts, - quanaggr[i, :], - color=colors[i], - linestyle="-", - linewidth=1.2, - label=labels[i], - ) + if monthly_plot: + for i in range(12): + ax.semilogx( + Ts, + quanaggr[i, :], + color=colors[i], + linestyle="-", + linewidth=1.2, + label=labels[i], + ) # Anual return periods if annualplot: From 49d360add9dfb163499462a6020b6bc59d536038 Mon Sep 17 00:00:00 2001 From: Colladito Date: Tue, 7 Oct 2025 13:37:13 +0200 Subject: [PATCH 10/18] [VCF] Update returnperiod_plot add optional parameter to compute annual confidence bands --- bluemath_tk/distributions/nonstat_gev.py | 77 ++++++++++++++---------- 1 file changed, 46 insertions(+), 31 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index 3781983..4f14692 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -6347,7 +6347,7 @@ def _CDFGEVt(self): return F - def returnperiod_plot(self, annualplot=True, monthly_plot=False, save=False): + def returnperiod_plot(self, annualplot=True, monthly_plot=False, save=False, conf_int=False): """ Funtion to plot the Aggregated Return period plot for each month and the annual Return period @@ -6372,30 +6372,32 @@ def returnperiod_plot(self, annualplot=True, monthly_plot=False, save=False): quanaggr[i, j] = self._aggquantile(1 - 1 / Ts[j], i / 12, (i + 1) / 12)[ 0 ] + # DO NOT COMPUTE CONFIDENCE INTERVAL IN MONTHLY RETURN PERIODS # stdQuan = self._ConfidInterQuanAggregate( # 1 - 1 / Ts[j], i / 12, (i + 1) / 12 # ) - stdQuan = 0.1 - stdDqX[i, j] = stdQuan * norm.ppf( - 1 - (1 - self.quanval) / 2, loc=0, scale=1 - ) + # # stdQuan = 0.1 + # stdDqX[i, j] = stdQuan * norm.ppf( + # 1 - (1 - self.quanval) / 2, loc=0, scale=1 + # ) # If annual data has to be plotted if annualplot: for j in range(nts): quanaggrA[j] = self._aggquantile(1 - 1 / Ts[j], 0, 1)[0] # Confidence intervals - stdup = np.zeros(nts) - stdlo = np.zeros(nts) - for i in range(nts): - # stdQuan = self._ConfidInterQuanAggregate(1 - 1 / Ts[i], 0, 1) - stdQuan = 0.1 - stdup[i] = quanaggrA[i] + stdQuan * norm.ppf( - 1 - (1 - self.quanval) / 2, loc=0, scale=1 - ) - stdlo[i] = quanaggrA[i] - stdQuan * norm.ppf( - 1 - (1 - self.quanval) / 2, loc=0, scale=1 - ) + if conf_int: + stdup = np.zeros(nts) + stdlo = np.zeros(nts) + for i in range(nts): + stdQuan = self._ConfidInterQuanAggregate(1 - 1 / Ts[i], 0, 1) + # stdQuan = 0.1 + stdup[i] = quanaggrA[i] + stdQuan * norm.ppf( + 1 - (1 - self.quanval) / 2, loc=0, scale=1 + ) + stdlo[i] = quanaggrA[i] - stdQuan * norm.ppf( + 1 - (1 - self.quanval) / 2, loc=0, scale=1 + ) ## Plot the return periods # datemax_mod = self.t % 1 @@ -6454,17 +6456,18 @@ def returnperiod_plot(self, annualplot=True, monthly_plot=False, save=False): Tapprox = 1 / (1 - ProHsmaxsort) # idx = np.where(Tapprox >= 2)[0] # ax.semilogx(Tapprox[idx], hmaxsort[idx], "ok", markersize=1.6) - ax.semilogx(Tapprox, hmaxsort, "ok", markersize=2) - ax.semilogx(Ts, stdlo, "--k", linewidth=1.1) - ax.semilogx(Ts, stdup, "--k", linewidth=1.1) - + ax.semilogx(Tapprox, hmaxsort, "+", markersize=7, label="Data", color="black") + if conf_int: + ax.semilogx(Ts, stdlo, linewidth=1.1, color = "gray", linestyle='dashed') + ax.semilogx(Ts, stdup, linewidth=1.1, color="gray", linestyle='dashed') + ax.set_xlabel("Return Period (Years)") ax.set_ylabel(f"{self.var_name}") ax.legend(loc="center left", bbox_to_anchor=(1, 0.5)) ax.set_xticks([1, 2, 5, 10, 20, 50, 100, 250, 500]) ax.set_xticklabels([1, 2, 5, 10, 20, 50, 100, 250, 500]) ax.set_xlim(left=0.9, right=Ts[-1] + 50) - ax.set_ylim(bottom=0) + # ax.set_ylim(bottom=0) # ax.set_title(f"Aggregate Quantiles ({self.var_name})") ax.grid(True) plt.tight_layout() @@ -6489,6 +6492,10 @@ def _aggquantile( beta_cov=None, alpha_cov=None, gamma_cov=None, + list_loc=None, + list_sc=None, + list_sh=None, + covariates=None ) -> np.ndarray: """ Function to compute the aggregated quantile between two time stamps @@ -6553,6 +6560,14 @@ def _aggquantile( alpha_cov = self.alpha_cov if gamma_cov is None: gamma_cov = self.gamma_cov + if list_loc is None: + list_loc = self.list_loc + if list_sc is None: + list_sc = self.list_sc + if list_sh is None: + list_sh = self.list_sh + if covariates is None: + covariates = self.covariates # Deal with scalars and vectors @@ -6584,15 +6599,15 @@ def _aggquantile( if pos.size: for i in range(beta_cov.size): cov_locint[i] = np.mean( - self.covariates.iloc[pos, self.list_loc[i]].values + covariates.iloc[pos, list_loc[i]].values ) for i in range(alpha_cov.size): cov_scint[i] = np.mean( - self.covariates.iloc[pos, self.list_sc[i]].values + covariates.iloc[pos, list_sc[i]].values ) for i in range(gamma_cov.size): cov_shint[i] = np.mean( - self.covariates.iloc[pos, self.list_sh[i]].values + covariates.iloc[pos, list_sh[i]].values ) else: cov_locint = np.zeros(beta_cov.size) @@ -6608,7 +6623,7 @@ def _aggquantile( beta, betaT, beta_cov, - self.covariates.iloc[:, self.list_loc].values, + covariates.iloc[:, list_loc].values, cov_locint, self.t, x, @@ -6934,8 +6949,8 @@ def _ConfidInterQuanAggregate(self, q, t0, t1) -> np.ndarray: beta_covlb[i] = self.beta_cov[i] * (1 + epsi) beta_covub[i] = self.beta_cov[i] * (1 - epsi) jacob[aux] = ( - self._aggquantile(q, t0, t1, beta_cov=beta_covlb[i])[0] - - self._aggquantile(q, t0, t1, beta_cov=beta_covlb[i])[0] + self._aggquantile(q, t0, t1, beta_cov=np.atleast_1d(beta_covlb[i]), list_loc=np.atleast_1d(self.list_loc[i]))[0] + - self._aggquantile(q, t0, t1, beta_cov=np.atleast_1d(beta_covlb[i]), list_loc=np.atleast_1d(self.list_loc[i]))[0] ) / (2 * self.beta_cov[i] * epsi) else: jacob[aux] = 0 @@ -6978,8 +6993,8 @@ def _ConfidInterQuanAggregate(self, q, t0, t1) -> np.ndarray: alpha_covlb[i] = self.alpha_cov[i] * (1 + epsi) alpha_covub[i] = self.alpha_cov[i] * (1 - epsi) jacob[aux] = ( - self._aggquantile(q, t0, t1, alpha_cov=alpha_covlb[i])[0] - - self._aggquantile(q, t0, t1, alpha_cov=alpha_covlb[i])[0] + self._aggquantile(q, t0, t1, alpha_cov=np.atleast_1d(alpha_covlb[i]), list_sc=np.atleast_1d(self.list_sc[i]))[0] + - self._aggquantile(q, t0, t1, alpha_cov=np.atleast_1d(alpha_covlb[i]), list_sc=np.atleast_1d(self.list_sc[i]))[0] ) / (2 * self.alpha_cov[i] * epsi) else: jacob[aux] = 0 @@ -7021,8 +7036,8 @@ def _ConfidInterQuanAggregate(self, q, t0, t1) -> np.ndarray: gamma_covlb[i] = self.gamma_cov[i] * (1 + epsi) gamma_covub[i] = self.gamma_cov[i] * (1 - epsi) jacob[aux] = ( - self._aggquantile(q, t0, t1, gamma_cov=gamma_covlb[i])[0] - - self._aggquantile(q, t0, t1, gamma_cov=gamma_covub[i])[0] + self._aggquantile(q, t0, t1, gamma_cov=gamma_covlb[i], list_sh=np.atleast_1d(self.list_sh[i]))[0] + - self._aggquantile(q, t0, t1, gamma_cov=gamma_covub[i], list_sh=np.atleast_1d(self.list_sh[i]))[0] ) / (2 * self.gamma_cov[i] * epsi) else: jacob[aux] = 0 From 2d00c08849b47031824a1c96f01dc03d43d7cb6a Mon Sep 17 00:00:00 2001 From: Colladito Date: Tue, 7 Oct 2025 17:15:37 +0200 Subject: [PATCH 11/18] [VCF] Add logging info and conf_int parameter in returnperiod_plot() --- bluemath_tk/distributions/nonstat_gev.py | 509 +++++++++++++++-------- 1 file changed, 326 insertions(+), 183 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index 4f14692..1f97158 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -71,6 +71,11 @@ def __init__( """ Initiliaze the Non-Stationary GEV. """ + debug=1 + self.set_logger_name( + name=self.__class__.__name__, level="DEBUG" if debug else "INFO" + ) + self.logger.debug("Initializing NonStatGEV") # Initialize arguments self.xt = xt @@ -1203,6 +1208,7 @@ def _fit( list_sh: list = [], ntrend_sh: int = 0, xini: Optional[np.ndarray] = None, + options: dict = None, ) -> dict: """ Auxiliar function to determine the optimal parameters of given Non-Stationary GEV @@ -1229,7 +1235,22 @@ def _fit( If trends in shape are included. xini : Optional[np.ndarray], default = None Initial parameter if given. + options : dict, default={ + "gtol": 1e-5, + "xtol": 1e-5, + "barrier_tol": 1e-4, + "maxiter": 200, + } + Dictionary with the optimization options, see scipy.minimize method "trust-constr" options """ + # Fitting options + if options is None: + options = { + "gtol": 1e-5, + "xtol": 1e-5, + "barrier_tol": 1e-4, + "maxiter": 500, + } # Total number of parameters to be estimated nmu = 2 * nmu @@ -1442,12 +1463,7 @@ def _fit( ntrend_sh, list_sh, ), - options={ - "gtol": 1e-5, - "xtol": 1e-5, - "barrier_tol": 1e-4, - "maxiter": 200, # Lower max iterations - }, # Options + options=options, # Options method="trust-constr", ) @@ -1492,12 +1508,7 @@ def _fit( ntrend_sh, list_sh, ), - options={ - "gtol": 1e-5, - "xtol": 1e-5, - "barrier_tol": 1e-4, - "maxiter": 200, # Lower max iterations - }, # Options + options=options, # Options method="trust-constr", ) @@ -2823,6 +2834,7 @@ def fit( list_sc: Optional[Union[list, str]] = "all", ntrend_sh: int = 1, list_sh: Optional[Union[list, str]] = "all", + options: dict = None, plot: bool = False, ) -> dict: """ @@ -2868,6 +2880,8 @@ def fit( - hessian: Hessian matrix of the log-likelihood function at the optimal solution - invI0: Inverse of Fisher information matrix """ + self.logger.debug("Fixed fit (.fit())") + if nmu % 2 != 0 or npsi % 2 != 0 or ngamma % 2 != 0: raise ValueError("Check the input of harmonics, must be even!") @@ -2908,6 +2922,7 @@ def fit( ntrend_sc=ntrend_sc, list_sh=list_sh, ntrend_sh=ntrend_sh, + options=options, ) # Update parameters in the class @@ -4594,7 +4609,9 @@ def _parametro( y = y + beta_cov[i] * indicesint[i] else: for i in range(nind): - indicesintaux = self._search(times, covariates[:, i], np.asarray(x).flatten()) + indicesintaux = self._search( + times, covariates[:, i], np.asarray(x).flatten() + ) y = y + beta_cov[i] * indicesintaux else: for i in range(nind): @@ -4799,7 +4816,7 @@ def _quantile(self, prob=None, harm=False) -> np.ndarray: return Q - def plot(self, return_plot: bool = False, save: bool = False, init_year=2000): + def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0): """ Plot the location, scale and shape parameters, also the PP plot and QQ plot @@ -4811,6 +4828,8 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year=2000): If True, return period plot is plotted save : bool, default=False If True, save all the figures in a "Figures/" + init_year : int, default=0 + Initial year for plotting purposes """ # Parameter Evaluation @@ -4948,8 +4967,8 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year=2000): # Location and Scale parameter plotting t_anual = np.mod(self.t, 1) quan95 = self._quantile() - rp10 = self._quantile(1-1/10) - rp100 = self._quantile(1-1/100) + rp10 = self._quantile(1 - 1 / 10) + rp100 = self._quantile(1 - 1 / 100) if ( self.betaT is None @@ -5018,7 +5037,7 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year=2000): label="Data", ) ax1.plot( - self.t+init_year, + self.t + init_year, mut, label="Location", linewidth=2, @@ -5027,67 +5046,57 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year=2000): alpha=1, ) ax1.fill_between( - self.t+init_year, + self.t + init_year, mut - psit, mut + psit, label=r"Location $\pm$ Scale", color="tab:blue", - alpha=0.4 + alpha=0.4, ) ax1.plot( - self.t+init_year, + self.t + init_year, rp10, label="10 years", linewidth=1, color="tab:red", linestyle="--", - alpha=.9, + alpha=0.9, ) ax1.plot( - self.t+init_year, + self.t + init_year, rp100, label="100 years", linewidth=1, color="tab:green", linestyle="--", - alpha=.9, + alpha=0.9, ) ax1.plot( - self.t+init_year, + self.t + init_year, epst, label="Shape", linewidth=1, color="tab:orange", - alpha=.9, + alpha=0.9, ) # TODO: Add aggregated return period lines - # rt_10 = np.zeros(40) - # for year in range(40): + # n_years = int(np.ceil(self.t[-1])) + # rt_10 = np.zeros(n_years) + # for year in range(n_years): # rt_10[year] = self._aggquantile(1-1/10, year, year+1) # 10-year return level at each year - # rt_20 = np.zeros(40) - # for year in range(40): - # rt_20[year] = self._aggquantile(1-1/20, year, year+1) - # rt_50 = np.zeros(40) - # for year in range(40): - # rt_50[year] = self._aggquantile(1-1/50, year, year+1) - # rt_100 = np.zeros(40) - # for year in range(40): + # rt_100 = np.zeros(n_years) + # for year in range(n_years): # rt_100[year] = self._aggquantile(1-1/100, year, year+1) - # l4 = ax1.plot( - # np.arange(0,40), rt_10, linestyle="-", linewidth=1, label="10 years", - # ) - # l5 = ax1.plot( - # np.arange(0,40), rt_20, linestyle="-", linewidth=1, label="20 years", - # ) - # l6 = ax1.plot( - # np.arange(0,40), rt_50, linestyle="-", linewidth=1, label="50 years", + + # ax1.plot( + # np.arange(init_year,init_year+n_years), rt_10, linestyle="-", linewidth=1, label="10 years", color="tab:red" # ) - # l7 = ax1.plot( - # np.arange(0,40), rt_100, linestyle="-", linewidth=1, label="100 years", + # ax1.plot( + # np.arange(init_year, init_year+n_years), rt_100, linestyle="-", linewidth=1, label="100 years", color="tab:green" # ) ax1.set_xlabel("Time (years)") @@ -5200,7 +5209,7 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year=2000): uppermaxs, color=self.colors[0], alpha=0.3, - label="Location +- scale" + label="Location +- scale", ) # l2 = ax2.plot(self.t[mask_month], psit[mask_month], label=r'$\psi_t$', linewidth=2, color=self.colors[1], alpha=1) # ax2.fill_between(self.t[mask_month], ci_low_psit[mask_month], ci_up_psit[mask_month], color=self.colors[1], @@ -5226,7 +5235,9 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year=2000): ax1.legend(handles=handles, loc="best") ax1.margins(x=0.01) if save: - plt.savefig(f"Figures/Evolution_Location_FirstYear{self.var_name}.png", dpi=300) + plt.savefig( + f"Figures/Evolution_Location_FirstYear{self.var_name}.png", dpi=300 + ) plt.show() ### Shape parameter plot @@ -5238,7 +5249,7 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year=2000): # epst - norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdepst # ) - # LOCATION + # LOCATION plt.figure(figsize=(20, 6)) plt.plot(self.t + init_year, mut, color="tab:blue") # plt.fill_between( @@ -5257,7 +5268,7 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year=2000): plt.savefig(f"Figures/Evolution_Location{self.var_name}.png", dpi=300) plt.show() - # SCALE + # SCALE plt.figure(figsize=(20, 6)) plt.plot(self.t + init_year, psit, color="tab:green") # plt.fill_between( @@ -5294,8 +5305,6 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year=2000): if save: plt.savefig(f"Figures/Evolution_Shape{self.var_name}.png", dpi=300) plt.show() - - # ### Harmonic Location parameter plot # if self.nmu > 0: @@ -5427,11 +5436,15 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year=2000): ) and return_plot: self.returnperiod_plot() - def paramplot(self, save: bool=False): + def paramplot(self, save: bool = False): """ Create a heatmap of parameter sensitivities for location, scale and shape. + + Parameters + ---------- + save : bool, False """ - # Get parameters + # Get parameters loc_params = [] param_names = [] if self.fit_result["beta0"] is not None: @@ -5440,7 +5453,7 @@ def paramplot(self, save: bool=False): if self.fit_result["beta"].size > 0: for i, b in enumerate(self.fit_result["beta"]): loc_params.append(b) - param_names.append(f"PC{i}") + param_names.append(f"Harm {i}") if self.fit_result["beta_cov"].size > 0: for i, b in enumerate(self.fit_result["beta_cov"]): loc_params.append(b) @@ -5448,7 +5461,7 @@ def paramplot(self, save: bool=False): if self.fit_result["betaT"].size > 0: loc_params.append(self.fit_result["betaT"]) param_names.append("Trend") - + # Scale parameters scale_params = [] if self.fit_result["alpha0"] is not None: @@ -5459,8 +5472,8 @@ def paramplot(self, save: bool=False): scale_params.extend(self.fit_result["alpha_cov"]) if self.fit_result["alphaT"].size > 0: scale_params.append(self.fit_result["alphaT"]) - - # Shape parameters + + # Shape parameters shape_params = [] if self.fit_result["gamma0"] is not None: shape_params.append(self.fit_result["gamma0"]) @@ -5474,44 +5487,49 @@ def paramplot(self, save: bool=False): # Create data matrix max_len = max(len(param_names), len(scale_params), len(shape_params)) data = np.zeros((3, len(param_names))) - data[0,:len(loc_params)] = loc_params - data[1,:len(scale_params)] = scale_params - data[2,:len(shape_params)] = shape_params + data[0, : len(loc_params)] = loc_params + data[1, : len(scale_params)] = scale_params + data[2, : len(shape_params)] = shape_params # Create figure - fig, ax = plt.subplots(figsize=(20,4)) - + fig, ax = plt.subplots(figsize=(20, 4)) + # Create heatmap - im = ax.imshow(data, cmap='RdBu', aspect='auto', vmin=-2, vmax=2) - + im = ax.imshow(data, cmap="RdBu", aspect="auto", vmin=-2, vmax=2) + # Add colorbar cbar = plt.colorbar(im) - cbar.set_label('Coefficient value') - + cbar.set_label("Coefficient value") + # Add text annotations for i in range(data.shape[0]): for j in range(data.shape[1]): - text = ax.text(j, i, f'{data[i, j]:.2f}', - ha="center", va="center", color="black") - + text = ax.text( + j, i, f"{data[i, j]:.2f}", ha="center", va="center", color="black" + ) + # Configure axes ax.set_xticks(np.arange(len(param_names))) ax.set_yticks(np.arange(3)) ax.set_xticklabels(param_names) - ax.set_yticklabels(['Location', 'Scale', 'Shape']) - + ax.set_yticklabels(["Location", "Scale", "Shape"]) + # Rotate x labels for better readability plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor") - + # Add title and adjust layout - ax.set_title('Time-dependent GEV Parameters') - ax.set_xlabel('Components') - + ax.set_title("Time-dependent GEV Parameters") + ax.set_xlabel("Components") + plt.tight_layout() - + if save: - plt.savefig(f"Figures/Parameters_Heatmap_{self.var_name}.png", dpi=300, bbox_inches='tight') - + plt.savefig( + f"Figures/Parameters_Heatmap_{self.var_name}.png", + dpi=300, + bbox_inches="tight", + ) + plt.show() def QQplot(self, save: bool = False): @@ -6347,7 +6365,13 @@ def _CDFGEVt(self): return F - def returnperiod_plot(self, annualplot=True, monthly_plot=False, save=False, conf_int=False): + def returnperiod_plot( + self, + annualplot: bool = True, + conf_int: bool = False, + monthly_plot: bool = False, + save: bool = False, + ): """ Funtion to plot the Aggregated Return period plot for each month and the annual Return period @@ -6355,9 +6379,39 @@ def returnperiod_plot(self, annualplot=True, monthly_plot=False, save=False, con ---------- annualplot : bool, default=True Whether to plot the annual return period plot + conf_int : bool, default=False + Whether to plot the confidence bands for annual return periods + Heavy computational time + monthly_plot : bool, default=False + Wheter to plot the return periods grouped by months + save : bool, default=False + Whether to save the plot """ - - Ts = np.array([1.1, 1.5, 2, 3, 4, 5, 7.5, 10, 15, 20, 30, 40, 50, 75, 100, 150, 200, 500, 1000]) + self.logger.debug("Calling: Return Period plot (.returnperiod_plot())") + + Ts = np.array( + [ + 1.1, + 1.5, + 2, + 3, + 4, + 5, + 7.5, + 10, + 15, + 20, + 30, + 40, + 50, + 75, + 100, + 150, + 200, + 500, + 1000, + ] + ) # Ts = np.concatenate( # (np.arange(2, 10, 1), np.arange(10, 100, 10), np.arange(100, 501, 100)) # ) @@ -6369,9 +6423,9 @@ def returnperiod_plot(self, annualplot=True, monthly_plot=False, save=False, con if monthly_plot: for i in range(12): for j in range(nts): - quanaggr[i, j] = self._aggquantile(1 - 1 / Ts[j], i / 12, (i + 1) / 12)[ - 0 - ] + quanaggr[i, j] = self._aggquantile( + 1 - 1 / Ts[j], i / 12, (i + 1) / 12 + )[0] # DO NOT COMPUTE CONFIDENCE INTERVAL IN MONTHLY RETURN PERIODS # stdQuan = self._ConfidInterQuanAggregate( # 1 - 1 / Ts[j], i / 12, (i + 1) / 12 @@ -6383,8 +6437,10 @@ def returnperiod_plot(self, annualplot=True, monthly_plot=False, save=False, con # If annual data has to be plotted if annualplot: + self.logger.debug("ReturnPeriodPlot: Annual return period.") for j in range(nts): quanaggrA[j] = self._aggquantile(1 - 1 / Ts[j], 0, 1)[0] + self.logger.debug(f"ReturnPeriodPlot: Annual return period {j} finished.") # Confidence intervals if conf_int: stdup = np.zeros(nts) @@ -6456,11 +6512,13 @@ def returnperiod_plot(self, annualplot=True, monthly_plot=False, save=False, con Tapprox = 1 / (1 - ProHsmaxsort) # idx = np.where(Tapprox >= 2)[0] # ax.semilogx(Tapprox[idx], hmaxsort[idx], "ok", markersize=1.6) - ax.semilogx(Tapprox, hmaxsort, "+", markersize=7, label="Data", color="black") + ax.semilogx( + Tapprox, hmaxsort, "+", markersize=7, label="Data", color="black" + ) if conf_int: - ax.semilogx(Ts, stdlo, linewidth=1.1, color = "gray", linestyle='dashed') - ax.semilogx(Ts, stdup, linewidth=1.1, color="gray", linestyle='dashed') - + ax.semilogx(Ts, stdlo, linewidth=1.1, color="gray", linestyle="dashed") + ax.semilogx(Ts, stdup, linewidth=1.1, color="gray", linestyle="dashed") + ax.set_xlabel("Return Period (Years)") ax.set_ylabel(f"{self.var_name}") ax.legend(loc="center left", bbox_to_anchor=(1, 0.5)) @@ -6495,7 +6553,7 @@ def _aggquantile( list_loc=None, list_sc=None, list_sh=None, - covariates=None + covariates=None, ) -> np.ndarray: """ Function to compute the aggregated quantile between two time stamps @@ -6569,9 +6627,8 @@ def _aggquantile( if covariates is None: covariates = self.covariates - # Deal with scalars and vectors - beta_cov = np.atleast_1d(beta_cov) + beta_cov = np.atleast_1d(beta_cov) alpha_cov = np.atleast_1d(alpha_cov) gamma_cov = np.atleast_1d(gamma_cov) @@ -6598,17 +6655,11 @@ def _aggquantile( cov_shint = np.zeros_like(gamma_cov) if pos.size: for i in range(beta_cov.size): - cov_locint[i] = np.mean( - covariates.iloc[pos, list_loc[i]].values - ) + cov_locint[i] = np.mean(covariates.iloc[pos, list_loc[i]].values) for i in range(alpha_cov.size): - cov_scint[i] = np.mean( - covariates.iloc[pos, list_sc[i]].values - ) + cov_scint[i] = np.mean(covariates.iloc[pos, list_sc[i]].values) for i in range(gamma_cov.size): - cov_shint[i] = np.mean( - covariates.iloc[pos, list_sh[i]].values - ) + cov_shint[i] = np.mean(covariates.iloc[pos, list_sh[i]].values) else: cov_locint = np.zeros(beta_cov.size) cov_scint = np.zeros(alpha_cov.size) @@ -6617,7 +6668,7 @@ def _aggquantile( # Require quantile zqout = np.zeros(m) - media,_ = quad( + media, _ = quad( lambda x: self._parametro( beta0, beta, @@ -6653,13 +6704,29 @@ def _aggquantile( def F(z): integ, _ = quad( lambda x: self._fzeroquanint( - x, z, q[il], - cov_locint, cov_scint, cov_shint, - beta0, beta, alpha0, alpha, - gamma0, gamma, betaT, alphaT, gammaT, - beta_cov, alpha_cov, gamma_cov, - self.t, self.kt), - float(t0[il]), float(t1[il]) + x, + z, + q[il], + cov_locint, + cov_scint, + cov_shint, + beta0, + beta, + alpha0, + alpha, + gamma0, + gamma, + betaT, + alphaT, + gammaT, + beta_cov, + alpha_cov, + gamma_cov, + self.t, + self.kt, + ), + float(t0[il]), + float(t1[il]), ) return integ + np.log(q[il]) / 12.0 @@ -6667,16 +6734,18 @@ def F(z): # root finding: start near `media`, bracket if possible sol = root_scalar( F, - x0=media, x1=media + 1.0, # secant starting points + x0=media, + x1=media + 1.0, # secant starting points method="secant", - xtol=1e-6, rtol=1e-6, - maxiter=200 + xtol=1e-6, + rtol=1e-6, + maxiter=200, ) if sol.converged: if abs(F(sol.root)) < 1e-2: zqout[il] = sol.root else: - zqout[il] = np.nan # "False zero" check + zqout[il] = np.nan # "False zero" check else: zqout[il] = np.nan except Exception: @@ -6705,7 +6774,7 @@ def _fzeroquanint( alpha_cov, gamma_cov, times=None, - ktold=None + ktold=None, ) -> np.ndarray: """ Auxiliar function to solve the quantile @@ -6715,7 +6784,7 @@ def _fzeroquanint( zn : np.ndarray """ if ktold is None: - ktold=self.kt + ktold = self.kt # Evaluate the location parameter at each time t as a function of the actual values of the parameters given by p mut1 = self._parametro( @@ -6761,7 +6830,9 @@ def _fzeroquanint( epst[posG] = 1 if times is not None: - kt2 = np.interp(t, np.asarray(self.t, float), np.asarray(ktold, float)).flatten() + kt2 = np.interp( + t, np.asarray(self.t, float), np.asarray(ktold, float) + ).flatten() else: kt2 = np.ones_like(mut1) @@ -6860,7 +6931,9 @@ def _fzeroderiquanint( epst[posG] = 1 if times is not None: - kt2 = np.interp(t, np.asarray(self.t, float), np.asarray(ktold, float)).flatten() + kt2 = np.interp( + t, np.asarray(self.t, float), np.asarray(ktold, float) + ).flatten() else: kt2 = np.ones_like(mut1) @@ -6949,8 +7022,20 @@ def _ConfidInterQuanAggregate(self, q, t0, t1) -> np.ndarray: beta_covlb[i] = self.beta_cov[i] * (1 + epsi) beta_covub[i] = self.beta_cov[i] * (1 - epsi) jacob[aux] = ( - self._aggquantile(q, t0, t1, beta_cov=np.atleast_1d(beta_covlb[i]), list_loc=np.atleast_1d(self.list_loc[i]))[0] - - self._aggquantile(q, t0, t1, beta_cov=np.atleast_1d(beta_covlb[i]), list_loc=np.atleast_1d(self.list_loc[i]))[0] + self._aggquantile( + q, + t0, + t1, + beta_cov=np.atleast_1d(beta_covlb[i]), + list_loc=np.atleast_1d(self.list_loc[i]), + )[0] + - self._aggquantile( + q, + t0, + t1, + beta_cov=np.atleast_1d(beta_covlb[i]), + list_loc=np.atleast_1d(self.list_loc[i]), + )[0] ) / (2 * self.beta_cov[i] * epsi) else: jacob[aux] = 0 @@ -6993,8 +7078,20 @@ def _ConfidInterQuanAggregate(self, q, t0, t1) -> np.ndarray: alpha_covlb[i] = self.alpha_cov[i] * (1 + epsi) alpha_covub[i] = self.alpha_cov[i] * (1 - epsi) jacob[aux] = ( - self._aggquantile(q, t0, t1, alpha_cov=np.atleast_1d(alpha_covlb[i]), list_sc=np.atleast_1d(self.list_sc[i]))[0] - - self._aggquantile(q, t0, t1, alpha_cov=np.atleast_1d(alpha_covlb[i]), list_sc=np.atleast_1d(self.list_sc[i]))[0] + self._aggquantile( + q, + t0, + t1, + alpha_cov=np.atleast_1d(alpha_covlb[i]), + list_sc=np.atleast_1d(self.list_sc[i]), + )[0] + - self._aggquantile( + q, + t0, + t1, + alpha_cov=np.atleast_1d(alpha_covlb[i]), + list_sc=np.atleast_1d(self.list_sc[i]), + )[0] ) / (2 * self.alpha_cov[i] * epsi) else: jacob[aux] = 0 @@ -7036,8 +7133,20 @@ def _ConfidInterQuanAggregate(self, q, t0, t1) -> np.ndarray: gamma_covlb[i] = self.gamma_cov[i] * (1 + epsi) gamma_covub[i] = self.gamma_cov[i] * (1 - epsi) jacob[aux] = ( - self._aggquantile(q, t0, t1, gamma_cov=gamma_covlb[i], list_sh=np.atleast_1d(self.list_sh[i]))[0] - - self._aggquantile(q, t0, t1, gamma_cov=gamma_covub[i], list_sh=np.atleast_1d(self.list_sh[i]))[0] + self._aggquantile( + q, + t0, + t1, + gamma_cov=gamma_covlb[i], + list_sh=np.atleast_1d(self.list_sh[i]), + )[0] + - self._aggquantile( + q, + t0, + t1, + gamma_cov=gamma_covub[i], + list_sh=np.atleast_1d(self.list_sh[i]), + )[0] ) / (2 * self.gamma_cov[i] * epsi) else: jacob[aux] = 0 @@ -7046,7 +7155,7 @@ def _ConfidInterQuanAggregate(self, q, t0, t1) -> np.ndarray: stdQuan = np.sqrt(jacob.T @ self.invI0 @ jacob) return stdQuan - + def summary(self): """ Print a summary of the fitted model, including parameter estimates, standard errors and fit statistics. @@ -7056,133 +7165,169 @@ def summary(self): z_norm = norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) print(f"\nFitted Time-Dependent GEV model for {self.var_name}") - print("="*70) - + print("=" * 70) + # Header format param_header = "Parameter".ljust(15) estimate_header = "Estimate".rjust(12) se_header = "Std. Error".rjust(15) - ci_header = f"{self.quanval*100:.0f}% CI".center(25) + ci_header = f"{self.quanval * 100:.0f}% CI".center(25) header = f"{param_header} {estimate_header} {se_header} {ci_header}" print("\nLocation Parameters") - print("-"*70) + print("-" * 70) print(header) - print("-"*70) - + print("-" * 70) + # Parameter line format def format_line(name, value, std_err): - ci_low = value - std_err*z_norm - ci_up = value + std_err*z_norm + ci_low = value - std_err * z_norm + ci_up = value + std_err * z_norm return f"{name:<15} {value:>12.4f} {std_err:>15.4f} {f'[{ci_low:>8.4f},{ci_up:>8.4f}]':>25}" - + # Location parameters print(format_line("Beta0", self.beta0, std_params[param_idx])) param_idx += 1 for i in range(self.nmu): - print(format_line(f"Beta{i+1} (sin)", self.beta[2*i], std_params[param_idx])) + print( + format_line( + f"Beta{i + 1} (sin)", self.beta[2 * i], std_params[param_idx] + ) + ) param_idx += 1 - print(format_line(f"Beta{i+1} (cos)", self.beta[2*i+1], std_params[param_idx])) + print( + format_line( + f"Beta{i + 1} (cos)", self.beta[2 * i + 1], std_params[param_idx] + ) + ) param_idx += 1 - + if self.ntrend_loc > 0: print(format_line("BetaT", self.betaT, std_params[param_idx])) param_idx += 1 - + for i in range(self.nind_loc): - print(format_line(f"Beta_cov{i+1}", self.beta_cov[i], std_params[param_idx])) + print( + format_line(f"Beta_cov{i + 1}", self.beta_cov[i], std_params[param_idx]) + ) param_idx += 1 - print("\nScale Parameters") - print("-"*70) + print("\nScale Parameters") + print("-" * 70) print(header) - print("-"*70) - + print("-" * 70) + print(format_line("Alpha0", self.alpha0, std_params[param_idx])) param_idx += 1 - + for i in range(self.npsi): - print(format_line(f"Alpha{i+1} (sin)", self.alpha[2*i], std_params[param_idx])) + print( + format_line( + f"Alpha{i + 1} (sin)", self.alpha[2 * i], std_params[param_idx] + ) + ) param_idx += 1 - print(format_line(f"Alpha{i+1} (cos)", self.alpha[2*i+1], std_params[param_idx])) + print( + format_line( + f"Alpha{i + 1} (cos)", self.alpha[2 * i + 1], std_params[param_idx] + ) + ) param_idx += 1 - + if self.ntrend_sc > 0: print(format_line("AlphaT", self.alphaT, std_params[param_idx])) param_idx += 1 - + for i in range(self.nind_sc): - print(format_line(f"Alpha_cov{i+1}", self.alpha_cov[i], std_params[param_idx])) + print( + format_line( + f"Alpha_cov{i + 1}", self.alpha_cov[i], std_params[param_idx] + ) + ) param_idx += 1 print("\nShape Parameters") - print("-"*70) + print("-" * 70) print(header) - print("-"*70) - + print("-" * 70) + if self.ngamma0 > 0: print(format_line("Gamma0", self.gamma0, std_params[param_idx])) param_idx += 1 - + for i in range(self.ngamma): - print(format_line(f"Gamma{i+1} (sin)", self.gamma[2*i], std_params[param_idx])) + print( + format_line( + f"Gamma{i + 1} (sin)", self.gamma[2 * i], std_params[param_idx] + ) + ) param_idx += 1 - print(format_line(f"Gamma{i+1} (cos)", self.gamma[2*i+1], std_params[param_idx])) + print( + format_line( + f"Gamma{i + 1} (cos)", self.gamma[2 * i + 1], std_params[param_idx] + ) + ) param_idx += 1 - + if self.ntrend_sh > 0: print(format_line("GammaT", self.gammaT, std_params[param_idx])) param_idx += 1 - + for i in range(self.nind_sh): - print(format_line(f"Gamma_cov{i+1}", self.gamma_cov[i], std_params[param_idx])) + print( + format_line( + f"Gamma_cov{i + 1}", self.gamma_cov[i], std_params[param_idx] + ) + ) param_idx += 1 - + print("\nFit Statistics") - print("-"*70) + print("-" * 70) stats_width = 30 - print(f"{'Log-likelihood:':<{stats_width}} {-self.fit_result['negloglikelihood']:>.4f}") + print( + f"{'Log-likelihood:':<{stats_width}} {-self.fit_result['negloglikelihood']:>.4f}" + ) print(f"{'AIC:':<{stats_width}} {self.fit_result['AIC']:>.4f}") print(f"{'Number of parameters:':<{stats_width}} {self.fit_result['n_params']}") - - + print(f"Success: {self.fit_result['success']}") + print(f"Message: {self.fit_result['message']}") def bsimp(fun, a, b, n=None, epsilon=1e-8, trace=0): """ BSIMP Numerically evaluate integral, low order method. - I = BSIMP('F',A,B) approximates the integral of F(X) from A to B + I = BSIMP('F',A,B) approximates the integral of F(X) from A to B within a relative error of 1e-3 using an iterative Simpson's rule. 'F' is a string containing the name of the function. Function F must return a vector of output values if given a vector of input values.% I = BSIMP('F',A,B,EPSILON) integrates to a total error of EPSILON. % - I = BSIMP('F',A,B,N,EPSILON,TRACE,TRACETOL) integrates to a - relative error of EPSILON, - beginning with n subdivisions of the interval [A,B],for non-zero - TRACE traces the function + I = BSIMP('F',A,B,N,EPSILON,TRACE,TRACETOL) integrates to a + relative error of EPSILON, + beginning with n subdivisions of the interval [A,B],for non-zero + TRACE traces the function evaluations with a point plot. [I,cnt] = BSIMP(F,a,b,epsilon) also returns a function evaluation count.% Roberto Minguez Solana% Copyright (c) 2001 by Universidad de Cantabria """ if n is None: - n = int(365*8*(b-a)) + n = int(365 * 8 * (b - a)) # The number of initial subintervals must be pair if n % 2 != 0: n = n + 1 - + # Step 1: - h = (b-a)/n + h = (b - a) / n # Step 3: - x = np.linspace(a, b, n+1) + x = np.linspace(a, b, n + 1) y = fun(x) # print("y size: ", y.size) # print("n: ", n) - # TODO: + # TODO: # if trace: auxI = 0 @@ -7196,40 +7341,38 @@ def bsimp(fun, a, b, n=None, epsilon=1e-8, trace=0): cnt = n # Step 4 - integral = (ainit + 4 * auxI1 + 2*auxI2) * h/3 + integral = (ainit + 4 * auxI1 + 2 * auxI2) * h / 3 auxtot = auxI1 + auxI2 # Step 5 - integral1 = integral + epsilon*2 + integral1 = integral + epsilon * 2 j = 2 # Step 6 - error = 1 + error = 1 while error > epsilon: - cnt = cnt + int(j*n//2) + cnt = cnt + int(j * n // 2) # Step 7 integral1 = integral # Step 8 - x = np.linspace(a+h/j, b-h/j, int(j*n//2)) + x = np.linspace(a + h / j, b - h / j, int(j * n // 2)) y = fun(x) - auxI = auxI + sum(y[0:int(j*n//2)]) + auxI = auxI + sum(y[0 : int(j * n // 2)]) # Step 9 - integral = (ainit + 4 * auxI + 2*auxtot) * h/(3*j) + integral = (ainit + 4 * auxI + 2 * auxtot) * h / (3 * j) # Step 10 - j = j*2 + j = j * 2 # Step 11 auxtot = auxtot + auxI auxI = 0 # Error if np.abs(integral1) > 1e-5: - error = np.abs((integral-integral1)/integral1) + error = np.abs((integral - integral1) / integral1) else: - error = np.abs(integral-integral1) - + error = np.abs(integral - integral1) return integral - \ No newline at end of file From f0716e82764789f40523f34c5c5bd1dab03722eb Mon Sep 17 00:00:00 2001 From: Colladito Date: Wed, 8 Oct 2025 13:45:40 +0200 Subject: [PATCH 12/18] [VCF] Add log info to debug --- bluemath_tk/distributions/nonstat_gev.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index 1f97158..e76de31 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -6702,6 +6702,7 @@ def _aggquantile( for il in range(m): # function of z whose root we want def F(z): + self.logger.debug("Inicio quad()") integ, _ = quad( lambda x: self._fzeroquanint( x, @@ -6728,6 +6729,7 @@ def F(z): float(t0[il]), float(t1[il]), ) + self.logger.debug("Fin quad()") return integ + np.log(q[il]) / 12.0 try: @@ -6783,6 +6785,7 @@ def _fzeroquanint( ------ zn : np.ndarray """ + self.logger.debug("Inicio _fzeroquanint()") if ktold is None: ktold = self.kt @@ -6852,6 +6855,7 @@ def _fzeroquanint( # GUMBEL case zn[posG] = np.exp(-xn[posG]) + self.logger.debug("Fin _fzeroquanint()") return zn def _fzeroderiquanint( From 369be6dbe6f5723d44fafde56663f70d345514aa Mon Sep 17 00:00:00 2001 From: Colladito Date: Wed, 8 Oct 2025 20:14:30 +0200 Subject: [PATCH 13/18] [VCF] Decrease errors in solving numerically the integral of aggregated quantiles --- bluemath_tk/distributions/nonstat_gev.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index e76de31..1efbaa3 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -6728,6 +6728,8 @@ def F(z): ), float(t0[il]), float(t1[il]), + epsabs = 1e-6, # Absolute error in the numerical integration + epsrel = 1e-6, # Relative error in the numerical integration ) self.logger.debug("Fin quad()") return integ + np.log(q[il]) / 12.0 From 0ec71855a17da7df7bc76406ae6f5708a8ccf379 Mon Sep 17 00:00:00 2001 From: Colladito Date: Thu, 9 Oct 2025 13:05:48 +0200 Subject: [PATCH 14/18] [VCF] Update plot and using numba to make faster the aggregated return periods --- bluemath_tk/distributions/nonstat_gev.py | 566 +++++++++++++++-------- 1 file changed, 382 insertions(+), 184 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index 1efbaa3..43918c5 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -10,6 +10,34 @@ from ..core.models import BlueMathModel from ..core.plotting.colors import default_colors +from numba import njit, prange + + +@njit +def search(times: np.ndarray, values: np.ndarray, xs) -> np.ndarray: + """ + Function to search the nearest value of certain time to use in self._parametro function + + Parameters + ---------- + times : np.ndarray + Times when covariates are known + values : np.ndarray + Values of the covariates at those times + """ + n = times.shape[0] + yin = np.zeros_like(xs) + pos = 0 + for j in range(xs.size): + found = 0 + while found == 0 and pos < n: + if xs[j] < times[pos]: + yin[j] = values[pos] + found = 1 + else: + pos += 1 + + return yin class NonStatGEV(BlueMathModel): @@ -4514,6 +4542,16 @@ def _update_params(self, **kwargs: dict) -> None: self.gamma_cov = kwargs.get("gamma_cov", np.empty(0)) self.xopt = kwargs.get("x", None) + + @staticmethod + @njit + def valva(x,y): + return x+y + + @staticmethod + @njit + def cuadrado(x): + return x * x def _parametro( self, @@ -4525,9 +4563,44 @@ def _parametro( indicesint: Optional[np.ndarray] = None, times: Optional[np.ndarray] = None, x: Optional[np.ndarray] = None, + ) -> np.ndarray: + + if beta is None: + beta = np.empty(0) + if betaT is None or betaT.size == 0: + betaT = np.empty(0) + ntend = 0 + else: + ntend = 1 + if beta_cov is None: + beta_cov = np.empty(0) + if covariates is None: + covariates = np.empty((0, 0)) + if indicesint is None: + indicesint = np.empty(0) + if times is None: + times = np.empty(0) + if x is not None: + x = np.asarray([x]) + else: + x = self.t + + return self.parametro(beta0, beta, betaT, beta_cov, covariates, indicesint, times,x,ntend) + + @staticmethod + @njit + def parametro( + beta0: Optional[float] = None, + beta: Optional[np.ndarray] = None, + betaT: Optional[float] = None, + beta_cov: Optional[np.ndarray] = None, + covariates: Optional[np.ndarray] = None, + indicesint: Optional[np.ndarray] = None, + times: Optional[np.ndarray] = None, + x: Optional[np.ndarray] = None, + ntend: Optional[int] = None, ) -> np.ndarray: - """ - This function computes the location, scale and shape parameters for given parameters. Expressions by (2)-(3) in the paper + """ This function computes the location, scale and shape parameters for given parameters. Expressions by (2)-(3) in the paper Parameters ---------- @@ -4545,7 +4618,7 @@ def _parametro( Covariate mean values in the integral interval times : np.ndarray, optional Times when covariates are known, used to find the nearest value using self._search - x : np.ndarray, optional + t : np.ndarray, optional Specific time point to evaluate the parameters at, if None, uses the times given Returns @@ -4553,28 +4626,8 @@ def _parametro( y : np.ndarray Values of the parameter """ - - if beta is None: - beta = np.empty(0) - if betaT is None or betaT.size == 0: - betaT = np.empty(0) - ntend = 0 - else: - ntend = 1 - if beta_cov is None: - beta_cov = np.empty(0) - if covariates is None: - covariates = np.empty((0, 0)) - if indicesint is None: - indicesint = np.empty(0) - if times is None: - times = np.empty(0) - if x is not None: - t = np.array([x]) - else: - t = self.t - - m = len(t) + + m = len(x) na, nind = covariates.shape nparam = beta.size @@ -4590,31 +4643,31 @@ def _parametro( # Adding the harmonic part if nparam > 0: - for i in range(nparam // 2): + for i in prange(nparam // 2): y = ( y - + beta[2 * i] * np.cos((i + 1) * 2 * np.pi * t) - + beta[2 * i + 1] * np.sin((i + 1) * 2 * np.pi * t) + + beta[2 * i] * np.cos((i + 1) * 2 * np.pi * x) + + beta[2 * i + 1] * np.sin((i + 1) * 2 * np.pi * x) ) # Adding the tendency part if ntend > 0: - y = y + betaT * t + y = y + betaT * x # Adding the covariate part if nind > 0: if indicesint.shape[0] > 0: if times.shape[0] == 0: - for i in range(nind): + for i in prange(nind): y = y + beta_cov[i] * indicesint[i] else: - for i in range(nind): - indicesintaux = self._search( - times, covariates[:, i], np.asarray(x).flatten() + for i in prange(nind): + indicesintaux = search( + times, covariates[:, i], x.flatten() ) y = y + beta_cov[i] * indicesintaux else: - for i in range(nind): + for i in prange(nind): y = y + beta_cov[i] * covariates[:, i] return y @@ -4644,6 +4697,7 @@ def _search(self, times: np.ndarray, values: np.ndarray, xs) -> np.ndarray: return yin + def _evaluate_params( self, beta0: Optional[float] = None, @@ -4971,14 +5025,16 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 rp100 = self._quantile(1 - 1 / 100) if ( - self.betaT is None - and self.alphaT is None + self.betaT.size == 0 + and self.alphaT.size == 0 + and self.gammaT.size == 0 and self.beta_cov.size == 0 and self.alpha_cov.size == 0 + and self.gamma_cov.size == 0 ): - t_ord = np.argsort(t_anual) + t_ord = np.argsort(t_anual) fig, ax1 = plt.subplots(figsize=(10, 6)) - l0 = ax1.plot( + """ax1.plot( t_anual[t_ord], self.xt[t_ord], marker="+", @@ -4987,7 +5043,6 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 markersize=5, label="Data", ) - ax2 = ax1.twinx() ax1.plot( t_anual[t_ord], mut[t_ord], @@ -4998,33 +5053,146 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 ) ax1.fill_between( t_anual[t_ord], - ci_low_mut[t_ord], - ci_up_mut[t_ord], - color=self.colors[0], + mut[t_ord] - psit[t_ord], + mut[t_ord] + psit[t_ord], + label=r"Location $\pm$ Scale", + color="tab:blue", alpha=0.3, ) - ax2.plot( + # ax1.fill_between( + # t_anual[t_ord], + # ci_low_mut[t_ord], + # ci_up_mut[t_ord], + # color=self.colors[0], + # alpha=0.3, + # ) + # ax2.plot( + # t_anual[t_ord], + # psit[t_ord], + # label="Scale", + # linewidth=2, + # color=self.colors[1], + # alpha=1, + # ) + # ax2.fill_between( + # t_anual[t_ord], + # ci_low_psit[t_ord], + # ci_up_psit[t_ord], + # color=self.colors[1], + # alpha=0.3, + # ) + + ax1.plot( t_anual[t_ord], - psit[t_ord], - label="Scale", + quan95[t_ord], + linestyle="dashed", + color=self.colors[2], + markersize=5, + label="95th Quan" + ) + + month_initials = ["J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"] + month_positions = [(i + 0.5) / 12 for i in range(12)] + + rt_10 = np.zeros(12) + rt_50 = np.zeros(12) + rt_100 = np.zeros(12) + for i in range(12): + rt_10[i] = self._aggquantile(1-1/10, i/12, (i+1)/12) # 10-year return level at each year + rt_50[i] = self._aggquantile(1-1/50, i/12, (i+1)/12) # 50-year return level at each year + rt_100[i] = self._aggquantile(1-1/100, i/12, (i+1)/12) # 100-year return level + + ax1.plot( + month_positions, rt_10, linestyle="-", linewidth=1, label="10 years", color="tab:red" + ) + ax1.plot( + month_positions, rt_50, linestyle="-", linewidth=1, label="50 years", color="tab:purple" + ) + ax1.plot( + month_positions, rt_100, linestyle="-", linewidth=1, label="100 years", color="tab:green" + )""" + + ############# + month_initials = ["J", "A", "S", "O", "N", "D", "J", "F", "M", "A", "M", "J"] + month_positions = [(i + 0.5) / 12 for i in range(12)] + + # Reorder the return periods to start from July + rt_10 = np.zeros(12) + rt_50 = np.zeros(12) + rt_100 = np.zeros(12) + for i in range(12): + # Map i to the correct month index (July = 0, June = 11) + month_idx = (i + 6) % 12 + rt_10[i] = self._aggquantile(1-1/10, month_idx/12, (month_idx+1)/12) + rt_50[i] = self._aggquantile(1-1/50, month_idx/12, (month_idx+1)/12) + rt_100[i] = self._aggquantile(1-1/100, month_idx/12, (month_idx+1)/12) + + # For the data points, shift the time values by 0.5 years + t_shifted = t_anual.copy() + t_shifted[t_anual < 0.5] += 0.5 + t_shifted[t_anual >= 0.5] -= 0.5 + + # Use t_shifted for plotting data points + ax1.plot( + t_shifted[t_ord], + self.xt[t_ord], + marker="+", + linestyle="None", + color="black", + markersize=5, + label="Data", + ) + + # Use t_shifted for other lines as well + ax1.plot( + t_shifted[t_ord], + mut[t_ord], + label="Location", linewidth=2, - color=self.colors[1], + color=self.colors[0], alpha=1, ) - ax2.fill_between( + + ax1.fill_between( t_anual[t_ord], - ci_low_psit[t_ord], - ci_up_psit[t_ord], - color=self.colors[1], + mut[t_ord] - psit[t_ord], + mut[t_ord] + psit[t_ord], + label=r"Location $\pm$ Scale", + color="tab:blue", alpha=0.3, ) + ax1.plot( - t_anual[t_ord], - quan95[t_ord], - linestyle="dashed", - color=self.colors[2], - markersize=5, + month_positions, rt_10, linestyle="-", linewidth=1, label="10 years", color="tab:red" ) + ax1.plot( + month_positions, rt_50, linestyle="-", linewidth=1, label="50 years", color="tab:purple" + ) + ax1.plot( + month_positions, rt_100, linestyle="-", linewidth=1, label="100 years", color="tab:green" + ) + + + + ############# + + ax1.set_title(f"Parameters Evolution ({self.var_name})") + ax1.set_xlabel("Time") + ax1.set_ylabel(r"$\mu_t$") + ax1.grid(True) + ax1.legend(loc="best") + ax1.margins(x=0.01) + ax1.set_xticks(month_positions, month_initials) + if save: + plt.savefig( + f"Figures/Adjustment_Evolution_{self.var_name}.png", + dpi=300, + ) + plt.show() + + + + else: fig, ax1 = plt.subplots(figsize=(20, 6)) ax1.plot( @@ -5054,25 +5222,25 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 alpha=0.4, ) - ax1.plot( - self.t + init_year, - rp10, - label="10 years", - linewidth=1, - color="tab:red", - linestyle="--", - alpha=0.9, - ) + # ax1.plot( + # self.t + init_year, + # rp10, + # label="10 years", + # linewidth=1, + # color="tab:red", + # linestyle="--", + # alpha=0.9, + # ) - ax1.plot( - self.t + init_year, - rp100, - label="100 years", - linewidth=1, - color="tab:green", - linestyle="--", - alpha=0.9, - ) + # ax1.plot( + # self.t + init_year, + # rp100, + # label="100 years", + # linewidth=1, + # color="tab:green", + # linestyle="--", + # alpha=0.9, + # ) ax1.plot( self.t + init_year, @@ -5088,14 +5256,20 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 # rt_10 = np.zeros(n_years) # for year in range(n_years): # rt_10[year] = self._aggquantile(1-1/10, year, year+1) # 10-year return level at each year + # rt_50 = np.zeros(n_years) + # for year in range(n_years): + # rt_50[year] = self._aggquantile(1-1/50, year, year+1) # 50-year return level at each year # rt_100 = np.zeros(n_years) # for year in range(n_years): - # rt_100[year] = self._aggquantile(1-1/100, year, year+1) + # rt_100[year] = self._aggquantile(1-1/100, year, year+1) # 100-year return level at each year # ax1.plot( # np.arange(init_year,init_year+n_years), rt_10, linestyle="-", linewidth=1, label="10 years", color="tab:red" # ) # ax1.plot( + # np.arange(init_year,init_year+n_years), rt_50, linestyle="-", linewidth=1, label="50 years", color="tab:purple" + # ) + # ax1.plot( # np.arange(init_year, init_year+n_years), rt_100, linestyle="-", linewidth=1, label="100 years", color="tab:green" # ) @@ -5115,13 +5289,15 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 plt.savefig(f"Figures/Adjustment_Evolution_{self.var_name}.png", dpi=300) plt.show() + + ###### 1st Year PLOT month_initials = ["J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"] month_positions = [(i + 0.5) / 12 for i in range(12)] #### Creating the first year plot mask_year = (self.t >= 0) & (self.t <= 1) fig, ax1 = plt.subplots(figsize=(10, 6)) - l0 = ax1.plot( + ax1.plot( self.t[mask_year], self.xt[mask_year], marker="+", @@ -5130,8 +5306,7 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 markersize=5, label=f"{self.var_name}", ) - # ax2 = ax1.twinx() - l1 = ax1.plot( + ax1.plot( self.t[mask_year], mut[mask_year], label=r"$\mu_t$", @@ -5139,7 +5314,6 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 color=self.colors[0], alpha=1, ) - # ax1.fill_between(self.t[mask_year], ci_low_mut[mask_year], ci_up_mut[mask_year], color=self.colors[0], alpha=0.3) uppermaxs = np.maximum( mut[mask_year] - psit[mask_year], mut[mask_year] + psit[mask_year] ) @@ -5149,9 +5323,7 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 ax1.fill_between( self.t[mask_year], lowermins, uppermaxs, color=self.colors[0], alpha=0.3 ) - # l2 = ax2.plot(self.t[mask_year], psit[mask_year], label=r'$\psi_t$', linewidth=2, color=self.colors[1], alpha=1) - # ax2.fill_between(self.t[mask_year], ci_low_psit[mask_year], ci_up_psit[mask_year], color=self.colors[1], alpha=0.3) - l3 = ax1.plot( + ax1.plot( self.t[mask_year], quan95[mask_year], linestyle="dashed", @@ -5161,11 +5333,8 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 ax1.set_title(f"Location and Scale Parameters (First Year) ({self.var_name})") ax1.set_xlabel("Time") ax1.set_ylabel(r"$\mu_t$") - # ax2.set_ylabel(r'$\psi_t$') ax1.grid(True) - # handles = [art for art in l0 + l1 + l2 + l3 if not art.get_label().startswith('_')] - handles = [art for art in l0 + l1 + l3 if not art.get_label().startswith("_")] - ax1.legend(handles=handles, loc="best") + ax1.legend(loc="best") ax1.margins(x=0.01) plt.xticks(month_positions, month_initials) if save: @@ -5175,72 +5344,81 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 ) plt.show() - # Creating the first monthly plot if not monthly or annual data - mask_month = (self.t >= 0) & (self.t <= 1 / 12) - if sum(mask_month) > 1: - fig, ax1 = plt.subplots(figsize=(10, 6)) - l0 = ax1.plot( - self.t[mask_month], - self.xt[mask_month], - marker="+", - linestyle="None", - color="black", - markersize=5, - label=f"{self.var_name}", - ) - # ax2 = ax1.twinx() - l1 = ax1.plot( - self.t[mask_month], - mut[mask_month], - label=r"$\mu_t$", - linewidth=2, - color=self.colors[0], - alpha=1, - ) - uppermaxs = np.maximum( - mut[mask_month] - psit[mask_month], mut[mask_month] + psit[mask_month] - ) - lowermins = np.minimum( - mut[mask_month] - psit[mask_month], mut[mask_month] + psit[mask_month] - ) - ax1.fill_between( - self.t[mask_month], - lowermins, - uppermaxs, - color=self.colors[0], - alpha=0.3, - label="Location +- scale", - ) - # l2 = ax2.plot(self.t[mask_month], psit[mask_month], label=r'$\psi_t$', linewidth=2, color=self.colors[1], alpha=1) - # ax2.fill_between(self.t[mask_month], ci_low_psit[mask_month], ci_up_psit[mask_month], color=self.colors[1], - # alpha=0.3) - l3 = ax1.plot( - self.t[mask_month], - quan95[mask_month], - linestyle="dashed", - color=self.colors[2], - label="95th Quantile", - ) - ax1.set_title( - f"Location and Scale Parameters (First Month) ({self.var_name})" - ) - ax1.set_xlabel("Time (yearly scale)") - ax1.set_ylabel(r"$\mu_t$") - # ax2.set_ylabel(r'$\psi_t$') - ax1.grid(True) - # handles = [art for art in l0 + l1 + l2 + l3 if not art.get_label().startswith('_')] - handles = [ - art for art in l0 + l1 + l3 if not art.get_label().startswith("_") - ] - ax1.legend(handles=handles, loc="best") - ax1.margins(x=0.01) - if save: - plt.savefig( - f"Figures/Evolution_Location_FirstYear{self.var_name}.png", dpi=300 - ) - plt.show() + ### FIRST MONTH PLOT Creating the first monthly plot if not monthly or annual data + # mask_month = (self.t >= 0) & (self.t <= 1 / 12) + # if sum(mask_month) > 1: + # fig, ax1 = plt.subplots(figsize=(10, 6)) + # ax1.plot( + # self.t[mask_month], + # self.xt[mask_month], + # marker="+", + # linestyle="None", + # color="black", + # markersize=5, + # label=f"{self.var_name}", + # ) + # # ax2 = ax1.twinx() + # ax1.plot( + # self.t[mask_month], + # mut[mask_month], + # label=r"$\mu_t$", + # linewidth=2, + # color=self.colors[0], + # alpha=1, + # ) + # uppermaxs = np.maximum( + # mut[mask_month] - psit[mask_month], mut[mask_month] + psit[mask_month] + # ) + # lowermins = np.minimum( + # mut[mask_month] - psit[mask_month], mut[mask_month] + psit[mask_month] + # ) + # ax1.fill_between( + # self.t[mask_month], + # lowermins, + # uppermaxs, + # color=self.colors[0], + # alpha=0.3, + # label="Location +- scale", + # ) + + # # ax1.plot( + # # self.t[mask_month], + # # rt_10[mask_month], + # # linestyle="dashed", + # # color=self.colors[2], + # # label="10 years", + # # ) + # # ax1.plot( + # # self.t[mask_month], + # # rt_50[mask_month], + # # linestyle="dashed", + # # color=self.colors[2], + # # label="50 years", + # # ) + # # ax1.plot( + # # self.t[mask_month], + # # rt_100[mask_month], + # # linestyle="dashed", + # # color=self.colors[2], + # # label="100 years", + # # ) + # ax1.set_title( + # f"Location and Scale Parameters (First Month) ({self.var_name})" + # ) + # ax1.set_xlabel("Time (yearly scale)") + # ax1.set_ylabel(r"$\mu_t$") + # # ax2.set_ylabel(r'$\psi_t$') + # ax1.grid(True) + # # handles = [art for art in l0 + l1 + l2 + l3 if not art.get_label().startswith('_')] + # ax1.legend(loc="best") + # ax1.margins(x=0.01) + # if save: + # plt.savefig( + # f"Figures/Evolution_Location_FirstYear{self.var_name}.png", dpi=300 + # ) + # plt.show() - ### Shape parameter plot + #### Shape parameter plot # Confidence interval for epst # ci_up = ( # epst + norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdepst @@ -5249,7 +5427,7 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 # epst - norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) * stdepst # ) - # LOCATION + ##### LOCATION ##### plt.figure(figsize=(20, 6)) plt.plot(self.t + init_year, mut, color="tab:blue") # plt.fill_between( @@ -5268,7 +5446,7 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 plt.savefig(f"Figures/Evolution_Location{self.var_name}.png", dpi=300) plt.show() - # SCALE + ##### SCALE ##### plt.figure(figsize=(20, 6)) plt.plot(self.t + init_year, psit, color="tab:green") # plt.fill_between( @@ -5287,7 +5465,7 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 plt.savefig(f"Figures/Evolution_Scale{self.var_name}.png", dpi=300) plt.show() - # SHAPE + ##### SHAPE ##### plt.figure(figsize=(20, 6)) plt.plot(self.t + init_year, epst, color="tab:orange") # plt.fill_between( @@ -6501,12 +6679,22 @@ def returnperiod_plot( # Anual return periods if annualplot: ax.semilogx(Ts, quanaggrA, color="black", linewidth=2, label="Annual") + # ny = int(np.ceil(self.t[-1])) + # hmax1 = np.zeros(ny) + # for j in range(ny): + # hmax1[j] = np.max( + # self.xt[np.where((self.t >= j) & (self.t < j + 1))[0]] + # ) + + # Vectorized way ny = int(np.ceil(self.t[-1])) - hmax1 = np.zeros(ny) - for j in range(ny): - hmax1[j] = np.max( - self.xt[np.where((self.t >= j) & (self.t < j + 1))[0]] - ) + # Create mask for each year's data + year_masks = [(self.t >= j) & (self.t < j + 1) for j in range(ny)] + # Calculate max values using masks and broadcasting + hmax1 = np.array([np.max(self.xt[mask]) if np.any(mask) else np.nan for mask in year_masks]) + # Remove any NaN values if present + hmax1 = hmax1[~np.isnan(hmax1)] + hmaxsort = np.sort(hmax1) ProHsmaxsort = np.arange(1, len(hmaxsort) + 1) / (len(hmaxsort) + 1) Tapprox = 1 / (1 - ProHsmaxsort) @@ -6682,22 +6870,25 @@ def _aggquantile( 0, 1, ) - # std = quad( - # lambda x: np.exp( - # self._parametro( - # alpha0, - # alpha, - # alphaT, - # alpha_cov, - # self.covariates.iloc[:, self.list_sc].values, - # cov_scint, - # self.t, - # x, - # ) - # ), - # 0, - # 1, - # )[0] + std, _ = quad( + lambda x: np.exp( + self._parametro( + alpha0, + alpha, + alphaT, + alpha_cov, + self.covariates.iloc[:, self.list_sc].values, + cov_scint, + self.t, + x, + ) + ), + 0, + 1, + ) + + a = media - 10 + b = media + 10 for il in range(m): # function of z whose root we want @@ -6728,22 +6919,29 @@ def F(z): ), float(t0[il]), float(t1[il]), - epsabs = 1e-6, # Absolute error in the numerical integration - epsrel = 1e-6, # Relative error in the numerical integration + epsabs=1e-5, epsrel=1e-5 ) self.logger.debug("Fin quad()") return integ + np.log(q[il]) / 12.0 - + try: - # root finding: start near `media`, bracket if possible + # sol = root_scalar( + # F, + # x0=media, + # x1=media + 1.0, # secant starting points + # method="secant", + # xtol=1e-6, + # rtol=1e-6, + # maxiter=200, + # ) + sol = root_scalar( F, - x0=media, - x1=media + 1.0, # secant starting points - method="secant", - xtol=1e-6, - rtol=1e-6, - maxiter=200, + bracket=(a, b), + method="toms748", + xtol=1e-6, + rtol=1e-6, + maxiter=100 ) if sol.converged: if abs(F(sol.root)) < 1e-2: From 43cc8f5c2e22796c22af3c83226be2fde4c4019a Mon Sep 17 00:00:00 2001 From: Colladito Date: Thu, 9 Oct 2025 16:59:33 +0200 Subject: [PATCH 15/18] [VCF] Add numba and make faster the return period plot. Small fixes in returnperiod_plot when only harmonics --- bluemath_tk/distributions/nonstat_gev.py | 125 +++++------------------ 1 file changed, 27 insertions(+), 98 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index 43918c5..ab8af99 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -13,7 +13,7 @@ from numba import njit, prange -@njit +@njit(fastmath=True) def search(times: np.ndarray, values: np.ndarray, xs) -> np.ndarray: """ Function to search the nearest value of certain time to use in self._parametro function @@ -25,19 +25,27 @@ def search(times: np.ndarray, values: np.ndarray, xs) -> np.ndarray: values : np.ndarray Values of the covariates at those times """ - n = times.shape[0] - yin = np.zeros_like(xs) - pos = 0 - for j in range(xs.size): - found = 0 - while found == 0 and pos < n: - if xs[j] < times[pos]: - yin[j] = values[pos] - found = 1 - else: - pos += 1 - - return yin + # n = times.shape[0] + # yin = np.zeros_like(xs) + # pos = 0 + # for j in range(xs.size): + # found = 0 + # while found == 0 and pos < n: + # if xs[j] < times[pos]: + # yin[j] = values[pos] + # found = 1 + # else: + # pos += 1 + + # return yin + + # Get insertion indices + idx = np.searchsorted(times, xs, side="right") + + # Clip to valid range (so no out-of-bounds) + idx = np.clip(idx, 0, len(times) - 1) + + return values[idx] class NonStatGEV(BlueMathModel): @@ -99,6 +107,8 @@ def __init__( """ Initiliaze the Non-Stationary GEV. """ + super().__init__() + debug=1 self.set_logger_name( name=self.__class__.__name__, level="DEBUG" if debug else "INFO" @@ -4588,7 +4598,7 @@ def _parametro( return self.parametro(beta0, beta, betaT, beta_cov, covariates, indicesint, times,x,ntend) @staticmethod - @njit + @njit(fastmath=True) def parametro( beta0: Optional[float] = None, beta: Optional[np.ndarray] = None, @@ -5032,85 +5042,7 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 and self.alpha_cov.size == 0 and self.gamma_cov.size == 0 ): - t_ord = np.argsort(t_anual) fig, ax1 = plt.subplots(figsize=(10, 6)) - """ax1.plot( - t_anual[t_ord], - self.xt[t_ord], - marker="+", - linestyle="None", - color="black", - markersize=5, - label="Data", - ) - ax1.plot( - t_anual[t_ord], - mut[t_ord], - label="Location", - linewidth=2, - color=self.colors[0], - alpha=1, - ) - ax1.fill_between( - t_anual[t_ord], - mut[t_ord] - psit[t_ord], - mut[t_ord] + psit[t_ord], - label=r"Location $\pm$ Scale", - color="tab:blue", - alpha=0.3, - ) - # ax1.fill_between( - # t_anual[t_ord], - # ci_low_mut[t_ord], - # ci_up_mut[t_ord], - # color=self.colors[0], - # alpha=0.3, - # ) - # ax2.plot( - # t_anual[t_ord], - # psit[t_ord], - # label="Scale", - # linewidth=2, - # color=self.colors[1], - # alpha=1, - # ) - # ax2.fill_between( - # t_anual[t_ord], - # ci_low_psit[t_ord], - # ci_up_psit[t_ord], - # color=self.colors[1], - # alpha=0.3, - # ) - - ax1.plot( - t_anual[t_ord], - quan95[t_ord], - linestyle="dashed", - color=self.colors[2], - markersize=5, - label="95th Quan" - ) - - month_initials = ["J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"] - month_positions = [(i + 0.5) / 12 for i in range(12)] - - rt_10 = np.zeros(12) - rt_50 = np.zeros(12) - rt_100 = np.zeros(12) - for i in range(12): - rt_10[i] = self._aggquantile(1-1/10, i/12, (i+1)/12) # 10-year return level at each year - rt_50[i] = self._aggquantile(1-1/50, i/12, (i+1)/12) # 50-year return level at each year - rt_100[i] = self._aggquantile(1-1/100, i/12, (i+1)/12) # 100-year return level - - ax1.plot( - month_positions, rt_10, linestyle="-", linewidth=1, label="10 years", color="tab:red" - ) - ax1.plot( - month_positions, rt_50, linestyle="-", linewidth=1, label="50 years", color="tab:purple" - ) - ax1.plot( - month_positions, rt_100, linestyle="-", linewidth=1, label="100 years", color="tab:green" - )""" ############# month_initials = ["J", "A", "S", "O", "N", "D", "J", "F", "M", "A", "M", "J"] @@ -5131,6 +5063,7 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 t_shifted = t_anual.copy() t_shifted[t_anual < 0.5] += 0.5 t_shifted[t_anual >= 0.5] -= 0.5 + t_ord = np.argsort(t_shifted) # Use t_shifted for plotting data points ax1.plot( @@ -5154,7 +5087,7 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 ) ax1.fill_between( - t_anual[t_ord], + t_shifted[t_ord], mut[t_ord] - psit[t_ord], mut[t_ord] + psit[t_ord], label=r"Location $\pm$ Scale", @@ -5172,10 +5105,6 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 month_positions, rt_100, linestyle="-", linewidth=1, label="100 years", color="tab:green" ) - - - ############# - ax1.set_title(f"Parameters Evolution ({self.var_name})") ax1.set_xlabel("Time") ax1.set_ylabel(r"$\mu_t$") From 113edf062ddb5f5364f43389eb08f8c320e7604c Mon Sep 17 00:00:00 2001 From: Colladito Date: Thu, 9 Oct 2025 17:54:24 +0200 Subject: [PATCH 16/18] [VCF] Update search function in NonStatGEV --- bluemath_tk/distributions/nonstat_gev.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index ab8af99..a252bb2 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -39,13 +39,10 @@ def search(times: np.ndarray, values: np.ndarray, xs) -> np.ndarray: # return yin - # Get insertion indices - idx = np.searchsorted(times, xs, side="right") - - # Clip to valid range (so no out-of-bounds) - idx = np.clip(idx, 0, len(times) - 1) - - return values[idx] + idx = np.searchsorted(times, xs, side='right') + mask = idx < len(times) + + return values[idx[mask]] class NonStatGEV(BlueMathModel): From 24a1749e44f18cc7b6241696fd3144ec558bf0d4 Mon Sep 17 00:00:00 2001 From: Colladito Date: Thu, 9 Oct 2025 18:28:48 +0200 Subject: [PATCH 17/18] [VCF] Update standard errors --- bluemath_tk/distributions/nonstat_gev.py | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index a252bb2..64f9668 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -1223,8 +1223,9 @@ def auto_adjust(self, max_iter: int = 1000, plot: bool = False) -> dict: self.invI0 = np.linalg.inv(-Hxx) fit_result["invI0"] = self.invI0 - std_param = np.sqrt(np.diag(self.invI0)) - fit_result["std_param"] = std_param + std_params = np.sqrt(np.diag(self.invI0)) + self.std_params = std_params + fit_result["std_params"] = std_params if plot: self.plot() @@ -2988,8 +2989,10 @@ def fit( self.invI0 = np.linalg.inv(-Hxx) fit_result["invI0"] = self.invI0 - std_param = np.sqrt(np.diag(self.invI0)) - fit_result["std_param"] = std_param + + std_params = np.sqrt(np.diag(self.invI0)) + self.std_params = std_params + fit_result["std_param"] = std_params if plot: self.plot() @@ -4550,16 +4553,6 @@ def _update_params(self, **kwargs: dict) -> None: self.xopt = kwargs.get("x", None) - @staticmethod - @njit - def valva(x,y): - return x+y - - @staticmethod - @njit - def cuadrado(x): - return x * x - def _parametro( self, beta0: Optional[float] = None, @@ -7290,7 +7283,7 @@ def summary(self): """ Print a summary of the fitted model, including parameter estimates, standard errors and fit statistics. """ - std_params = np.sqrt(np.diag(self.invI0)) + std_params = self.std_params param_idx = 0 z_norm = norm.ppf(1 - (1 - self.quanval) / 2, loc=0, scale=1) From 993437694668d2ca1fab398422d593dae9729abd Mon Sep 17 00:00:00 2001 From: Colladito Date: Thu, 9 Oct 2025 20:58:42 +0200 Subject: [PATCH 18/18] [VCF] Stable version NonStatGEV. Vectorize parametro function and remove _search function. --- bluemath_tk/distributions/nonstat_gev.py | 121 ++++++++++------------- 1 file changed, 54 insertions(+), 67 deletions(-) diff --git a/bluemath_tk/distributions/nonstat_gev.py b/bluemath_tk/distributions/nonstat_gev.py index 64f9668..03d40ac 100644 --- a/bluemath_tk/distributions/nonstat_gev.py +++ b/bluemath_tk/distributions/nonstat_gev.py @@ -13,36 +13,36 @@ from numba import njit, prange -@njit(fastmath=True) -def search(times: np.ndarray, values: np.ndarray, xs) -> np.ndarray: - """ - Function to search the nearest value of certain time to use in self._parametro function - - Parameters - ---------- - times : np.ndarray - Times when covariates are known - values : np.ndarray - Values of the covariates at those times - """ - # n = times.shape[0] - # yin = np.zeros_like(xs) - # pos = 0 - # for j in range(xs.size): - # found = 0 - # while found == 0 and pos < n: - # if xs[j] < times[pos]: - # yin[j] = values[pos] - # found = 1 - # else: - # pos += 1 - - # return yin - - idx = np.searchsorted(times, xs, side='right') - mask = idx < len(times) - - return values[idx[mask]] +# @njit(fastmath=True) +# def search(times: np.ndarray, values: np.ndarray, xs) -> np.ndarray: +# """ +# Function to search the nearest value of certain time to use in self._parametro function + +# Parameters +# ---------- +# times : np.ndarray +# Times when covariates are known +# values : np.ndarray +# Values of the covariates at those times +# """ +# # n = times.shape[0] +# # yin = np.zeros_like(xs) +# # pos = 0 +# # for j in range(xs.size): +# # found = 0 +# # while found == 0 and pos < n: +# # if xs[j] < times[pos]: +# # yin[j] = values[pos] +# # found = 1 +# # else: +# # pos += 1 + +# # return yin + +# idx = np.searchsorted(times, xs, side='right') +# mask = idx < len(times) + +# return values[idx[mask]] class NonStatGEV(BlueMathModel): @@ -4617,7 +4617,7 @@ def parametro( indicesint : np.ndarray, optional Covariate mean values in the integral interval times : np.ndarray, optional - Times when covariates are known, used to find the nearest value using self._search + Times when covariates are known, used to find the nearest value t : np.ndarray, optional Specific time point to evaluate the parameters at, if None, uses the times given @@ -4658,45 +4658,32 @@ def parametro( if nind > 0: if indicesint.shape[0] > 0: if times.shape[0] == 0: - for i in prange(nind): - y = y + beta_cov[i] * indicesint[i] + # for i in prange(nind): + # y = y + beta_cov[i] * indicesint[i] + y = y + indicesint @ beta_cov else: - for i in prange(nind): - indicesintaux = search( - times, covariates[:, i], x.flatten() - ) - y = y + beta_cov[i] * indicesintaux + # for i in prange(nind): + # indicesintaux = search( + # times, covariates[:, i], x.flatten() + # ) + # y = y + beta_cov[i] * indicesintaux + idx = np.searchsorted(times, x, side='right') + valid = idx < times.size + + y_add = np.zeros_like(x, dtype=np.float64) + if np.any(valid): + # pick rows from covariates and do one matvec + A = covariates[idx[valid], :] # (k, nind) + y_add[valid] = A @ beta_cov # (k,) + + y = y + y_add else: - for i in prange(nind): - y = y + beta_cov[i] * covariates[:, i] + # for i in prange(nind): + # y = y + beta_cov[i] * covariates[:, i] + y = y + covariates @ beta_cov - return y - - def _search(self, times: np.ndarray, values: np.ndarray, xs) -> np.ndarray: - """ - Function to search the nearest value of certain time to use in self._parametro function - - Parameters - ---------- - times : np.ndarray - Times when covariates are known - values : np.ndarray - Values of the covariates at those times - """ - n = times.shape[0] - yin = np.zeros_like(xs) - pos = 0 - for j in range(xs.size): - found = 0 - while found == 0 and pos < n: - if xs[j] < times[pos]: - yin[j] = values[pos] - found = 1 - else: - pos += 1 - - return yin + return y def _evaluate_params( self, @@ -5170,7 +5157,7 @@ def plot(self, return_plot: bool = False, save: bool = False, init_year: int = 0 alpha=0.9, ) - # TODO: Add aggregated return period lines + # Aggregated return period lines # n_years = int(np.ceil(self.t[-1])) # rt_10 = np.zeros(n_years) # for year in range(n_years):