19090154521

Committed 05 Nov 2025 03:19AM UTC coverage: 72.783%. Remained the same

Build # 19090154521

Build Type

push

github

Committed by

pavelkomarov

Commit Message

fixed bug where vector lengths might not agree

Run Details

0 of 1 new or added line in 1 file covered. (0.0%)

837 of 1150 relevant lines covered (72.78%)

0.73 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

11.54

/pynumdiff/utils/evaluate.py

"""Some tools to help evaluate and plot performance, used in optimization and in jupyter notebooks"""
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

from pynumdiff.utils import utility

# pylint: disable-msg=too-many-locals, too-many-arguments
def plot(x, dt, x_hat, dxdt_hat, x_truth, dxdt_truth, xlim=None, show_error=True, markersize=5):
    """Make comparison plots of 'x (blue) vs x_truth (black) vs x_hat (red)' and 'dxdt_truth
    (black) vs dxdt_hat (red)'

    :param np.array[float] x: array of noisy data
    :param float dt: a float number representing the step size
    :param np.array[float] x_hat: array of smoothed estimation of x
    :param np.array[float] dxdt_hat: array of estimated derivative
    :param np.array[float] x_truth: array of noise-free time series
    :param np.array[float] dxdt_truth: array of true derivative
    :param list[int] xlim: a list specifying range of x
    :param bool show_error: whether to show the rmse
    :param int markersize: marker size of noisy observations

    :return: Display two plots
    """
    t = np.arange(len(x))*dt
    if xlim is None:
        xlim = [t[0], t[-1]]

    fig, axes = plt.subplots(1, 2, figsize=(18, 6), constrained_layout=True)
    
    axes[0].plot(t, x_truth, '--', color='black', linewidth=3, label=r"true $x$")
    axes[0].plot(t, x, '.', color='blue', zorder=-100, markersize=markersize, label=r"noisy data")
    axes[0].plot(t, x_hat, color='red', label=r"estimated $\hat{x}$")
    axes[0].set_ylabel('Position', fontsize=18)
    axes[0].set_xlabel('Time', fontsize=18)
    axes[0].set_xlim(xlim[0], xlim[-1])
    axes[0].tick_params(axis='x', labelsize=15)
    axes[0].tick_params(axis='y', labelsize=15)
    axes[0].legend(loc='lower right', fontsize=12)

    axes[1].plot(t, dxdt_truth, '--', color='black', linewidth=3, label=r"true  $\frac{dx}{dt}$")
    axes[1].plot(t, dxdt_hat, color='red', label=r"est. $\hat{\frac{dx}{dt}}$")
    axes[1].set_ylabel('Velocity', fontsize=18)
    axes[1].set_xlabel('Time', fontsize=18)
    axes[1].set_xlim(xlim[0], xlim[-1])
    axes[1].tick_params(axis='x', labelsize=15)
    axes[1].tick_params(axis='y', labelsize=15)
    axes[1].legend(loc='lower right', fontsize=12)

    if show_error:
        _, _, rms_dxdt = rmse(x, dt, x_hat, dxdt_hat, x_truth, dxdt_truth)
        R_sqr = error_correlation(dxdt_hat, dxdt_truth)
        axes[1].text(0.05, 0.95, f"RMSE = {rms_dxdt:.2f}\n$R^2$ = {R_sqr:.2g}",
                     transform=axes[1].transAxes, fontsize=15, verticalalignment='top')
    
    return fig, axes


def plot_comparison(dt, dxdt_truth, dxdt_hat1, title1, dxdt_hat2, title2, dxdt_hat3, title3):
    """This is intended to show method performances with different choices of parameter"""
    t = np.arange(0, dt*len(dxdt_truth), dt)
    fig, axes = plt.subplots(1, 3, figsize=(22,6))

    for i,(dxdt_hat,title) in enumerate(zip([dxdt_hat1, dxdt_hat2, dxdt_hat3], [title1, title2, title3])):
        axes[i].plot(t, dxdt_truth, '--', color='black', linewidth=3, label=r"true  $\frac{dx}{dt}$")
        axes[i].plot(t, dxdt_hat, color='red', label=r"est. $\hat{\frac{dx}{dt}}$")
        if i==0: axes[i].set_ylabel('Velocity', fontsize=18)
        axes[i].set_xlabel('Time', fontsize=18)
        axes[i].tick_params(axis='x', labelsize=15)
        axes[i].tick_params(axis='y', labelsize=15)
        axes[i].set_title(title, fontsize=18)
        if i==2: axes[i].legend(loc='lower right', fontsize=12)
        _, _, rms_dxdt = rmse(np.zeros(dxdt_hat.shape), dt, np.zeros(dxdt_hat.shape), dxdt_hat, np.zeros(dxdt_hat.shape), dxdt_truth)
        R_sqr = error_correlation(dxdt_hat, dxdt_truth)
        axes[i].text(0.05, 0.95, f"RMSE = {rms_dxdt:.2f}\n$R^2$ = {R_sqr:.2g}",
                     transform=axes[i].transAxes, fontsize=15, verticalalignment='top')

    fig.tight_layout()


def rmse(x, dt, x_hat, dxdt_hat, x_truth=None, dxdt_truth=None, padding=0):
    """Evaluate x_hat based on RMSE, calculating different ones depending on whether :code:`dxdt_truth`
    and :code:`x_truth` are known.

    :param np.array[float] x: data that was differentiated
    :param float dt: step size
    :param np.array[float] x_hat: estimated (smoothed) x
    :param np.array[float] dxdt_hat: estimated xdot
    :param np.array[float] x_truth: true value of x, if known
    :param np.array[float] dxdt_truth: true value of dxdt, if known, optional
    :param int padding: number of snapshots on either side of the array to ignore when calculating
        the metric. If :code:`'auto'`, defaults to 2.5% of the size of x

    :return: tuple[float, float, float] containing\n
            - **rms_rec_x** -- RMS error between the integral of dxdt_hat and x
            - **rms_x** -- RMS error between x_hat and x_truth, returns None if x_truth is None
            - **rms_dxdt** -- RMS error between dxdt_hat and dxdt_truth, returns None if dxdt_hat is None
    """
    if np.isnan(x_hat).any():
        return np.nan, np.nan, np.nan
    if padding == 'auto':
        padding = int(0.025*len(x))
        padding = max(padding, 1)
    s = slice(padding, len(x)-padding) # slice out data we want to measure

    # RMS of dxdt and x_hat
    root = np.sqrt(s.stop - s.start)
    rms_dxdt = np.linalg.norm(dxdt_hat[s] - dxdt_truth[s]) / root if dxdt_truth is not None else None
    rms_x = np.linalg.norm(x_hat[s] - x_truth[s]) / root if x_truth is not None else None

    # RMS reconstructed x from integrating dxdt vs given raw x, available even in the absence of ground truth
    rec_x = utility.integrate_dxdt_hat(dxdt_hat, dt)
    x0 = utility.estimate_integration_constant(x, rec_x)
    rec_x = rec_x + x0
    rms_rec_x = np.linalg.norm(rec_x[s] - x[s]) / root

    return rms_rec_x, rms_x, rms_dxdt


def error_correlation(dxdt_hat, dxdt_truth, padding=0):
    """Calculate the error correlation (pearsons correlation coefficient) between the estimated
    dxdt and true dxdt

    :param np.array[float] dxdt_hat: estimated xdot
    :param np.array[float] dxdt_truth: true value of dxdt, if known, optional
    :param int padding: number of snapshots on either side of the array to ignore when calculating
        the metric. If :code:`'auto'`, defaults to 2.5% of the size of x

    :return: (float) -- r-squared correlation coefficient
    """
    if padding == 'auto':
        padding = int(0.025*len(dxdt_hat))
        padding = max(padding, 1)
    s = slice(padding, len(dxdt_hat)-padding) # slice out data we want to measure
    errors = (dxdt_hat[s] - dxdt_truth[s])
    r = stats.linregress(dxdt_truth[s] - np.mean(dxdt_truth[s]), errors)
    return r.rvalue**2


def total_variation(x, padding=0):
    """Calculate the total variation of an array. Used by optimizer.

    :param np.array[float] x: data
    :param int padding: number of snapshots on either side of the array to ignore when calculating
        the metric. If :code:`'auto'`, defaults to 2.5% of the size of x

    :return: (float) -- total variation
    """
    if np.isnan(x).any():
        return np.nan
    if padding == 'auto':
        padding = int(0.025*len(x))
        padding = max(padding, 1)
    x = x[padding:len(x)-padding]
    
    return np.linalg.norm(x[1:]-x[:-1], 1)/len(x) # normalized version of what cvxpy.tv does

1	"""Some tools to help evaluate and plot performance, used in optimization and in jupyter notebooks"""
2	import numpy as np	1✔
3	import matplotlib.pyplot as plt	1✔
4	from scipy import stats	1✔
5
6	from pynumdiff.utils import utility	1✔
7
8	# pylint: disable-msg=too-many-locals, too-many-arguments
9	def plot(x, dt, x_hat, dxdt_hat, x_truth, dxdt_truth, xlim=None, show_error=True, markersize=5):	1✔
10	"""Make comparison plots of 'x (blue) vs x_truth (black) vs x_hat (red)' and 'dxdt_truth
11	(black) vs dxdt_hat (red)'
12
13	:param np.array[float] x: array of noisy data
14	:param float dt: a float number representing the step size
15	:param np.array[float] x_hat: array of smoothed estimation of x
16	:param np.array[float] dxdt_hat: array of estimated derivative
17	:param np.array[float] x_truth: array of noise-free time series
18	:param np.array[float] dxdt_truth: array of true derivative
19	:param list[int] xlim: a list specifying range of x
20	:param bool show_error: whether to show the rmse
21	:param int markersize: marker size of noisy observations
22
23	:return: Display two plots
24	"""
NEW 25	t = np.arange(len(x))*dt	×
26	if xlim is None:	×
27	xlim = [t[0], t[-1]]	×
28
29	fig, axes = plt.subplots(1, 2, figsize=(18, 6), constrained_layout=True)	×
30
31	axes[0].plot(t, x_truth, '--', color='black', linewidth=3, label=r"true $x$")	×
32	axes[0].plot(t, x, '.', color='blue', zorder=-100, markersize=markersize, label=r"noisy data")	×
33	axes[0].plot(t, x_hat, color='red', label=r"estimated $\hat{x}$")	×
34	axes[0].set_ylabel('Position', fontsize=18)	×
35	axes[0].set_xlabel('Time', fontsize=18)	×
36	axes[0].set_xlim(xlim[0], xlim[-1])	×
37	axes[0].tick_params(axis='x', labelsize=15)	×
38	axes[0].tick_params(axis='y', labelsize=15)	×
39	axes[0].legend(loc='lower right', fontsize=12)	×
40
41	axes[1].plot(t, dxdt_truth, '--', color='black', linewidth=3, label=r"true $\frac{dx}{dt}$")	×
42	axes[1].plot(t, dxdt_hat, color='red', label=r"est. $\hat{\frac{dx}{dt}}$")	×
43	axes[1].set_ylabel('Velocity', fontsize=18)	×
44	axes[1].set_xlabel('Time', fontsize=18)	×
45	axes[1].set_xlim(xlim[0], xlim[-1])	×
46	axes[1].tick_params(axis='x', labelsize=15)	×
47	axes[1].tick_params(axis='y', labelsize=15)	×
48	axes[1].legend(loc='lower right', fontsize=12)	×
49
50	if show_error:	×
51	_, _, rms_dxdt = rmse(x, dt, x_hat, dxdt_hat, x_truth, dxdt_truth)	×
52	R_sqr = error_correlation(dxdt_hat, dxdt_truth)	×
53	axes[1].text(0.05, 0.95, f"RMSE = {rms_dxdt:.2f}\n$R^2$ = {R_sqr:.2g}",	×
54	transform=axes[1].transAxes, fontsize=15, verticalalignment='top')
55
56	return fig, axes	×
57
58
59	def plot_comparison(dt, dxdt_truth, dxdt_hat1, title1, dxdt_hat2, title2, dxdt_hat3, title3):	1✔
60	"""This is intended to show method performances with different choices of parameter"""
61	t = np.arange(0, dt*len(dxdt_truth), dt)	×
62	fig, axes = plt.subplots(1, 3, figsize=(22,6))	×
63
64	for i,(dxdt_hat,title) in enumerate(zip([dxdt_hat1, dxdt_hat2, dxdt_hat3], [title1, title2, title3])):	×
65	axes[i].plot(t, dxdt_truth, '--', color='black', linewidth=3, label=r"true $\frac{dx}{dt}$")	×
66	axes[i].plot(t, dxdt_hat, color='red', label=r"est. $\hat{\frac{dx}{dt}}$")	×
67	if i==0: axes[i].set_ylabel('Velocity', fontsize=18)	×
68	axes[i].set_xlabel('Time', fontsize=18)	×
69	axes[i].tick_params(axis='x', labelsize=15)	×
70	axes[i].tick_params(axis='y', labelsize=15)	×
71	axes[i].set_title(title, fontsize=18)	×
72	if i==2: axes[i].legend(loc='lower right', fontsize=12)	×
73	_, _, rms_dxdt = rmse(np.zeros(dxdt_hat.shape), dt, np.zeros(dxdt_hat.shape), dxdt_hat, np.zeros(dxdt_hat.shape), dxdt_truth)	×
74	R_sqr = error_correlation(dxdt_hat, dxdt_truth)	×
75	axes[i].text(0.05, 0.95, f"RMSE = {rms_dxdt:.2f}\n$R^2$ = {R_sqr:.2g}",	×
76	transform=axes[i].transAxes, fontsize=15, verticalalignment='top')
77
78	fig.tight_layout()	×
79
80
81	def rmse(x, dt, x_hat, dxdt_hat, x_truth=None, dxdt_truth=None, padding=0):	1✔
82	"""Evaluate x_hat based on RMSE, calculating different ones depending on whether :code:`dxdt_truth`
83	and :code:`x_truth` are known.
84
85	:param np.array[float] x: data that was differentiated
86	:param float dt: step size
87	:param np.array[float] x_hat: estimated (smoothed) x
88	:param np.array[float] dxdt_hat: estimated xdot
89	:param np.array[float] x_truth: true value of x, if known
90	:param np.array[float] dxdt_truth: true value of dxdt, if known, optional
91	:param int padding: number of snapshots on either side of the array to ignore when calculating
92	the metric. If :code:`'auto'`, defaults to 2.5% of the size of x
93
94	:return: tuple[float, float, float] containing\n
95	- rms_rec_x -- RMS error between the integral of dxdt_hat and x
96	- rms_x -- RMS error between x_hat and x_truth, returns None if x_truth is None
97	- rms_dxdt -- RMS error between dxdt_hat and dxdt_truth, returns None if dxdt_hat is None
98	"""
99	if np.isnan(x_hat).any():	×
100	return np.nan, np.nan, np.nan	×
101	if padding == 'auto':	×
102	padding = int(0.025*len(x))	×
103	padding = max(padding, 1)	×
104	s = slice(padding, len(x)-padding) # slice out data we want to measure	×
105
106	# RMS of dxdt and x_hat
107	root = np.sqrt(s.stop - s.start)	×
108	rms_dxdt = np.linalg.norm(dxdt_hat[s] - dxdt_truth[s]) / root if dxdt_truth is not None else None	×
109	rms_x = np.linalg.norm(x_hat[s] - x_truth[s]) / root if x_truth is not None else None	×
110
111	# RMS reconstructed x from integrating dxdt vs given raw x, available even in the absence of ground truth
112	rec_x = utility.integrate_dxdt_hat(dxdt_hat, dt)	×
113	x0 = utility.estimate_integration_constant(x, rec_x)	×
114	rec_x = rec_x + x0	×
115	rms_rec_x = np.linalg.norm(rec_x[s] - x[s]) / root	×
116
117	return rms_rec_x, rms_x, rms_dxdt	×
118
119
120	def error_correlation(dxdt_hat, dxdt_truth, padding=0):	1✔
121	"""Calculate the error correlation (pearsons correlation coefficient) between the estimated
122	dxdt and true dxdt
123
124	:param np.array[float] dxdt_hat: estimated xdot
125	:param np.array[float] dxdt_truth: true value of dxdt, if known, optional
126	:param int padding: number of snapshots on either side of the array to ignore when calculating
127	the metric. If :code:`'auto'`, defaults to 2.5% of the size of x
128
129	:return: (float) -- r-squared correlation coefficient
130	"""
131	if padding == 'auto':	×
132	padding = int(0.025*len(dxdt_hat))	×
133	padding = max(padding, 1)	×
134	s = slice(padding, len(dxdt_hat)-padding) # slice out data we want to measure	×
135	errors = (dxdt_hat[s] - dxdt_truth[s])	×
136	r = stats.linregress(dxdt_truth[s] - np.mean(dxdt_truth[s]), errors)	×
137	return r.rvalue**2	×
138
139
140	def total_variation(x, padding=0):	1✔
141	"""Calculate the total variation of an array. Used by optimizer.
142
143	:param np.array[float] x: data
144	:param int padding: number of snapshots on either side of the array to ignore when calculating
145	the metric. If :code:`'auto'`, defaults to 2.5% of the size of x
146
147	:return: (float) -- total variation
148	"""
149	if np.isnan(x).any():	×
150	return np.nan	×
151	if padding == 'auto':	×
152	padding = int(0.025*len(x))	×
153	padding = max(padding, 1)	×
154	x = x[padding:len(x)-padding]	×
155
156	return np.linalg.norm(x[1:]-x[:-1], 1)/len(x) # normalized version of what cvxpy.tv does	×

florisvb / PyNumDiff / 19090154521

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous