From 0c61540fb347b0db85514f4d9b606f1ff5dc4675 Mon Sep 17 00:00:00 2001 From: thurinj Date: Mon, 19 Aug 2024 21:42:28 +1000 Subject: [PATCH] Added parabola fit to the depth plot MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit I added a parabola to the depth plot, with an estimate of the best fitting depth at the minimum of the parabola. It will print the best value and the ± error from a 5% uncertainty threshold. --- mtuq/graphics/uq/_matplotlib.py | 74 +++++++++++++++++++++++++++++++-- 1 file changed, 70 insertions(+), 4 deletions(-) diff --git a/mtuq/graphics/uq/_matplotlib.py b/mtuq/graphics/uq/_matplotlib.py index 9daac087..6c6cddd8 100644 --- a/mtuq/graphics/uq/_matplotlib.py +++ b/mtuq/graphics/uq/_matplotlib.py @@ -9,8 +9,9 @@ from mtuq.graphics.uq import _nothing_to_plot from mtuq.graphics._gmt import read_cpt, _cpt_path from mtuq.graphics.uq._gmt import _parse_vw, _parse_lune_array, _parse_force -from matplotlib.colors import BoundaryNorm from matplotlib import ticker +from matplotlib.colors import BoundaryNorm +import matplotlib.patheffects as path_effects from mtuq.util.math import wrap_180, to_delta_gamma, to_mij from mpl_toolkits.axes_grid1 import make_axes_locatable @@ -391,10 +392,75 @@ def _plot_depth_matplotlib(filename, depths, values, # String formating to display only x.xx values magnitudes = ['{:.2f}'.format(m) for m in magnitudes] for i, magnitude in enumerate(magnitudes): - pyplot.text(depths[i], + text = pyplot.text(depths[i], values[i]+0.045*(values.max()-values.min()), - magnitude, fontsize=fontsize-6, ha='center') + magnitude, fontsize=fontsize-6, ha='center', zorder=2000) + # Adding a white stroke around the text + text.set_path_effects([path_effects.Stroke(linewidth=3, foreground='white'), + path_effects.Normal()]) + + # Fit a parabola to the best-fit depth + imin = np.argmin(values) + + if imin == 0: + # Best fit at the left edge, use first three points + x1, y1 = depths[0], values[0] + x2, y2 = depths[1], values[1] + x3, y3 = depths[2], values[2] + elif imin == len(depths) - 1: + # Best fit at the right edge, use last three points + x1, y1 = depths[-3], values[-3] + x2, y2 = depths[-2], values[-2] + x3, y3 = depths[-1], values[-1] + else: + # Best fit not at an edge, use the point and its neighbors + x1, y1 = depths[imin-1], values[imin-1] + x2, y2 = depths[imin], values[imin] + x3, y3 = depths[imin+1], values[imin+1] + + # Calculate coefficients for the parabola (y = ax^2 + bx + c) + denom = (x1-x2)*(x1-x3)*(x2-x3) + a = (x3*(y2-y1) + x2*(y1-y3) + x1*(y3-y2)) / denom + b = (x3**2*(y1-y2) + x2**2*(y3-y1) + x1**2*(y2-y3)) / denom + c = (x2*x3*(x2-x3)*y1 + x3*x1*(x3-x1)*y2 + x1*x2*(x1-x2)*y3) / denom + + + xlim = ax.get_xlim() + ylim = ax.get_ylim() + # Generate parabola points + x_parabola = np.linspace(min(depths)-0.1*(max(depths)-min(depths)), + max(depths)+0.1*(max(depths)-min(depths)), 100) + + y_parabola = a*x_parabola**2 + b*x_parabola + c + + best_fit_depth = -b / (2*a) + best_fit_value = a * best_fit_depth**2 + b * best_fit_depth + c + print('Best fitting depth:', best_fit_depth) + + percentage = 0.05 # Define the percentage of the best fit value to use as a threshold + delta_y = percentage * best_fit_value + y_threshold = best_fit_value + delta_y + + # Step 3: Get the roots where the parabola equals the threshold + discriminant = np.sqrt(b**2 - 4*a*(c - y_threshold)) + + x1 = (-b + discriminant) / (2 * a) + x2 = (-b - discriminant) / (2 * a) + + # Step 4: Calculate the uncertainty k + k = abs(x2 - x1) / 2 + print('Uncertainty: ±', k) + + ax.hlines(y=y_threshold, xmin=x1, xmax=x2, colors='gray', linestyles='dashed', label='5% Threshold') + ax.plot(best_fit_depth, ylim[0], marker='v', color='white', markersize=10, + markeredgecolor='black', markeredgewidth=2, clip_on=False) + + # Plot the parabola + ax.plot(x_parabola, y_parabola, '--', color='gray') + # ax.legend() + ax.set_xlim(xlim) + ax.set_ylim(ylim) if lune_array is not None: from mtuq.graphics.beachball import _plot_beachball_matplotlib @@ -541,7 +607,7 @@ def _add_marker(axis, coords, size=250): edgecolors=[0,1,0], linewidths=1.75, clip_on=False, - zorder=100, + zorder=200, )