diff --git a/content/GA_1_4/Analysis_solution.ipynb b/content/GA_1_4/Analysis_solution.ipynb
index 7815eac8a8e2d17c1afa2f246ce200f225f2198f..f5fbf7b428055dab440c367eca9e036bdc77982a 100644
--- a/content/GA_1_4/Analysis_solution.ipynb
+++ b/content/GA_1_4/Analysis_solution.ipynb
@@ -50,7 +50,7 @@
},
{
"cell_type": "code",
- "execution_count": 25,
+ "execution_count": 1,
"id": "b381f517",
"metadata": {},
"outputs": [],
@@ -97,7 +97,7 @@
},
{
"cell_type": "code",
- "execution_count": 26,
+ "execution_count": 2,
"id": "b8be3730",
"metadata": {},
"outputs": [],
@@ -116,7 +116,7 @@
},
{
"cell_type": "code",
- "execution_count": 27,
+ "execution_count": 3,
"id": "e02d5566",
"metadata": {},
"outputs": [
@@ -154,6 +154,47 @@
"model_summary(m2_blue)"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "cabfa723",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "data_type\n",
+ "model_type\n",
+ "times\n",
+ "y\n",
+ "days\n",
+ "groundwater\n",
+ "groundwater_data\n",
+ "A\n",
+ "std\n",
+ "Sigma_Y\n",
+ "Sigma_X_hat\n",
+ "x_hat\n",
+ "y_hat\n",
+ "e_hat\n",
+ "Sigma_Y_hat\n",
+ "std_y\n",
+ "Sigma_e_hat\n",
+ "std_e_hat\n",
+ "k\n",
+ "CI_Y\n",
+ "CI_res\n",
+ "CI_Y_hat\n",
+ "alpha\n"
+ ]
+ }
+ ],
+ "source": [
+ "for key in m1_blue.keys():\n",
+ " print(key)"
+ ]
+ },
{
"cell_type": "markdown",
"id": "0b4e08a3",
@@ -171,14 +212,14 @@
},
{
"cell_type": "code",
- "execution_count": 28,
+ "execution_count": 5,
"id": "73af3940",
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
- "model_id": "3d28d06251d7448881d0883ed9e54c8a",
+ "model_id": "f86d3002a0d346f3a55f55c0fe788062",
"version_major": 2,
"version_minor": 0
},
@@ -195,7 +236,8 @@
"x1_slider = widgets.FloatSlider(value=0, min=-0.1, max=0.1, step=0.001, description='x1')\n",
"x2_slider = widgets.FloatSlider(value=1, min=-1, max=1, step=0.01, description='x2')\n",
"\n",
- "interact(update_plot, x0=x0_slider, x1=x1_slider, x2=x2_slider, x3=fixed(None),\n",
+ "interact(model_widget,\n",
+ " x0=x0_slider, x1=x1_slider, x2=x2_slider, x3=fixed(None),\n",
" m=[('Data Type 1', m1_blue), ('Data Type 2', m2_blue)]);\n"
]
},
@@ -234,7 +276,7 @@
},
{
"cell_type": "code",
- "execution_count": 29,
+ "execution_count": 6,
"id": "712e4114",
"metadata": {},
"outputs": [],
@@ -248,6 +290,7 @@
" d['Sigma_Y'] = d_old['Sigma_Y']\n",
" d['days'] = d_old['days']\n",
" d['groundwater'] = d_old['groundwater']\n",
+ " d['groundwater_data'] = d_old['groundwater_data']\n",
" return d\n",
"\n",
"m1 = initialize_new_dict(m1_blue)\n",
@@ -255,6 +298,60 @@
"\n"
]
},
+ {
+ "cell_type": "markdown",
+ "id": "dee36d36",
+ "metadata": {},
+ "source": [
+ "<div style=\"background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
+ "<p>\n",
+ "<b>Task 0.4: </b> \n",
+ " \n",
+ "Confirm that the stochastic model is transferred properly by printing the appropriate key-value pair in the dictionary.\n",
+ " \n",
+ "</p>\n",
+ "</div>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "id": "9aeda35f",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[4. 0. 0. ... 0. 0. 0.]\n",
+ " [0. 4. 0. ... 0. 0. 0.]\n",
+ " [0. 0. 4. ... 0. 0. 0.]\n",
+ " ...\n",
+ " [0. 0. 0. ... 4. 0. 0.]\n",
+ " [0. 0. 0. ... 0. 4. 0.]\n",
+ " [0. 0. 0. ... 0. 0. 4.]]\n",
+ "[[225. 0. 0. ... 0. 0. 0.]\n",
+ " [ 0. 225. 0. ... 0. 0. 0.]\n",
+ " [ 0. 0. 225. ... 0. 0. 0.]\n",
+ " ...\n",
+ " [ 0. 0. 0. ... 225. 0. 0.]\n",
+ " [ 0. 0. 0. ... 0. 225. 0.]\n",
+ " [ 0. 0. 0. ... 0. 0. 225.]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "# YOUR_CODE_HERE\n",
+ "\n",
+ "# SOLUTION:\n",
+ "# First uncomment and run this to quickly view the keys:\n",
+ "# for key in m1.keys():\n",
+ "# print(key)\n",
+ "# Clearly Sigma_Y is what we need:\n",
+ "print(m1['Sigma_Y'])\n",
+ "print(m2['Sigma_Y'])\n"
+ ]
+ },
{
"cell_type": "markdown",
"id": "80a9b60f",
@@ -281,29 +378,35 @@
"4. $k$, as 'groundwater factor' that can be seen as the response of the layers due to changes in the groundwater level.\n",
"\n",
"It can be seen that the model is non-linear. We will therefore use non-linear least-squares to solve for the unknown parameters. \n",
- "\n"
+ "\n",
+ "### Functional Model: Python Implementation\n",
+ "\n",
+ "We will use a Python function to evaluate the non-linear function. Recall that the general form of the functional model is $\\mathrm{Y}=q(\\mathrm{X})$, where $\\mathrm{X}$ represents the vector of model parameters. Note also that to implement this model in practice often additional (deterministic!) parameters are required. This can be written in a Python function with the following form of input arguments:\n",
+ "\n",
+ "```\n",
+ "Y = functional_model(X, <auxiliary_arguments>)\n",
+ "```\n",
+ "\n",
+ "Where `<auxiliary_arguments>` will be different in type and/or number on a case-by-case basis. Your code will generally be more compact and adaptable to other cases if the parameters are specified in a list, array or tuple. "
]
},
{
"cell_type": "markdown",
- "id": "6653c958-0686-448f-8ecc-f21946db82bf",
+ "id": "e28b32ed",
"metadata": {},
"source": [
"<div style=\"background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
"<p>\n",
- "<b>Task 1.1: </b> \n",
- " \n",
- "Choose initial values for the model parameters. Use the code and Markdown cells below to justify your decision. We suggest two possible approaches: a) use the forward model and make a plot to see if you can get it in the right order of magnitude, or b) make an inference about what the values might be using knowledge about each term in the model.\n",
- " \n",
- "<i>Note: it may be useful at this point to define a function for your forward model to check the values. You will be able to re-use it in later Tasks as well.<i>\n",
- " \n",
+ "<b>Task 1.1:</b> \n",
+ "\n",
+ "Read the code below to understand the functional model (and fill in the missing parameter values). Be sure to understand which objects are the parameters of interest, and what the auxiliary arguments are (type and number).\n",
"</p>\n",
"</div>"
]
},
{
"cell_type": "code",
- "execution_count": 30,
+ "execution_count": 7,
"id": "3ba5bf67",
"metadata": {
"id": "3ba5bf67",
@@ -311,31 +414,124 @@
},
"outputs": [],
"source": [
- "def forward_model(time, do_i, R_i, a_i, k_i, GW):\n",
- " \"\"\"Compute the displacements based on our initial guess.\n",
+ "def functional_model(X, d):\n",
+ " \"\"\"Functional model, q: ground surface displacement.\n",
+ "\n",
+ " Inputs:\n",
+ " X: tuple, list or array of parameters\n",
+ " (d, R, a, k)\n",
+ " d: dictionary of model parameters\n",
+ "\n",
+ " Outputs: ndarray of ground level\n",
+ "\n",
+ " Note: \"times\" is not included as an argument because\n",
+ " the model is only configured to compute the functional\n",
+ " model at times where the groundwater measurements have\n",
+ " been interpolated (i.e., the times of the satellite\n",
+ " observations). To use this model for prediction, the\n",
+ " interpolation function must be incorporated in the\n",
+ " dictionary and/or the function. \n",
" \"\"\"\n",
" \n",
- " y_comp = do_i + R_i*(1-np.exp(-time/a_i))+k_i*GW\n",
+ " # Hint: use d['days'] and d['groundwater'] for\n",
+ " # the deterministic parameters\n",
+ "\n",
+ " # y = (YOUR_CODE_HERE\n",
+ " # + YOUR_CODE_HERE*(1 - np.exp(-d['days']/YOUR_CODE_HERE))\n",
+ " # + YOUR_CODE_HEREX[3]*d['groundwater']\n",
+ " # )\n",
+ " \n",
+ " # SOLUTION:\n",
+ " y = (X[0]\n",
+ " + X[1]*(1 - np.exp(-d['days']/X[2]))\n",
+ " + X[3]*d['groundwater']\n",
+ " )\n",
" \n",
- " return y_comp"
+ " return y\n",
+ "\n",
+ "m1['functional_model'] = functional_model\n",
+ "m2['functional_model'] = functional_model"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fbd55e2b",
+ "metadata": {},
+ "source": [
+ "<div style=\"background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
+ "<p>\n",
+ "<b>Task 1.2: </b> \n",
+ " \n",
+ "Use the widget below to gain understanding of the non-linear model and its parameters, as well as the role that the deterministic parameters (observation times and groundwater measurements) play in the functional model.\n",
+ " \n",
+ "</p>\n",
+ "</div>"
]
},
{
"cell_type": "code",
- "execution_count": 31,
+ "execution_count": 8,
+ "id": "b20c1e61",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "9fd705b21f954ab2bde969ecf1c7267e",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "interactive(children=(FloatSlider(value=0.0, description='x0', max=10.0, min=-10.0), FloatSlider(value=0.0, de…"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x0_slider = widgets.FloatSlider(value=0, min=-10, max=10, step=0.1, description='x0')\n",
+ "x1_slider = widgets.FloatSlider(value=0, min=-50, max=50, step=1, description='x1')\n",
+ "x2_slider = widgets.FloatSlider(value=1, min=10, max=1000, step=10, description='x2')\n",
+ "x3_slider = widgets.FloatSlider(value=1, min=-1, max=1, step=0.01, description='x3')\n",
+ "\n",
+ "interact(model_widget,\n",
+ " x0=x0_slider, x1=x1_slider, x2=x2_slider, x3=x3_slider,\n",
+ " m=[('Data Type 1', m1), ('Data Type 2', m2)]);\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c9bfc4d7",
+ "metadata": {},
+ "source": [
+ "<div style=\"background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
+ "<p>\n",
+ "<b>Task 1.3: </b> \n",
+ " \n",
+ "Choose initial values for the model parameters. Use the code and Markdown cells below to justify your decision. We suggest two possible approaches: a) use the functional model and make a plot to see if you can get it in the right order of magnitude, or b) make an inference about what the values might be using knowledge about each term in the model.\n",
+ " \n",
+ "</p>\n",
+ "</div>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
"id": "66ddf324-11fe-4516-a37c-b8a90ce84c18",
"metadata": {},
"outputs": [],
"source": [
- "# k_i = \n",
- "# R_i = \n",
- "# a_i = \n",
- "# do_i = \n",
- "# SOLUTION\n",
- "k_i = 0.15\n",
+ "# d_i = YOUR_CODE_HERE\n",
+ "# R_i = YOUR_CODE_HERE\n",
+ "# a_i = YOUR_CODE_HERE\n",
+ "# k_i = YOUR_CODE_HERE\n",
+ "\n",
+ "# SOLUTION (values don't need to be exactly the same)\n",
+ "d_i = 9\n",
"R_i = -25\n",
"a_i = 300\n",
- "do_i = 9"
+ "k_i = 0.15"
]
},
{
@@ -373,7 +569,7 @@
"<p>\n",
"<b>Task 8.2: </b> \n",
" \n",
- "Set up the Jacobian matrix for the non-linear least-squares. Using the initial values from Task 8.1, print the first 5 rows of the Jacobian matrix to confirm that your function works correctly (i.e., confirm that code runs and initial values give acceptable results).\n",
+ "Set up the Jacobian matrix for the non-linear least-squares by completing the function below.\n",
" \n",
"\n",
"</p>\n",
@@ -382,7 +578,7 @@
},
{
"cell_type": "code",
- "execution_count": 32,
+ "execution_count": 21,
"id": "4aef41e9",
"metadata": {
"id": "4aef41e9",
@@ -390,32 +586,69 @@
},
"outputs": [],
"source": [
- "def jacobian(time, do_i, R_i, a_i, k_i, GW):\n",
- " \"\"\"Create the Jacobian matrix.\n",
+ "def jacobian(X, d):\n",
+ " \"\"\"Compute Jacobian of the functional model.\n",
+ "\n",
+ " Functional model, q: ground surface displacement.\n",
+ "\n",
+ " Inputs:\n",
+ " X: tuple, list or array of parameters\n",
+ " (d, R, a, k)\n",
+ " d: dictionary of model parameters\n",
+ "\n",
+ " Outputs: The Jacobian\n",
+ " (partial derivatives w.r.t. parameters)\n",
+ "\n",
+ " Note: \"times\" is not included as an argument because\n",
+ " the model is only configured to compute the functional\n",
+ " model at times where the groundwater measurements have\n",
+ " been interpolated (i.e., the times of the satellite\n",
+ " observations). To use this model for prediction, the\n",
+ " interpolation function must be incorporated in the\n",
+ " dictionary and/or the function. \n",
+ " \"\"\"\n",
" \n",
- " The columns represent the linearized system of equations.\n",
- " \"\"\" \n",
+ " # Hint: use d['days'] and d['groundwater'] for\n",
+ " # the deterministic parameters\n",
+ "\n",
+ "\n",
+ " # J1 = \n",
+ " # J2 = \n",
+ " # J3 = \n",
+ " # J4 = \n",
+ " # j = YOUR_CODE_HERE\n",
+ "\n",
+ " # SOLUTION\n",
+ " J1 = np.ones(len(d['days']))\n",
+ " J2 = 1 - np.exp(-d['days']/X[2])\n",
+ " J3 = -X[1]*d['days']/X[2]**2 * np.exp(-d['days']/X[2])\n",
+ " J4 = np.ones(len(d['days']))*d['groundwater']\n",
+ " J = np.column_stack((J1, J2, J3, J4))\n",
" \n",
- " # J_c1 = \n",
- " # J_c2 = \n",
- " # J_c3 = \n",
- " # J_c4 = \n",
+ " return J"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "91709ccb",
+ "metadata": {},
+ "source": [
+ "<div style=\"background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
+ "<p>\n",
+ "<b>Task 8.2: </b> \n",
" \n",
- " # SOLUTION\n",
- " J_c1 = np.ones(len(time))\n",
- " J_c2 = 1 - np.exp(-time/a_i)\n",
- " J_c3 = -R_i*time/a_i**2 * np.exp(-time/a_i)\n",
- " J_c4 = np.ones(len(time))*GW\n",
+ "Confirm that the Jacobian function works properly by using it with the initial parameter values established above. If the function is called successfully, check the result by printing the first few rows of the Jacobian matrix and checking the dimensions.\n",
"\n",
+ "<em>Tip: it would be useful for the Report if you also printed the dimensions and/or redundancy of the model.</em>\n",
" \n",
- " J = np.column_stack((J_c1, J_c2, J_c3, J_c4))\n",
- " \n",
- " return J"
+ "\n",
+ "</p>\n",
+ "</div>"
]
},
{
"cell_type": "code",
- "execution_count": 33,
+ "execution_count": 22,
"id": "c57e5e26",
"metadata": {
"colab": {
@@ -440,7 +673,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
- "The first 5 rows of the Jacobian matrix for InSAR are:\n",
+ "The first 5 rows of the Jacobian matrix (INSAR):\n",
"[[ 1.000e+00 0.000e+00 0.000e+00 -1.097e+02]\n",
" [ 1.000e+00 3.921e-02 3.203e-03 -1.067e+02]\n",
" [ 1.000e+00 7.688e-02 6.154e-03 -1.038e+02]\n",
@@ -454,55 +687,20 @@
}
],
"source": [
- "Jacob_insar = jacobian(m1['days'], do_i, R_i, a_i, k_i, m1['groundwater'])\n",
+ "# YOUR_CODE_HERE\n",
"\n",
- "print ('The first 5 rows of the Jacobian matrix for InSAR are:')\n",
- "print (Jacob_insar[0:5,:])\n",
+ "# SOLUTION:\n",
+ "test_J = jacobian((d_i, R_i, a_i, k_i), m1)\n",
"\n",
- "n_2 = np.shape(Jacob_insar)[1]\n",
+ "print ('The first 5 rows of the Jacobian matrix (INSAR):')\n",
+ "print (test_J[0:5,:])\n",
+ "\n",
+ "n_2 = np.shape(test_J)[1]\n",
"print(f'\\nThe number of unknowns is {n_2}')\n",
"print(f'The redundancy (InSAR) is {m1[\"y\"].shape[0] - n_2}')\n",
"print(f'The redundancy (GNSS) is {m2[\"y\"].shape[0] - n_2}')"
]
},
- {
- "cell_type": "markdown",
- "id": "5a640428",
- "metadata": {},
- "source": [
- "<div style=\"background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
- "<p>\n",
- "<b>Task 8.3:</b> \n",
- " \n",
- "What is the redundancy of the model with InSAR and GNSS, respectively?\n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "468b9145-b05a-4c06-90be-dbd85145dd12",
- "metadata": {},
- "source": [
- "**Write your answer in this Markdown cell.**"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "193256db",
- "metadata": {},
- "source": [
- "<div style=\"background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
- "<p>\n",
- "\n",
- "$\\textbf{Solution}$\n",
- "\n",
- "- The number of unknowns is 4\n",
- "- The redundancy (InSAR) is 57\n",
- "- The redundancy (GNSS) is 726\n",
- " \n"
- ]
- },
{
"cell_type": "markdown",
"id": "b4633d84",
@@ -523,16 +721,23 @@
"<p>\n",
"<b>Task 9: </b> \n",
" \n",
- "Set up a Gauss-Newton iteration algorithm (complete the code below). Choose the criterion to stop the iteration. \n",
- "Explain below how and why you define the stop-criterium. \n",
+ "Set up a Gauss-Newton iteration algorithm by completing the function below.\n",
"\n",
"</p>\n",
"</div>"
]
},
+ {
+ "cell_type": "markdown",
+ "id": "ef9e4ad4",
+ "metadata": {},
+ "source": [
+ "<div style=\"background-color:#facb8e; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%\"> <p>Note that the iteration scheme reuses the function <code>BLUE()</code> from last week, which has been included in the files <code>functions.py</code> (imported at the top of this notebook). In addition, note that during the iteration scheme, specific values of the dictionary are defined based on the the _difference_ between the inital guess and the computed value. This creates a _linear_ system of equations and facilitates the use BLUE, which otherwise can only solve for <em>linear</em> functional models.</p></div>"
+ ]
+ },
{
"cell_type": "code",
- "execution_count": 34,
+ "execution_count": 12,
"id": "FbsD75dR0w5b",
"metadata": {
"id": "FbsD75dR0w5b",
@@ -540,31 +745,28 @@
},
"outputs": [],
"source": [
- "def nlsq_iteration(x0, d):\n",
- " \"\"\"Use Gauss-Newton iteration to find non-linear parameters.\"\"\"\n",
- "\n",
+ "def gauss_newton_iteration(X0, d):\n",
+ " \"\"\"Use Gauss-Newton iteration to find non-linear parameters.\n",
" \n",
- " day = d['days']\n",
- " y_obs = d['y']\n",
- " Sigma_Y = d['Sigma_Y']\n",
- " GW = d['groundwater']\n",
+ " Inputs:\n",
+ " X0: initial guess for the parameters (d, R, a, k)\n",
+ " d: dictionary of model parameters\n",
"\n",
- " xnrm = 1000 # initialize stop criteria\n",
+ " Outputs: dictionary with the non-linear model results.\n",
+ " \"\"\"\n",
"\n",
- " xhat_i = np.zeros((50, 4))\n",
- " xhat_i[0,:] = x0\n",
+ " x_norm = 1000 # initialize stop criteria\n",
"\n",
- " do_i = xhat_i[0,0]\n",
- " R_i = xhat_i[0,1]\n",
- " a_i = xhat_i[0,2]\n",
- " k_i = xhat_i[0,3]\n",
+ " x_hat_i = np.zeros((50, 4))\n",
+ " x_hat_i[0,:] = X0\n",
+ " X_i = X0\n",
"\n",
- " iter_ind = 0\n",
+ " iteration = 0\n",
"\n",
- " while xnrm >= 1e-12 and iter_ind < 49:\n",
+ " while x_norm >= 1e-12 and iteration < 49:\n",
"\n",
" # computed deformation yi based on 'estimates' \n",
- " y_i = forward_model(day, do_i, R_i, a_i, k_i, GW)\n",
+ " y_i = functional_model(X_i, d)\n",
" \n",
" # dy = YOUR_CODE_HERE\n",
" # SOLTUION\n",
@@ -573,7 +775,7 @@
"\n",
" # J = YOUR_CODE_HERE\n",
" # SOLTUION\n",
- " J = jacobian(day, do_i, R_i, a_i, k_i, GW)\n",
+ " J = jacobian(X_i, d)\n",
" d['A'] = J\n",
"\n",
"\n",
@@ -583,31 +785,24 @@
" d = BLUE(d)\n",
"\n",
"\n",
- " xhat_i[iter_ind+1,:] = xhat_i[iter_ind,:] + d['x_hat'].T\n",
+ " x_hat_i[iteration+1,:] = x_hat_i[iteration,:] + d['x_hat'].T\n",
"\n",
- " do_i = xhat_i[iter_ind+1,0]\n",
- " R_i = xhat_i[iter_ind+1,1]\n",
- " a_i = xhat_i[iter_ind+1,2]\n",
- " k_i = xhat_i[iter_ind+1,3]\n",
- " \n",
+ " X_i = x_hat_i[iteration+1,:]\n",
"\n",
- " # xnrm = YOUR_CODE_HERE\n",
+ " # x_norm = YOUR_CODE_HERE\n",
" # SOLUTION\n",
- " xnrm = d['x_hat'].T @ np.linalg.inv(d['Sigma_X_hat']) @ d['x_hat']\n",
+ " x_norm = d['x_hat'].T @ np.linalg.inv(d['Sigma_X_hat']) @ d['x_hat']\n",
"\n",
" # Update the iteration number\n",
- " iter_ind += 1\n",
+ " iteration += 1\n",
"\n",
- " if iter_ind==49:\n",
+ " if iteration==49:\n",
" print(\"Number of iterations too large, check initial values.\")\n",
"\n",
- " xhat = xhat_i[iter_ind,:]\n",
- " xhat_i = xhat_i[0:iter_ind+1, :]\n",
- "\n",
- " d['xhat_i'] = xhat_i\n",
+ " d['x_hat_i'] = x_hat_i[0:iteration+1, :]\n",
" d['dx'] = d['x_hat']\n",
- " d['x_hat'] = d['xhat_i'][iter_ind,:]\n",
- " d['niter'] = iter_ind\n",
+ " d['x_hat'] = d['x_hat_i'][iteration,:]\n",
+ " d['niter'] = iteration\n",
" d['dy'] = d['y']\n",
" d['y'] = y_obs\n",
"\n",
@@ -669,58 +864,40 @@
},
{
"cell_type": "code",
- "execution_count": 35,
+ "execution_count": 19,
"id": "313380a0",
"metadata": {
"tags": []
},
"outputs": [
{
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "\n",
- " InSAR Reults for each iteration (#Interations = 7 )\n",
- "[[ 9.000e+00 -2.500e+01 3.000e+02 1.500e-01]\n",
- " [ 1.273e+01 -1.982e+01 1.407e+02 1.726e-01]\n",
- " [ 1.296e+01 -2.162e+01 1.768e+02 1.719e-01]\n",
- " [ 1.297e+01 -2.192e+01 1.793e+02 1.713e-01]\n",
- " [ 1.297e+01 -2.192e+01 1.794e+02 1.713e-01]\n",
- " [ 1.297e+01 -2.192e+01 1.794e+02 1.713e-01]\n",
- " [ 1.297e+01 -2.192e+01 1.794e+02 1.713e-01]\n",
- " [ 1.297e+01 -2.192e+01 1.794e+02 1.713e-01]]\n",
- "\n",
- " GNSS Reults for each iteration (#Interations = 10 )\n",
- "[[ 9.000e+00 -2.500e+01 3.000e+02 1.500e-01]\n",
- " [ 3.747e+00 -1.777e+01 2.435e+02 1.425e-01]\n",
- " [ 3.895e+00 -1.812e+01 2.269e+02 1.418e-01]\n",
- " [ 3.939e+00 -1.815e+01 2.247e+02 1.416e-01]\n",
- " [ 3.945e+00 -1.815e+01 2.242e+02 1.415e-01]\n",
- " [ 3.946e+00 -1.815e+01 2.241e+02 1.415e-01]\n",
- " [ 3.946e+00 -1.815e+01 2.241e+02 1.415e-01]\n",
- " [ 3.946e+00 -1.815e+01 2.241e+02 1.415e-01]\n",
- " [ 3.946e+00 -1.815e+01 2.241e+02 1.415e-01]\n",
- " [ 3.946e+00 -1.815e+01 2.241e+02 1.415e-01]\n",
- " [ 3.946e+00 -1.815e+01 2.241e+02 1.415e-01]]\n"
+ "ename": "NameError",
+ "evalue": "name 'd_i' is not defined",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
+ "\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)",
+ "Cell \u001b[1;32mIn[19], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m x0 \u001b[38;5;241m=\u001b[39m [d_i, R_i, a_i,k_i]\n\u001b[0;32m 3\u001b[0m m1 \u001b[38;5;241m=\u001b[39m nlsq_iteration(x0, m1)\n\u001b[0;32m 4\u001b[0m m2 \u001b[38;5;241m=\u001b[39m nlsq_iteration(x0, m2)\n",
+ "\u001b[1;31mNameError\u001b[0m: name 'd_i' is not defined"
]
}
],
"source": [
- "x0 = [do_i, R_i, a_i,k_i]\n",
+ "x0 = [d_i, R_i, a_i,k_i]\n",
"\n",
"m1 = nlsq_iteration(x0, m1)\n",
"m2 = nlsq_iteration(x0, m2)\n",
"\n",
"print('\\n InSAR Reults for each iteration (#Interations =', m1['niter'], ')')\n",
- "print(m1['xhat_i'])\n",
+ "print(m1['x_hat_i'])\n",
"\n",
"print('\\n GNSS Reults for each iteration (#Interations =', m2['niter'], ')')\n",
- "print(m2['xhat_i'])"
+ "print(m2['x_hat_i'])"
]
},
{
"cell_type": "code",
- "execution_count": 36,
+ "execution_count": 12,
"id": "7d52c8b9",
"metadata": {
"id": "7d52c8b9",
@@ -728,7 +905,7 @@
},
"outputs": [],
"source": [
- "# def plot_fit_iteration(xhat_i):\n",
+ "# def plot_fit_iteration(x_hat_i):\n",
"# \"\"\"Plot value of each parameter, each iteration.\"\"\"\n",
"# plt.figure(figsize = (15,4))\n",
"# plt.subplots_adjust(top = 2)\n",
@@ -745,31 +922,31 @@
"# plt.subplot(2,2,4)\n",
"# YOUR_CODE_HERE\n",
" \n",
- "def plot_fit_iteration(xhat_i):\n",
+ "def plot_fit_iteration(x_hat_i):\n",
" \"\"\"Plot value of each parameter, each iteration.\"\"\"\n",
" plt.figure(figsize = (15,4))\n",
" plt.subplots_adjust(top = 2)\n",
"\n",
" plt.subplot(2,2,1)\n",
- " plt.plot(xhat_i[:,0], linewidth=4)\n",
+ " plt.plot(x_hat_i[:,0], linewidth=4)\n",
" plt.title('Estimated offset')\n",
" plt.ylabel('Offset [mm]')\n",
" plt.xlabel('Number of iterations [-]')\n",
"\n",
" plt.subplot(2,2,2)\n",
- " plt.plot(xhat_i[:,1], linewidth=4)\n",
+ " plt.plot(x_hat_i[:,1], linewidth=4)\n",
" plt.title('Estimated R value')\n",
" plt.ylabel('Estimated R value [mm]')\n",
" plt.xlabel('Number of iterations [-]')\n",
"\n",
" plt.subplot(2,2,3)\n",
- " plt.plot(xhat_i[:,2], linewidth=4)\n",
+ " plt.plot(x_hat_i[:,2], linewidth=4)\n",
" plt.title('Estimated $a$ value')\n",
" plt.ylabel('a value [days]')\n",
" plt.xlabel('Number of iterations [-]')\n",
"\n",
" plt.subplot(2,2,4)\n",
- " plt.plot(xhat_i[:,3], linewidth=4)\n",
+ " plt.plot(x_hat_i[:,3], linewidth=4)\n",
" plt.title('Estimated GW factor')\n",
" plt.ylabel('Estimated GW factor [-]')\n",
" plt.xlabel('Number of iterations [-]')"
@@ -777,7 +954,7 @@
},
{
"cell_type": "code",
- "execution_count": 37,
+ "execution_count": 13,
"id": "OtSGy-5NfSP6",
"metadata": {
"colab": {
@@ -811,12 +988,12 @@
}
],
"source": [
- "plot_fit_iteration(m1['xhat_i'])"
+ "plot_fit_iteration(m1['x_hat_i'])"
]
},
{
"cell_type": "code",
- "execution_count": 38,
+ "execution_count": 14,
"id": "aELTXpENfvx_",
"metadata": {
"colab": {
@@ -850,7 +1027,7 @@
}
],
"source": [
- "plot_fit_iteration(m2['xhat_i'])"
+ "plot_fit_iteration(m2['x_hat_i'])"
]
},
{
@@ -923,7 +1100,7 @@
},
{
"cell_type": "code",
- "execution_count": 39,
+ "execution_count": 15,
"id": "5af9e513",
"metadata": {
"id": "5af9e513",
@@ -948,7 +1125,7 @@
},
{
"cell_type": "code",
- "execution_count": 40,
+ "execution_count": 16,
"id": "fJA-BwcchWga",
"metadata": {
"colab": {
@@ -1061,7 +1238,7 @@
},
{
"cell_type": "code",
- "execution_count": 41,
+ "execution_count": 17,
"id": "pDJfjxgs6veD",
"metadata": {
"colab": {
@@ -1090,7 +1267,7 @@
" <Axes: title={'center': 'InSAR Residuals'}, xlabel='Time', ylabel='InSAR residual [mm]'>)"
]
},
- "execution_count": 41,
+ "execution_count": 17,
"metadata": {},
"output_type": "execute_result"
},
@@ -1116,7 +1293,7 @@
}
],
"source": [
- "m1['y_hat'] = forward_model(m1['days'], m1['x_hat'][0], m1['x_hat'][1], m1['x_hat'][2],\\\n",
+ "m1['y_hat'] = functional_model(m1['days'], m1['x_hat'][0], m1['x_hat'][1], m1['x_hat'][2],\\\n",
" m1['x_hat'][3], m1['groundwater'])\n",
"m1 = get_CI(m1, 0.04)\n",
"plot_model(m1)\n",
@@ -1125,7 +1302,7 @@
},
{
"cell_type": "code",
- "execution_count": 42,
+ "execution_count": 18,
"id": "ngZMBQM87QMr",
"metadata": {
"colab": {
@@ -1154,7 +1331,7 @@
" <Axes: title={'center': 'GNSS Residuals'}, xlabel='Time', ylabel='GNSS residual [mm]'>)"
]
},
- "execution_count": 42,
+ "execution_count": 18,
"metadata": {},
"output_type": "execute_result"
},
@@ -1180,7 +1357,7 @@
}
],
"source": [
- "m2['y_hat'] = forward_model(m2['days'], m2['x_hat'][0], m2['x_hat'][1], m2['x_hat'][2],\\\n",
+ "m2['y_hat'] = functional_model(m2['days'], m2['x_hat'][0], m2['x_hat'][1], m2['x_hat'][2],\\\n",
" m2['x_hat'][3], m2['groundwater'])\n",
"m2 = get_CI(m2, 0.04)\n",
"plot_model(m2)\n",
@@ -1336,7 +1513,7 @@
},
{
"cell_type": "code",
- "execution_count": 43,
+ "execution_count": 19,
"id": "d431d4f5",
"metadata": {
"tags": []
@@ -1361,7 +1538,7 @@
},
{
"cell_type": "code",
- "execution_count": 44,
+ "execution_count": 20,
"id": "92b5b797-093d-46fb-a2c0-2b9ccdae4052",
"metadata": {
"colab": {
@@ -1399,7 +1576,7 @@
},
{
"cell_type": "code",
- "execution_count": 45,
+ "execution_count": 21,
"id": "ae6038da",
"metadata": {
"tags": []
diff --git a/content/GA_1_4/auxiliary_files/m1_blue.pickle b/content/GA_1_4/auxiliary_files/m1_blue.pickle
index 9bdca1a2772a719aa5bc0b0b5ad8902a454cf4d3..4ec71533197e9c7d843f44706da8b5f5e4f8778a 100644
Binary files a/content/GA_1_4/auxiliary_files/m1_blue.pickle and b/content/GA_1_4/auxiliary_files/m1_blue.pickle differ
diff --git a/content/GA_1_4/auxiliary_files/m2_blue.pickle b/content/GA_1_4/auxiliary_files/m2_blue.pickle
index dc58bd13d51c06bb466d4d7238682c92c475db76..c672521e5203c355c21855153c531903ff064fc4 100644
Binary files a/content/GA_1_4/auxiliary_files/m2_blue.pickle and b/content/GA_1_4/auxiliary_files/m2_blue.pickle differ
diff --git a/content/GA_1_4/functions.py b/content/GA_1_4/functions.py
index 01c7f6a0bd1180488d6c6bf8602d5e86280a1518..a03c78100b4734ec5b17e2c9bcfab435f1737821 100644
--- a/content/GA_1_4/functions.py
+++ b/content/GA_1_4/functions.py
@@ -110,7 +110,7 @@ def plot_residual_histogram(d):
return fig, ax
-def update_plot(x0, x1, x2, x3, m):
+def model_widget(x0, x1, x2, x3, m):
plt.figure(figsize=(15,5))
plt.plot(m['times'], m['y'], 'o', label=m['data_type'])
plt.ylabel('Displacement [mm]')
@@ -119,20 +119,30 @@ def update_plot(x0, x1, x2, x3, m):
if x3 is None:
y_fit = m['A'] @ [x0, x1, x2]
if (x0 == 0) & (x1 == 0) & (x2 == 1):
- plt.plot(m['times'], y_fit, 'r', label='Groundwater data', linewidth=2)
+ plt.plot(m['times'], y_fit, 'r', label='Groundwater Interpolation', linewidth=2)
else:
plt.plot(m['times'], y_fit, 'r', label='Fit', linewidth=2)
else:
- y_fit = m['A'] @ [x0, x1, x2, x3]
- if (x0 == 0) & (x1 == 0) & (x2 == 0) & (x3 == 1):
- plt.plot(m['times'], y_fit, 'r', label='Groundwater data', linewidth=2)
+ functional_model = m['functional_model']
+ x_hat = (x0, x1, x2, x3)
+ y_fit = functional_model(x_hat, m)
+ if (x0 == 0) & (x1 == 0) & (x2 == 10) & (x3 == 1):
+ plt.plot(m['times'], y_fit, 'r', label='Groundwater Interpolation', linewidth=2)
else:
plt.plot(m['times'], y_fit, 'r', label='Fit', linewidth=2)
-
- W = np.linalg.inv(m['Sigma_Y'])
- ss_res = (m['y'] - y_fit).T @ W @ (m['y'] - y_fit)
- plt.title(f'Mean of squared residuals: {ss_res:.0f}')
+ plt.plot(m['groundwater_data']['times'],
+ m['groundwater_data']['y'],
+ 'ro', label='Groundwater Data',
+ markeredgecolor='black', markeredgewidth=1)
+
+ if 'Sigma_Y' in m:
+ W = np.linalg.inv(m['Sigma_Y'])
+ ss_res = (m['y'] - y_fit).T @ W @ (m['y'] - y_fit)
+ plt.title(f'Mean of squared residuals: {ss_res:.0f}')
+ else:
+ plt.title('Sigma_Y not yet defined in dictionary')
+
plt.ylim(-150, 30)
plt.grid()
plt.legend()
diff --git a/content/GA_1_4/setup.py b/content/GA_1_4/setup.py
index 9e2ccab293ae5134c47bb2e06d2c84933e1d7731..ae5fa19ffc89fd04324bb0fb952c82de69b6f8c1 100644
--- a/content/GA_1_4/setup.py
+++ b/content/GA_1_4/setup.py
@@ -36,14 +36,19 @@ m2['times'] = pd.to_datetime(gnss['times'])
m2['y'] = (gnss['observations[m]']).to_numpy()*1000
m2['days'], years_gnss = to_days_years(m2['times'])
-times_gw = pd.to_datetime(gw['times'])
-y_gw = (gw['observations[mm]']).to_numpy()
-days_gw, years_gw = to_days_years(times_gw)
+groundwater_data = {}
+groundwater_data['times'] = pd.to_datetime(gw['times'])
+groundwater_data['y'] = (gw['observations[mm]']).to_numpy()
+groundwater_data['days'], _ = to_days_years(groundwater_data['times'])
-interp = interpolate.interp1d(days_gw, y_gw)
+
+interp = interpolate.interp1d(groundwater_data['days'],
+ groundwater_data['y'])
m1['groundwater'] = interp(m1['days'])
m2['groundwater'] = interp(m2['days'])
+m1['groundwater_data'] = groundwater_data
+m2['groundwater_data'] = groundwater_data
# INSAR