From b163fcc03a90ad14bf08d2ffef13e9f1c3877da6 Mon Sep 17 00:00:00 2001 From: Julien Thurin Date: Tue, 17 Feb 2026 14:40:29 +0900 Subject: [PATCH] Reimplement omega angle calculation and confidence curve from clean master. Follows PR #341 --- mtuq/graphics/__init__.py | 2 +- mtuq/graphics/uq/_matplotlib.py | 49 +++++++++++++++++++++ mtuq/graphics/uq/omega.py | 75 +++++++++++++++++++++++++++++---- 3 files changed, 116 insertions(+), 10 deletions(-) diff --git a/mtuq/graphics/__init__.py b/mtuq/graphics/__init__.py index e7e21bfb..75ea3d9b 100644 --- a/mtuq/graphics/__init__.py +++ b/mtuq/graphics/__init__.py @@ -37,7 +37,7 @@ plot_data_greens1, plot_data_greens2, plot_data_greens3 from mtuq.graphics.uq.omega import\ - plot_cdf, plot_pdf, plot_screening_curve + plot_cdf, plot_pdf, plot_screening_curve, plot_confidence_curve from mtuq.graphics.uq import\ likelihood_analysis diff --git a/mtuq/graphics/uq/_matplotlib.py b/mtuq/graphics/uq/_matplotlib.py index 3873a046..6b98268e 100644 --- a/mtuq/graphics/uq/_matplotlib.py +++ b/mtuq/graphics/uq/_matplotlib.py @@ -607,6 +607,55 @@ def _plot_omega_matplotlib(filename, omega, values, pyplot.savefig(filename) pyplot.close() +def _plot_confidence_curve_matplotlib(filename, fractional_volume, values, + title=None, xlabel='V', ylabel=r'$\mathcal{P}(V)$', figsize=(4., 4.), fontsize=18.): + + fractional_volume = np.insert(fractional_volume, 0, 0) + values = np.insert(values, 0, 0) + fractional_volume = np.append(fractional_volume, 1) + values = np.append(values, 1) + fig, ax = pyplot.subplots(figsize=figsize) + + ax.plot(fractional_volume, values, 'r-', linewidth=2.5, clip_on=False) + + ax.plot(fractional_volume, fractional_volume, linestyle='--', color='gray', linewidth=1.5) + + ax.fill_between(fractional_volume, values, 0, color='gray', alpha=0.3) + + # Make sure the plot starts at 0,0 + fractional_volume = np.insert(fractional_volume, 0, 0) + values = np.insert(values, 0, 0) + + ax.plot([1], [1], 'ko', clip_on=False) # Marker at top-right corner + ax.plot([0], [0], 'ko', clip_on=False) # Bottom-left corner marker + + + # Display in text the average value of P(V) + average = np.mean(values) + ax.text(0.72, 0.10, r'$\mathcal{{P}}_{{AV}} = {:.2f}$'.format(average), # Use raw string (r'') and LaTeX math formatting + fontsize=fontsize+2, ha='center', va='center', transform=ax.transAxes) + + # Customize tick labels + ax.set_xticks(np.linspace(0, 1, 11)) # Keep original ticks + ax.set_yticks(np.linspace(0, 1, 11)) + ax.set_xticklabels(['0' if tick == 0 else '1' if tick == 1 else '' for tick in np.linspace(0, 1, 11)], fontsize=fontsize-2) + ax.set_yticklabels(['0' if tick == 0 else '1' if tick == 1 else '' for tick in np.linspace(0, 1, 11)], fontsize=fontsize-2) + + + if title: + ax.set_title(title, fontsize=fontsize) + + if xlabel: + ax.set_xlabel(xlabel, fontsize=fontsize) + + if ylabel: + ax.set_ylabel(ylabel, fontsize=fontsize) + + ax.set_xlim(0, 1) + ax.set_ylim(0, 1) + ax.margins(x=0.02, y=0.02) + + pyplot.savefig(filename, bbox_inches='tight', dpi=300) # # utility functions # diff --git a/mtuq/graphics/uq/omega.py b/mtuq/graphics/uq/omega.py index 90629d41..2920b2e9 100644 --- a/mtuq/graphics/uq/omega.py +++ b/mtuq/graphics/uq/omega.py @@ -7,7 +7,7 @@ import numpy as np from mtuq import Force, MomentTensor -from mtuq.graphics.uq._matplotlib import _plot_omega_matplotlib +from mtuq.graphics.uq._matplotlib import _plot_omega_matplotlib, _plot_confidence_curve_matplotlib from mtuq.grid_search import DataArray, DataFrame from mtuq.util import warn from mtuq.util.math import to_mij @@ -47,7 +47,7 @@ def plot_pdf(filename, df, var, m0=None, nbins=50, normalized=False, **kwargs): -def plot_cdf(filename, df, var, nbins=50, normalized=False, **kwargs): +def plot_cdf(filename, df, var, m0=None, nbins=50, normalized=False, **kwargs): """ Plots cumulative distribution function over angular distance .. rubric :: Input arguments @@ -77,7 +77,37 @@ def plot_cdf(filename, df, var, nbins=50, normalized=False, **kwargs): _plot_omega(filename, omega, np.cumsum(pdf), **kwargs) +def plot_confidence_curve(filename, df, var, m0=None, nbins=50, normalized=False, **kwargs): + """ Plots confidence curve over fractional volume + .. rubric :: Input arguments + + ``filename`` (`str`): + Name of output image file + + ``df`` (`DataFrame`): + Data structure containing moment tensors and corresponding misfit values + + ``var`` (`float` or `array`): + Data variance + + ``nbins`` (`int`): + Number of angular distance bins + + ``normalized`` (`bool`): + Normalize each angular distance bin by volume of corresponding shell? + + """ + if not isuniform(df): + warn('plot_confidence_curve requires randomly-drawn grid') + return + + omega, pdf = _calculate_pdf(df, var, m0=m0, nbins=nbins, + normalized=normalized) + _, pdf_homo = _calculate_pdf(df*0, var, m0=m0, nbins=nbins, + normalized=normalized) + + _plot_omega(filename, np.cumsum(pdf_homo/np.sum(pdf_homo)), np.cumsum(pdf/np.sum(pdf)), backend=_plot_confidence_curve_matplotlib, **kwargs) def plot_screening_curve(filename, ds, var, nbins=50, **kwargs): """ Plots explosion screening curve (maximum likelihood versus angular @@ -164,7 +194,7 @@ def _calculate_omega(df, m0=None): # extract reference vector if type(m0)==MomentTensor: # convert from lune to mij parameters - m0 = m0.as_vector() + m0 = m0.as_matrix() elif type(m0)==Force: raise NotImplementedError @@ -173,14 +203,41 @@ def _calculate_omega(df, m0=None): # assume df holds likelihoods, try maximum likelihood estimate idx = _argmax(df) m0 = m[idx,:] - - # vectorized dot product - dp = np.dot(m, m0) - dp /= np.sum(m0**2)**0.5 - dp /= np.sum(m**2, axis=1)**0.5 + m0 = np.array([[m0[0], m0[3], m0[4]], + [m0[3], m0[1], m0[5]], + [m0[4], m0[5], m0[2]]]) + + + m_tensors = np.zeros((m.shape[0], 3, 3)) + m_tensors[:, 0, 0] = m[:, 0] # Mrr + m_tensors[:, 1, 1] = m[:, 1] # Mtt + m_tensors[:, 2, 2] = m[:, 2] # Mpp + m_tensors[:, 0, 1] = m[:, 3] # Mrt + m_tensors[:, 1, 0] = m[:, 3] # Mrt (symmetric) + m_tensors[:, 0, 2] = m[:, 4] # Mrp + m_tensors[:, 2, 0] = m[:, 4] # Mrp (symmetric) + m_tensors[:, 1, 2] = m[:, 5] # Mtp + m_tensors[:, 2, 1] = m[:, 5] # Mtp (symmetric) + + # Compute the dot product of the tensors M and N + dot_product = np.einsum('...ij,...ij->...', m0, m_tensors) + + # Compute the norms of the tensors M and N + norm_M = np.sqrt(np.einsum('...ij,...ij->...', m0, m0)) + norm_N = np.sqrt(np.einsum('...ij,...ij->...', m_tensors, m_tensors)) + + # Compute the cosine of the angle between the tensors + cos_angle = dot_product / (norm_M * norm_N) + + # Clip values to the valid range for arccos (prevent errors from numerical precision) + cos_angle = cos_angle.clip(-1, 1) # return angles as NumPy array - omega = 180./np.pi * np.arccos(dp) + omega = 180./np.pi * np.arccos(cos_angle) + + # Fix nan values for identical vectors + omega[np.isnan(omega)] = 0. + return omega