GA 1.4 setup

content/GA_1_4/data/groundwater_levels.csv

+,times,observations[mm]
+0,2017-01-01,-109.698
+1,2017-02-01,-102.044
+2,2017-03-01,-127.198
+3,2017-04-01,-151.899
+4,2017-05-01,-149.629
+5,2017-06-01,-132.828
+6,2017-07-01,-122.058
+7,2017-08-01,-142.505
+8,2017-09-01,-122.97
+9,2017-10-01,-123.722
+10,2017-11-01,-111.894
+11,2017-12-01,-118.856
+12,2018-01-01,-109.28
+13,2018-02-01,-114.459
+14,2018-03-01,-115.165
+15,2018-04-01,-139.411
+16,2018-05-01,-156.938
+17,2018-06-01,-166.784
+18,2018-07-01,-142.645
+19,2018-08-01,-133.115
+20,2018-09-01,-139.957
+21,2018-10-01,-112.534
+22,2018-11-01,-107.121
+23,2018-12-01,-122.64
+24,2019-01-01,-117.268

content/GA_1_4/data/insar_observations.csv

+,times,observations[m]
+0,2017-01-01,-0.0037915269917409627
+1,2017-01-13,-0.005998562574055986
+2,2017-01-25,-0.011402835952374462
+3,2017-02-06,-0.009919855648270475
+4,2017-02-18,-0.011282623963035473
+5,2017-03-02,-0.016512756482362813
+6,2017-03-14,-0.016996845493251048
+7,2017-03-26,-0.019696104362885405
+8,2017-04-07,-0.021837098696008004
+9,2017-04-19,-0.022895436362601166
+10,2017-05-01,-0.022327056485246743
+11,2017-05-13,-0.020821197338077402
+12,2017-05-25,-0.024863340904648194
+13,2017-06-06,-0.02056783048161778
+14,2017-06-18,-0.02213817613870604
+15,2017-06-30,-0.02163775448792241
+16,2017-07-12,-0.02560414369326818
+17,2017-07-24,-0.02473621201120934
+18,2017-08-05,-0.023681466237484187
+19,2017-08-17,-0.0282966692824492
+20,2017-08-29,-0.029358013071156923
+21,2017-09-10,-0.02852166622379215
+22,2017-09-22,-0.024896329017228425
+23,2017-10-04,-0.021029517374774372
+24,2017-10-16,-0.023354103774718086
+25,2017-10-28,-0.021855914891603143
+26,2017-11-09,-0.024765897955232773
+27,2017-11-21,-0.022823151059776652
+28,2017-12-03,-0.026800418671976987
+29,2017-12-15,-0.02509500935712605
+30,2017-12-27,-0.026505283952593273
+31,2018-01-08,-0.02670526176685112
+32,2018-01-20,-0.025821579299208672
+33,2018-02-01,-0.027508046803947844
+34,2018-02-13,-0.024288527653980718
+35,2018-02-25,-0.02670351434206529
+36,2018-03-09,-0.026471702793063708
+37,2018-03-21,-0.029256128085217777
+38,2018-04-02,-0.027596164768845584
+39,2018-04-14,-0.032995994476835454
+40,2018-04-26,-0.031334671938387675
+41,2018-05-08,-0.03447316193762846
+42,2018-05-20,-0.03646288082564127
+43,2018-06-01,-0.03622067189983248
+44,2018-06-13,-0.0373391550961804
+45,2018-06-25,-0.03482925880558608
+46,2018-07-07,-0.03375116486465379
+47,2018-07-19,-0.03153714180290244
+48,2018-07-31,-0.03053964391509817
+49,2018-08-12,-0.03146927023464485
+50,2018-08-24,-0.03159840378773475
+51,2018-09-05,-0.03141629804036172
+52,2018-09-17,-0.032855680003182834
+53,2018-09-29,-0.02881159629275797
+54,2018-10-11,-0.026992604435845863
+55,2018-10-23,-0.02767485725544704
+56,2018-11-04,-0.023014576182087923
+57,2018-11-16,-0.02589125901511123
+58,2018-11-28,-0.028883565745858694
+59,2018-12-10,-0.030314178558457035
+60,2018-12-22,-0.03027546561762635

content/GA_1_4/functions.py

+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.stats import norm
+import ipywidgets as widgets
+from ipywidgets import interact
+
+def plot_model(d, alt_model=None):
+    """Time Series of observations with model and CI.
+    
+    Uses dict as input defined from existing values in dict
+    Outputs figure and axes objects.
+
+    alt_model: tuple to add a line to the plot, 3 elements:
+      - string: array with the times
+      - list or ndarray with times (x-axis)
+      - list or ndarray with model values (y-axis)
+      - e.g., alt_model=('model 2', [...], [...])
+    """
+    
+    times = d['times']
+    y = d['y']
+    y_hat = d['y_hat']
+    CI_Y_hat = d['CI_Y_hat']
+    data_type = d['data_type']
+
+    fig, ax = plt.subplots(figsize = (15,5))
+
+    ax.plot(times, y, 'k+',  label = 'Observations')
+    ax.plot(times, y_hat,  label = 'Fitted model')
+    ax.fill_between(times,
+                    (y_hat - CI_Y_hat), 
+                    (y_hat + CI_Y_hat),
+                    facecolor='orange',
+                    alpha=0.4,
+                    label = f'{(1-d['alpha'])*100:.1f}%' + ' Confidence Region')
+    
+    if alt_model is not None:
+        ax.plot(alt_model[1],
+                alt_model[2],
+                label = alt_model[0])
+    
+    ax.legend()
+    ax.set(xlabel='Time',
+           ylabel=f'{data_type}' + ' Displacement [mm]',
+           title=f'{data_type}' + ' Observations and Fitted Model')
+
+    
+    return fig, ax
+
+
+def plot_residual(d):
+    """Time Series of residuals and CI.
+    
+    Uses dict as input defined from existing values in dict
+    Outputs figure and axes objects.
+    """
+    
+    times = d['times']
+    e_hat = d['e_hat']
+    CI_res = d['CI_res']
+    data_type = d['data_type']
+
+    fig, ax = plt.subplots(figsize = (15,5))
+
+    ax.plot(times, e_hat, 'o', markeredgecolor='black', label='Residual')
+    ax.plot(times, -CI_res,
+             '--', color='orange', 
+             label = f'{(1-d['alpha'])*100:.1f}%' + ' Confidence Region')
+    ax.plot(times, +CI_res,
+             '--', color='orange')
+    
+    ax.legend()
+    ax.set(xlabel='Time',
+           ylabel=f'{data_type}' + ' residual [mm]',
+           title=f'{data_type}' + ' Residuals')
+
+    return fig, ax
+
+
+def plot_residual_histogram(d):
+    """Histogram of residuals with Normal PDF.
+    
+    Uses dict as input defined from existing values in dict
+    Outputs figure and axes objects.
+    """
+    
+    e_hat = d['e_hat']
+    data_type = d['data_type']
+
+    fig, ax = plt.subplots(figsize = (7,5))
+    
+    ax.hist(e_hat, density=True,  edgecolor='black')
+    x = np.linspace(np.min(e_hat), np.max(e_hat), num=100);
+    ax.plot(x,
+            norm.pdf(x, loc=0.0, scale = np.std(e_hat)),
+            linewidth=4.0)
+
+    ax.set(xlabel='Residuals [mm]',
+           ylabel=f'{data_type}' + ' Density [-]',
+           title=f'{data_type}' + ' Residuals Histogram')
+    
+    print ('The mean value of the', data_type, 'residuals is',
+           np.around(np.mean(e_hat),5), 'mm')
+    print ('The standard deviation of the', data_type, 'residuals is',
+        np.around(np.std(e_hat),3), 'mm')
+
+
+    return fig, ax
+
+def xhat_slider_plot(A, y, t, Sigma_y=None):
+    """Interactive plot of the solution space for x_hat."""
+    n, k = A.shape
+
+    if Sigma_y is None:
+        Sigma_y = np.eye(n)
+    xhat = np.linalg.inv(A.T @ np.linalg.inv(Sigma_y) @ A) @ A.T @ np.linalg.inv(Sigma_y) @ y
+    Sigma_xhat = np.linalg.inv(A.T @ np.linalg.inv(Sigma_y) @ A)
+    std_xhat = np.sqrt(np.diag(Sigma_xhat))
+
+    sliders = {}
+    for i in range(k):
+        sliders[f'xhat_{i}'] = widgets.FloatSlider(value=xhat[i],
+                                                   min=xhat[i] - 10*std_xhat[i],
+                                                   max=xhat[i] + 10*std_xhat[i],
+                                                   step=0.1*std_xhat[i],
+                                                   description=f'xhat_{i}')
+
+    def update_plot(**kwargs):
+        xhat_values = np.array([kwargs[f'xhat_{i}'] for i in range(k)])
+        y_fit = A @ xhat_values
+        W = np.linalg.inv(Sigma_y)
+        ss_res = (y - y_fit).T @ W @ (y - y_fit)
+
+        plt.figure(figsize=(10, 5))
+        plt.plot(t, y, 'o', label='data')
+        plt.plot(t, y_fit, label='y_fit')
+        plt.title(f'Mean of squared residuals: {ss_res:.2f}')
+        plt.ylabel('y')
+        plt.xlabel('t')
+        plt.grid()
+        plt.legend()
+        plt.show()
+
+    interact(update_plot, **sliders)
+
+# Example usage
+# A, y, t should be defined before calling this function
+# xhat_slider_plot(A, y, t)
\ No newline at end of file