diff --git a/synced_files/GA_1_7/Solution/Discharge/Analysis_discharge_solution.ipynb b/synced_files/GA_1_7/Solution/Discharge/Analysis_discharge_solution.ipynb
index e3f5d7f696379cc3e121abca5aafa9072adbcdf5..f45a35da3a9436db1850e484781001c3d5b4aa70 100644
--- a/synced_files/GA_1_7/Solution/Discharge/Analysis_discharge_solution.ipynb
+++ b/synced_files/GA_1_7/Solution/Discharge/Analysis_discharge_solution.ipynb
@@ -606,7 +606,7 @@
},
{
"cell_type": "markdown",
- "id": "b7129e8e",
+ "id": "2aba5273",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
diff --git a/synced_files/GA_1_7/Solution/Emissions/Analysis_emissions_solution.ipynb b/synced_files/GA_1_7/Solution/Emissions/Analysis_emissions_solution.ipynb
index 0533fbd802a6f2a4cdc6fe0ecda74c89bf682111..9b48d70d47b2384953fc01e38f5ea0c28d7e1a7a 100644
--- a/synced_files/GA_1_7/Solution/Emissions/Analysis_emissions_solution.ipynb
+++ b/synced_files/GA_1_7/Solution/Emissions/Analysis_emissions_solution.ipynb
@@ -159,7 +159,7 @@
},
{
"cell_type": "markdown",
- "id": "8cfea392",
+ "id": "86061aa8",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
@@ -262,7 +262,7 @@
},
{
"cell_type": "markdown",
- "id": "417d1089",
+ "id": "1bb7c655",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
@@ -374,7 +374,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "8b0f7732",
+ "id": "a93719e4",
"metadata": {},
"outputs": [],
"source": [
@@ -437,7 +437,7 @@
},
{
"cell_type": "markdown",
- "id": "2774ec22",
+ "id": "424cf35a",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
@@ -538,7 +538,7 @@
},
{
"cell_type": "markdown",
- "id": "923f6481",
+ "id": "4cc75538",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
@@ -552,7 +552,7 @@
},
{
"cell_type": "markdown",
- "id": "be078334",
+ "id": "d7552937",
"metadata": {},
"source": [
"If you run the code in the cell below, you will obtain a scatter plot of both variables. Explore the relationship between both variables and answer the following questions:\n",
@@ -573,7 +573,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "15c62261",
+ "id": "00efa69d",
"metadata": {},
"outputs": [],
"source": [
@@ -589,7 +589,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "f506d847",
+ "id": "8c9e2eeb",
"metadata": {},
"outputs": [],
"source": [
@@ -602,7 +602,7 @@
},
{
"cell_type": "markdown",
- "id": "d01c9497",
+ "id": "c57e8e7c",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
diff --git a/synced_files/GA_1_7/Solution/Force/Analysis_force_solution.ipynb b/synced_files/GA_1_7/Solution/Force/Analysis_force_solution.ipynb
index 0fb50c38996a524b924408f2e5e5060348c6a44c..ea2403e27960fed8bdfc90e31ae645c996c34da3 100644
--- a/synced_files/GA_1_7/Solution/Force/Analysis_force_solution.ipynb
+++ b/synced_files/GA_1_7/Solution/Force/Analysis_force_solution.ipynb
@@ -164,7 +164,7 @@
},
{
"cell_type": "markdown",
- "id": "ab28e757",
+ "id": "096ccdb1",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
@@ -269,7 +269,7 @@
},
{
"cell_type": "markdown",
- "id": "0ebf7fed",
+ "id": "048967da",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
@@ -382,7 +382,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "ebb8fa61",
+ "id": "49ec3949",
"metadata": {},
"outputs": [],
"source": [
@@ -445,7 +445,7 @@
},
{
"cell_type": "markdown",
- "id": "1d352601",
+ "id": "500948cb",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
@@ -546,7 +546,7 @@
},
{
"cell_type": "markdown",
- "id": "64197dc5",
+ "id": "e100f74b",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
@@ -561,7 +561,7 @@
},
{
"cell_type": "markdown",
- "id": "ab892dfa",
+ "id": "bb0c2545",
"metadata": {},
"source": [
"If you run the code in the cell below, you will obtain a scatter plot of both variables. Explore the relationship between both variables and answer the following questions:\n",
@@ -582,7 +582,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "78a28820",
+ "id": "13722902",
"metadata": {},
"outputs": [],
"source": [
@@ -598,7 +598,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "aa746d9d",
+ "id": "0e8d2a14",
"metadata": {},
"outputs": [],
"source": [
@@ -611,7 +611,7 @@
},
{
"cell_type": "markdown",
- "id": "43263510",
+ "id": "db232608",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
diff --git a/synced_files/GA_1_7/Student/Discharge/Analysis_discharge.ipynb b/synced_files/GA_1_7/Student/Discharge/Analysis_discharge.ipynb
index ebd2eacb012af4e54b21dcce234980712c9aee34..31aa299cc07f34c515b79dbcfd599ca9f2b6d7ea 100644
--- a/synced_files/GA_1_7/Student/Discharge/Analysis_discharge.ipynb
+++ b/synced_files/GA_1_7/Student/Discharge/Analysis_discharge.ipynb
@@ -384,7 +384,7 @@
},
{
"cell_type": "markdown",
- "id": "a77e757b",
+ "id": "09d3bdef",
"metadata": {},
"source": [
"If you run the code in the cell below, you will obtain a scatter plot of both variables. Explore the relationship between both variables and answer the following questions:\n",
@@ -405,7 +405,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "d92acf77",
+ "id": "5ffa38b1",
"metadata": {},
"outputs": [],
"source": [
@@ -421,7 +421,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "0b7f02b3",
+ "id": "b6ac7c8a",
"metadata": {},
"outputs": [],
"source": [
diff --git a/synced_files/GA_1_7/Student/Emissions/Analysis_emissions.ipynb b/synced_files/GA_1_7/Student/Emissions/Analysis_emissions.ipynb
index db52cd423c6709a14eaab7c0b4b60e24348b5ece..332f3f700d5b74c117de84620c38be1ebe960ebf 100644
--- a/synced_files/GA_1_7/Student/Emissions/Analysis_emissions.ipynb
+++ b/synced_files/GA_1_7/Student/Emissions/Analysis_emissions.ipynb
@@ -378,7 +378,7 @@
},
{
"cell_type": "markdown",
- "id": "9c309df1",
+ "id": "efdf2897",
"metadata": {},
"source": [
"If you run the code in the cell below, you will obtain a scatter plot of both variables. Explore the relationship between both variables and answer the following questions:\n",
@@ -399,7 +399,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "a7c914c5",
+ "id": "3ea28524",
"metadata": {},
"outputs": [],
"source": [
@@ -415,7 +415,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "d0d94667",
+ "id": "812ff8af",
"metadata": {},
"outputs": [],
"source": [
diff --git a/synced_files/GA_1_7/Student/Force/Analysis_force.ipynb b/synced_files/GA_1_7/Student/Force/Analysis_force.ipynb
index 4756a662afb568c24098c1bd828a7bdf75388dc6..c3474d97eb4c149c355228afd49b11d8614ed01d 100644
--- a/synced_files/GA_1_7/Student/Force/Analysis_force.ipynb
+++ b/synced_files/GA_1_7/Student/Force/Analysis_force.ipynb
@@ -384,7 +384,7 @@
},
{
"cell_type": "markdown",
- "id": "faaad0f2",
+ "id": "14f0731a",
"metadata": {},
"source": [
"If you run the code in the cell below, you will obtain a scatter plot of both variables. Explore the relationship between both variables and answer the following questions:\n",
@@ -405,7 +405,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "24488ae4",
+ "id": "2547841f",
"metadata": {},
"outputs": [],
"source": [
@@ -421,7 +421,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "a2d21fe8",
+ "id": "211edbaf",
"metadata": {},
"outputs": [],
"source": [
diff --git a/synced_files/GA_1_7/Unused/Temp/Distribution_Fitting_T.ipynb b/synced_files/GA_1_7/Unused/Temp/Distribution_Fitting_T.ipynb
index 3eb5bc832c22971b8ba9ab5942a4421b3b90651b..21d918cbb31cb4a3d016a2eee579563b3c933693 100644
--- a/synced_files/GA_1_7/Unused/Temp/Distribution_Fitting_T.ipynb
+++ b/synced_files/GA_1_7/Unused/Temp/Distribution_Fitting_T.ipynb
@@ -163,7 +163,7 @@
},
{
"cell_type": "markdown",
- "id": "3b11a8c8",
+ "id": "005b7091",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
@@ -263,7 +263,7 @@
},
{
"cell_type": "markdown",
- "id": "0deb9cca",
+ "id": "6337c5ab",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
@@ -372,7 +372,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "1a8fe54b",
+ "id": "10b1c5a2",
"metadata": {},
"outputs": [],
"source": [
@@ -434,7 +434,7 @@
},
{
"cell_type": "markdown",
- "id": "bbfb882b",
+ "id": "454cbbae",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
@@ -533,7 +533,7 @@
},
{
"cell_type": "markdown",
- "id": "72a21985",
+ "id": "97f0b74d",
"metadata": {},
"source": [
"<div style=\"background-color:#FAE99E; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
@@ -546,7 +546,7 @@
},
{
"cell_type": "markdown",
- "id": "af8f794d",
+ "id": "847abf2d",
"metadata": {},
"source": [
"If you run the code in the cell below, you will obtain a scatter plot of both variables. Explore the relationship between both variables and answer the following questions:\n",
@@ -567,7 +567,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "010bb291",
+ "id": "5cc3cedc",
"metadata": {},
"outputs": [],
"source": [
diff --git a/synced_files/Week_2_2/old/PA_2_1_solution_sympy.ipynb b/synced_files/Week_2_2/old/PA_2_1_solution_sympy.ipynb
index 7515e05a990ef6fce3b0dc26ac09cd79139410ae..149cd64f2202d05f4ecffc747e4c512ad30599bd 100644
--- a/synced_files/Week_2_2/old/PA_2_1_solution_sympy.ipynb
+++ b/synced_files/Week_2_2/old/PA_2_1_solution_sympy.ipynb
@@ -404,7 +404,7 @@
},
{
"cell_type": "markdown",
- "id": "624efe8f",
+ "id": "373e8503",
"metadata": {},
"source": [
"You could also solve this problem using sympy! What would be the benefit of doing this? Check below how long it will take! You won't need timeit for this one..."
@@ -413,7 +413,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "93726bb0",
+ "id": "0e6b5a13",
"metadata": {},
"outputs": [],
"source": [
@@ -423,7 +423,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "e9ac4516",
+ "id": "8308b7f9",
"metadata": {},
"outputs": [],
"source": [
@@ -433,7 +433,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "e70f950a",
+ "id": "6bfc7a47",
"metadata": {},
"outputs": [],
"source": [
@@ -445,7 +445,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "329c37c6",
+ "id": "1b562e5b",
"metadata": {},
"outputs": [],
"source": [
@@ -456,7 +456,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "0e0516e0",
+ "id": "0db38f94",
"metadata": {},
"outputs": [],
"source": [
@@ -465,7 +465,7 @@
},
{
"cell_type": "markdown",
- "id": "59a4084e",
+ "id": "f250b6ef",
"metadata": {},
"source": [
"What is the result for the 501th value using float values?"
@@ -474,7 +474,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "3f1418f5",
+ "id": "03e6e369",
"metadata": {},
"outputs": [],
"source": [
@@ -514,7 +514,7 @@
},
{
"cell_type": "markdown",
- "id": "91bae574",
+ "id": "aa2d4fd3",
"metadata": {},
"source": []
}
diff --git a/synced_files/Week_2_2/old/old_PA10_solution_sympy.ipynb b/synced_files/Week_2_2/old/old_PA10_solution_sympy.ipynb
index 5825760df5d0b24341eb3c05c4ff14943ad17c65..bfd1003805f0435cc0d202143475e32413196d21 100644
--- a/synced_files/Week_2_2/old/old_PA10_solution_sympy.ipynb
+++ b/synced_files/Week_2_2/old/old_PA10_solution_sympy.ipynb
@@ -404,7 +404,7 @@
},
{
"cell_type": "markdown",
- "id": "7b165f69",
+ "id": "d7b98994",
"metadata": {},
"source": [
"You could also solve this problem using sympy! What would be the benefit of doing this? Check below how long it will take! You won't need timeit for this one..."
@@ -413,7 +413,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "2a506f9f",
+ "id": "7631b1aa",
"metadata": {},
"outputs": [],
"source": [
@@ -423,7 +423,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "43ba02f7",
+ "id": "74bfd6f5",
"metadata": {},
"outputs": [],
"source": [
@@ -433,7 +433,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "f477d554",
+ "id": "b1283a67",
"metadata": {},
"outputs": [],
"source": [
@@ -445,7 +445,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "cb890c60",
+ "id": "f33e4806",
"metadata": {},
"outputs": [],
"source": [
@@ -456,7 +456,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "ef9e5aed",
+ "id": "8389b3dc",
"metadata": {},
"outputs": [],
"source": [
@@ -465,7 +465,7 @@
},
{
"cell_type": "markdown",
- "id": "3a22ca72",
+ "id": "4126904f",
"metadata": {},
"source": [
"What is the result for the 501th value using float values?"
@@ -474,7 +474,7 @@
{
"cell_type": "code",
"execution_count": null,
- "id": "ff902516",
+ "id": "3eb834dd",
"metadata": {},
"outputs": [],
"source": [
@@ -514,7 +514,7 @@
},
{
"cell_type": "markdown",
- "id": "7d2b74af",
+ "id": "d6b16644",
"metadata": {},
"source": []
}
diff --git a/synced_files/test.ipynb b/synced_files/test.ipynb
deleted file mode 100644
index 586d353ad6ab10ccf2fcf1b546c55b9e27232402..0000000000000000000000000000000000000000
--- a/synced_files/test.ipynb
+++ /dev/null
@@ -1,1483 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "id": "785b4e1d",
- "metadata": {},
- "source": [
- "# GA 1.3: test\n",
- "\n",
- "<h1 style=\"position: absolute; display: flex; flex-grow: 0; flex-shrink: 0; flex-direction: row-reverse; top: 60px;right: 30px; margin: 0; border: 0\">\n",
- " <style>\n",
- " .markdown {width:100%; position: relative}\n",
- " article { position: relative }\n",
- " </style>\n",
- " <img src=\"https://gitlab.tudelft.nl/mude/public/-/raw/main/tu-logo/TU_P1_full-color.png\" style=\"width:100px\" />\n",
- " <img src=\"https://gitlab.tudelft.nl/mude/public/-/raw/main/mude-logo/MUDE_Logo-small.png\" style=\"width:100px\" />\n",
- "</h1>\n",
- "<h2 style=\"height: 10px\">\n",
- "</h2>\n",
- "\n",
- "*[CEGM1000 MUDE](http://mude.citg.tudelft.nl/): Week 1.3. Due: Friday, September 20, 2024.*"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "4ad8b9cb",
- "metadata": {},
- "source": [
- "<div style=\"background-color:#ffa6a6; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%\"><p><b>Note:</b> don't forget to read the \"Assignment Context\" section of the README, it contains important information to understand this analysis.</p></div>"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "181ccfd5",
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "import numpy as np\n",
- "from scipy import interpolate\n",
- "from scipy.stats import norm\n",
- "import pandas as pd\n",
- "import matplotlib.pyplot as plt\n",
- "\n",
- "\n",
- "from functions import *\n",
- "\n",
- "np.set_printoptions(precision=3)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "5ca94e0e",
- "metadata": {},
- "source": [
- "## Part 0: Dictionary Review\n",
- "\n",
- "As described above, several functions in this assignment require the use of a Python dictionary to make it easier to keep track of important data, variables and results for the various _models_ we will be constructing and validating.\n",
- "\n",
- "_It may be useful to revisit PA 1.1, where there was a brief infroduction to dictionaires. That PA contains all the dictionary info you need for GA 1.3. A [read-only copy is here](https://mude.citg.tudelft.nl/2024/files/Week_1_1/PA_1_1_Catch_Them_All.html) and [the source code (notebook) is here](https://gitlab.tudelft.nl/mude/2024-week-1-1)._"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d8c39791",
- "metadata": {},
- "source": [
- "<div style=\"background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
- "<p>\n",
- "\n",
- "$\\textbf{Task 0.1}$ \n",
- " \n",
- "Read and run the cell below to make sure you remember how to use a dictionary.\n",
- "\n",
- "Modify the function to print some of the other key-value pairs of the dictionary.\n",
- "\n",
- "<em>It may also be useful to use the cell below when working on later tasks in this assignment.</em>\n",
- " \n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "8683b5ed",
- "metadata": {},
- "outputs": [],
- "source": [
- "my_dictionary = {'key1': 'value1',\n",
- " 'key2': 'value2',\n",
- " 'name': 'Dictionary Example',\n",
- " 'a_list': [1, 2, 3],\n",
- " 'an_array': np.array([1, 2, 3]),\n",
- " 'a_string': 'hello'\n",
- " }\n",
- "\n",
- "def function_that_uses_my_dictionary(d):\n",
- " print(d['key1'])\n",
- "\n",
- " # SOLUTION:\n",
- " print(d['name'])\n",
- " print(d['a_list'])\n",
- " print(d['an_array'])\n",
- " print(d['a_string'])\n",
- "\n",
- " if 'new_key' in d:\n",
- " print('new_key exists and has value:', d['new_key'])\n",
- " return\n",
- "\n",
- "function_that_uses_my_dictionary(my_dictionary)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "86bc7f97",
- "metadata": {},
- "source": [
- "<div style=\"background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
- "<p>\n",
- "\n",
- "$\\textbf{Task 0.2}$ \n",
- "\n",
- "Test your knowledge by adding a new key <code>new_key</code> and then executing the function to print the value.\n",
- " \n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "41c56f43",
- "metadata": {},
- "outputs": [],
- "source": [
- "# YOUR_CODE_HERE\n",
- "# function_that_uses_my_dictionary(my_dictionary)\n",
- "\n",
- "# SOLUTION:\n",
- "my_dictionary['new_key'] = 'new_value'\n",
- "function_that_uses_my_dictionary(my_dictionary)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "160d6250",
- "metadata": {
- "id": "160d6250"
- },
- "source": [
- "## Task 1: Preparing the data\n",
- "\n",
- "Within this assignment you will work with two types of data: InSAR data and GNSS data. The cell below will load the data and visualize the observed displacements time. In this task we use the package `pandas`, which is really useful for handling time series. We will learn how to use it later in the quarter; for now, you only need to recognize that it imports the data as a `dataframe` object, which we then convert into a numpy array using the code below."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "02b12781",
- "metadata": {},
- "source": [
- "<div style=\"background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\"> <p>Tip: note that we have converted all observations to millimeters.</p></div>"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "1f28eba3",
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "gnss = pd.read_csv('./data/gnss_observations.csv')\n",
- "times_gnss = pd.to_datetime(gnss['times'])\n",
- "y_gnss = (gnss['observations[m]']).to_numpy()*1000\n",
- "\n",
- "insar = pd.read_csv('./data/insar_observations.csv')\n",
- "times_insar = pd.to_datetime(insar['times'])\n",
- "y_insar = (insar['observations[m]']).to_numpy()*1000\n",
- "\n",
- "gw = pd.read_csv('./data/groundwater_levels.csv')\n",
- "times_gw = pd.to_datetime(gw['times'])\n",
- "y_gw = (gw['observations[mm]']).to_numpy()"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "aa8906a9-2ebe-4432-b4f3-8bf074d6b181",
- "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",
- "Once you have used the cell above to import the data, investigate the data sets using the code cell below. Then provide some relevant summary information in the Markdown cell.\n",
- "\n",
- "<em>Hint: at the least, you should be able to tell how many data points are in each data set and get an understanding of the mean and standard deviation of each. Make sure you compare the different datasets and use consistent units.</em>\n",
- " \n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "12321724-dc58-42ee-9bd5-44265d3bc921",
- "metadata": {
- "id": "0491cc69"
- },
- "source": [
- "<div style=\"background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\"> <p>The code below gives some examples of the quantitative and qualitative ways you could have looked at the data. It is more than you were expected to do; the important thing is that you showed the ability to learn something about the data and describe aspects that are relevant to our problem. We use a dictionary to easily access the different data series using their names, which are entered as the dictionary keys (also not expected of you, but it's hopefully fun to learn useful tricks).</div>"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "9f025cfc-4f89-492d-ac26-f5b6381d0c70",
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "# YOUR_CODE_HERE\n",
- "\n",
- "# SOLUTION:\n",
- "data_list = ['y_gnss', 'y_insar', 'y_gw']\n",
- "data_dict = {data_list[0]: y_gnss,\n",
- " data_list[1]: y_insar,\n",
- " data_list[2]: y_gw}\n",
- "def print_summary(data):\n",
- " '''Summarize an array with simple print statements.'''\n",
- " print('Minimum = ', data.min())\n",
- " print('Maximum = ', data.max())\n",
- " print('Mean = ', data.mean())\n",
- " print('Std dev = ', data.std())\n",
- " print('Shape = ', data.shape)\n",
- " print('First value = ', data[0])\n",
- " print('Last value = ', data[-1])\n",
- " print('\\n')\n",
- " \n",
- "for item in data_list:\n",
- " print('Summary for array: ', item)\n",
- " print('------------------------------------------------')\n",
- " print_summary(data_dict[item])"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "71cf2133-37a6-4536-82a6-42a46b8a1c66",
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "# SOLUTION:\n",
- "times_dict = {data_list[0]: times_gnss,\n",
- " data_list[1]: times_insar,\n",
- " data_list[2]: times_gw}\n",
- "def plot_data(times, data, label):\n",
- " plt.figure(figsize=(15,4))\n",
- " plt.plot(times, data, 'co', mec='black')\n",
- " plt.title(label)\n",
- " plt.xlabel('Times')\n",
- " plt.ylabel('Data [mm]')\n",
- " plt.show()\n",
- "\n",
- "plt.figure(figsize=(15,4))\n",
- "for i in range(3):\n",
- " plot_data(times_dict[data_list[i]],\n",
- " data_dict[data_list[i]],\n",
- " data_list[i])"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "a9c02e8f-81f9-41a3-b894-c23dd9617207",
- "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",
- "There are a lot more GNSS data points than InSAR or groundwater. The GNSS observations also have more noise, and what seem to be outliers. In this case the mean and standard deviation do not mean much, because there is clearly a trend with time. We can at least confirm that the time periods of measurements overlap, although the intervals between measurements is certainly not uniform (note that you don't need to do anything with the times, since they are pandas time series and we have not covered them yet).\n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "9fe5a729",
- "metadata": {},
- "source": [
- "You may have noticed that the groundwater data is available for different times than the GNSS and InSAR data. You will therefore have to *interpolate* the data to the same times for a further analysis. You can use the SciPy function ```interpolate.interp1d``` (read its [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html)).\n",
- "\n",
- "The cells below do the following:\n",
- "1. Define a function to convert the time unit\n",
- "2. Convert the time stamps for all data\n",
- "3. Use `interp1d` to interpolate the groundwater measurements at the time of the satellite measurements"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "f02ed4c4",
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "def to_days_years(times):\n",
- " '''Convert the observation times to days and years.'''\n",
- " \n",
- " times_datetime = pd.to_datetime(times)\n",
- " time_diff = (times_datetime - times_datetime[0])\n",
- " days_diff = (time_diff / np.timedelta64(1,'D')).astype(int)\n",
- " \n",
- " days = days_diff.to_numpy()\n",
- " years = days/365\n",
- " \n",
- " return days, years"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "edf14892",
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "days_gnss, years_gnss = to_days_years(times_gnss)\n",
- "days_insar, years_insar = to_days_years(times_insar)\n",
- "days_gw, years_gw = to_days_years(times_gw)\n",
- "\n",
- "interp = interpolate.interp1d(days_gw, y_gw)\n",
- "\n",
- "GW_at_GNSS_times = interp(days_gnss)\n",
- "GW_at_InSAR_times = interp(days_insar)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "16827704-4fc5-4afb-879e-0ceeac45eb18",
- "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",
- "Answer/complete the code and Markdown cells below:\n",
- "<ol>\n",
- " <li>What is <code>interp</code>? (what kind of object is it, and how does it work?)</li>\n",
- " <li>How did the groundwater observation array change? Be quantitative. </li>\n",
- "</ol>\n",
- " \n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "6cdfb46b-1324-4c2b-8148-5a6a102ede2e",
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "# YOUR_CODE_HERE\n",
- "\n",
- "# SOLUTION:\n",
- "print('array size of GW_at_GNSS_times', len(GW_at_GNSS_times))\n",
- "print('array size of GW_at_InSAR_times', len(GW_at_InSAR_times))\n",
- "print('array size of GW before interpolation', len(y_gw))\n",
- "\n",
- "print('\\nFirst values of times_gw:')\n",
- "print(times_gw[0:2])\n",
- "print('\\nFirst values of y_gw:')\n",
- "print(y_gw[0:2])\n",
- "print('\\nFirst values of times_gnss:')\n",
- "print(times_gnss[0:2])\n",
- "print('\\nFirst values of GW_at_GNSS_times:')\n",
- "print(GW_at_GNSS_times[0:2])"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "3b8c68eb-9774-4c3c-91da-f29f035b178c",
- "metadata": {},
- "source": [
- "**Write your answer in this Markdown cell.**"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "dcc0e2a3-3f34-4bde-af54-37923cb803cd",
- "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",
- "<ol>\n",
- " <li><code>interp</code> is a function that will return a value (gw level) for the input(s) (date(s)). The interpolated value is found by linearly interpolating between the two nearest times in the gw observations.</li>\n",
- " <li>The observation arrays of <code>GW_at_GNSS_times</code> and <code>GW_at_INSAR_times</code> changed in size to match the size of the GNSS and InSAR observations, respectively.</li>\n",
- "</ol> \n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "349ebd38-d9e1-49b4-b6bd-f5dc399727e0",
- "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",
- "Create a single plot to compare observed displacement for the GNSS and InSAR data sets.\n",
- " \n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "e868e488",
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "# plt.figure(figsize=(15,5))\n",
- "# plt.plot(YOUR_CODE_HERE, YOUR_CODE_HERE,\n",
- "# 'o', mec='black', label = 'GNSS')\n",
- "# plt.plot(YOUR_CODE_HERE, YOUR_CODE_HERE,\n",
- "# 'o', mec='black', label = 'InSAR')\n",
- "# plt.legend()\n",
- "# plt.ylabel('Displacement [mm]')\n",
- "# plt.xlabel('Time')\n",
- "# plt.show()\n",
- "\n",
- "# SOLUTION:\n",
- "plt.figure(figsize=(15,5))\n",
- "plt.plot(times_gnss, y_gnss, 'o', mec='black', label = 'GNSS')\n",
- "plt.plot(times_insar, y_insar, 'o', mec='black', label = 'InSAR')\n",
- "plt.legend()\n",
- "plt.ylabel('Displacement [mm]')\n",
- "plt.xlabel('Time')\n",
- "plt.show()"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "c9b45b8e",
- "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.4:</b> \n",
- "Describe the datasets based on the figure above and your observations from the previous tasks. What kind of deformation do you see? And what are the differences between both datasets? Be quantitative.\n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "48bae179",
- "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 points obviously show subsidence, the displacement shows a similar pattern for both datasets. The GNSS data is much noisier than InSAR (range is around 60 mm versus only a few mm), but has a higher sampling rate. Also there seem to be more outliers in the GNSS data compared to InSAR, especially at the start of the observation period. InSAR has only observations every 6 days but is less noisy. \n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f9a7bdd4",
- "metadata": {},
- "source": [
- "Before we move on, it is time to do a little bit of housekeeping.\n",
- "\n",
- "Have you found it confusing to keep track of two sets of variables---one for each data type? Let's use a dictionary to store relevant information about each model. We will use this in the plotting functions for this task (and again next week), so make sure you take the time to see what is happening. Review also Part 0 at the top of this notebook if you need a refresher on dictionaries."
- ]
- },
- {
- "cell_type": "markdown",
- "id": "cdcf20f7",
- "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.5:</b> \n",
- " \n",
- "Run the cell below to define a dictionary for storing information about the two (future) models.\n",
- "\n",
- "</p>\n",
- "</div>\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "2c27b4f3",
- "metadata": {},
- "outputs": [],
- "source": [
- "model_insar = {'data_type': 'InSAR',\n",
- " 'y':y_insar,\n",
- " 'times':times_insar,\n",
- " 'groundwater': GW_at_InSAR_times\n",
- " }\n",
- "\n",
- "model_gnss = {'data_type': 'GNSS',\n",
- " 'y':y_gnss,\n",
- " 'times':times_gnss,\n",
- " 'groundwater': GW_at_GNSS_times\n",
- " }"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "76c9115b",
- "metadata": {
- "id": "76c9115b"
- },
- "source": [
- "## Task 2: Set-up linear functional model\n",
- "\n",
- "We want to investigate how we could model the observed displacements of the road. Because the road is built in the Green Heart we expect that the observed displacements are related to the groundwater level. Furthermore, we assume that the displacements can be modeled using a constant velocity. The model is defined as \n",
- "$$\n",
- "d = d_0 + vt + k \\ \\textrm{GW},\n",
- "$$\n",
- "where $d$ is the displacement, $t$ is time and $\\textrm{GW}$ is the groundwater level (that we assume to be deterministic). \n",
- "\n",
- "Therefore, the model has 3 unknowns:\n",
- "1. $d_0$, as the initial displacement at $t_0$;\n",
- "2. $v$, as the displacement velocity;\n",
- "3. $k$, as the 'groundwater factor', which can be seen as the response of the soil to changes in the groundwater level.\n",
- "\n",
- "\n",
- "As a group you will construct the **functional model** that is defined as \n",
- "$$\n",
- "\\mathbb{E}(Y) = \\mathrm{A x}.\n",
- "$$\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f6aca691",
- "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 2.1:</b> \n",
- " \n",
- "Construct the design matrix $A$ (for both InSAR and GNSS observations), then show the first 5 observations and confirm the dimensions of $A$.\n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "3a3eb1a1",
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "# YOUR_CODE_HERE\n",
- "\n",
- "# SOLUTION:\n",
- "A_insar = np.ones((len(times_insar), 3))\n",
- "A_insar[:,1] = days_insar\n",
- "A_insar[:,2] = GW_at_InSAR_times\n",
- "\n",
- "print ('The first 5 rows of the A matrix (InSAR) are:')\n",
- "print (A_insar[0:5, :])\n",
- "\n",
- "print ('The first 5 observations [mm] of y_insar are:')\n",
- "print (y_insar[0:5])\n",
- "\n",
- "m_insar = np.shape(A_insar)[0]\n",
- "n_insar = np.shape(A_insar)[1]\n",
- "print(f'm = {m_insar} and n = {n_insar}')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "4bcd395d",
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "# YOUR_CODE_HERE\n",
- "\n",
- "# SOLUTION:\n",
- "A_gnss = np.ones((len(times_gnss), 3))\n",
- "A_gnss[:,1] = days_gnss\n",
- "A_gnss[:,2] = GW_at_GNSS_times\n",
- "\n",
- "print ('The first 5 rows of the A matrix (GNSS) are:')\n",
- "print (A_gnss[0:5, :])\n",
- "\n",
- "print ('\\nThe first 5 observations [mm] of y_gnss are:')\n",
- "print (y_gnss[0:5])\n",
- "\n",
- "m_gnss = np.shape(A_gnss)[0]\n",
- "n_gnss = np.shape(A_gnss)[1]\n",
- "print(f'm = {m_gnss} and n = {n_gnss}')"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d390f466",
- "metadata": {},
- "source": [
- "<div style=\"background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
- "<p>\n",
- "\n",
- "$\\textbf{Task 2.2}$\n",
- " \n",
- "Answer the following questions:\n",
- "\n",
- "- What is the dimension of the observables' vector $Y$?\n",
- "- What are the unknowns of the functional model?\n",
- "- What is the redundancy for this model?\n",
- "\n",
- "</p>\n",
- "</div>\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "d40a0ecf",
- "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",
- "For InSAR:\n",
- "<ol>\n",
- " <li>The number of observations is 61.</li>\n",
- " <li>The number of unknowns is 3.</li>\n",
- " <li>The redundancy is 58.</li>\n",
- "</ol>\n",
- " \n",
- "For GNSS:\n",
- "<ol>\n",
- " <li>The number of observations is 730.</li>\n",
- " <li>The number of unknowns is 3.</li>\n",
- " <li>The redundancy is 727.</li>\n",
- "</ol> \n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "cde2f4db",
- "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 2.3:</b> \n",
- " \n",
- "Add the A matrix to the dictionaries for each model. This will be used to plot results later in the notebook.\n",
- "\n",
- "</p>\n",
- "</div>\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "396ac3a5",
- "metadata": {},
- "outputs": [],
- "source": [
- "# model_insar['A'] = YOUR_CODE_HERE\n",
- "# model_gnss['A'] = YOUR_CODE_HERE\n",
- "\n",
- "# SOLUTION:\n",
- "model_insar['A'] = A_insar\n",
- "model_gnss['A'] = A_gnss\n",
- "\n",
- "print(\"Keys and Values (type) for model_insar:\")\n",
- "for key, value in model_insar.items():\n",
- " print(f\"{key:16s} --> {type(value)}\")\n",
- "print(\"\\nKeys and Values (type) for model_gnss:\")\n",
- "for key, value in model_gnss.items():\n",
- " print(f\"{key:16s} --> {type(value)}\")"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "9325d32b",
- "metadata": {
- "id": "9325d32b",
- "tags": []
- },
- "source": [
- "## 3. Set-up stochastic model\n",
- "\n",
- "We will use the Best Linear Unbiased Estimator (BLUE) to solve for the unknown parameters. Therefore we also need a stochastic model, which is defined as\n",
- "$$\n",
- "\\mathbb{D}(Y) = \\Sigma_{Y}.\n",
- "$$\n",
- "where $\\Sigma_{Y}$ is the covariance matrix of the observables' vector. \n",
- "\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "dc3aec4c",
- "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 3.1:</b> \n",
- " \n",
- "Construct the covariance matrix for each type of data and assume that \n",
- "\n",
- "- the observables are independent\n",
- "\n",
- "- the observables are normally distributed\n",
- "\n",
- "- the observables' standard deviation is\n",
- " \n",
- " - $\\sigma_\\textrm{InSAR} = 2$ mm \n",
- " - $\\sigma_\\textrm{GNSS} = 15$ mm \n",
- " \n",
- "</p>\n",
- "</div>\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "163acdb3",
- "metadata": {
- "colab": {
- "base_uri": "https://localhost:8080/"
- },
- "executionInfo": {
- "elapsed": 40,
- "status": "ok",
- "timestamp": 1664699881875,
- "user": {
- "displayName": "C Yin",
- "userId": "14075875094781565898"
- },
- "user_tz": -120
- },
- "id": "163acdb3",
- "outputId": "8bc99da3-a61e-4a25-8a54-c90f33299a57",
- "tags": []
- },
- "outputs": [],
- "source": [
- "# YOUR_CODE_HERE\n",
- "\n",
- "# SOLUTION:\n",
- "std_insar = 2 #mm\n",
- "\n",
- "Sigma_Y_insar = np.identity(len(times_insar))*std_insar**2\n",
- "\n",
- "print ('Sigma_Y (InSAR) is defined as:')\n",
- "print (Sigma_Y_insar)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "5d583bd8",
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "# YOUR_CODE_HERE\n",
- "\n",
- "# SOLUTION:\n",
- "std_gnss = 15 #mm (corrected from original value of 5 mm)\n",
- "\n",
- "Sigma_Y_gnss = np.identity(len(times_gnss))*std_gnss**2\n",
- "\n",
- "print ('\\nSigma_Y (GNSS) is defined as:')\n",
- "print (Sigma_Y_gnss)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "b2665071",
- "metadata": {},
- "source": [
- "<div style=\"background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\">\n",
- "<p>\n",
- "\n",
- "$\\textbf{Task 3.2}$\n",
- " \n",
- "Answer the following questions:\n",
- "\n",
- "- What information is contained in the covariance matrix?\n",
- "- How do you implement the assumption that all observations are independent?\n",
- "- What is the dimension of $\\Sigma_{Y}$?\n",
- "- How do you create $\\Sigma_{Y}$?\n",
- "\n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "4322a5b5",
- "metadata": {},
- "source": [
- "_Write your answer in this cell._"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "124ddd77",
- "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 covariance matrix contains information on the quality of the observations, where an entry on the diagonal represents the variance of one observation at a particular epoch. If there is an indication that for instance the quality for a particular time interval differs, different $\\sigma$ values can be put in the stochastic model for these epochs. \n",
- "- The off-diagonal terms in the matrix are related to the correlation between observations at different epochs, where a zero value on the off-diagonal indicates zero correlation.\n",
- "- The dimension of the matrix is 61x61 for InSAR and 730x730 for GNSS.\n",
- "- See code.\n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8e87d1fe",
- "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 3.3:</b> \n",
- " \n",
- "Add <code>Sigma_Y</code> to the dictionaries for each model.\n",
- "\n",
- "</p>\n",
- "</div>\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "13cffb4c",
- "metadata": {},
- "outputs": [],
- "source": [
- "# model_insar['Sigma_Y] = YOUR_CODE_HERE\n",
- "# model_gnss['Sigma_Y'] = YOUR_CODE_HERE\n",
- "\n",
- "# SOLUTION:\n",
- "model_insar['Sigma_Y'] = Sigma_Y_insar\n",
- "model_gnss['Sigma_Y'] = Sigma_Y_gnss"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "09e965bf",
- "metadata": {
- "id": "09e965bf"
- },
- "source": [
- "## 4. Apply best linear unbiased estimation\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "83ec4a48",
- "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 4.1:</b> \n",
- " \n",
- "Write a function to apply BLUE in the cell below and use the function to estimate the unknowns for the model using the data.\n",
- "\n",
- "Compute the modeled displacements ($\\hat{\\mathrm{y}}$), and corresponding residuals ($\\hat{\\mathrm{\\epsilon}}$), as well as associated values (as requested by the blank code lines).\n",
- "</p>\n",
- "</div>\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "936d6b0c",
- "metadata": {},
- "source": [
- "<div style=\"background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\"> <p><strong>Note on code implementation</strong>: you'll see that the functions in this assignment use a dictionary; this greatly reduces the number of input/output variables needed in a function. However, it can make the code inside the function more difficult to read due to the key syntax (e.g., <code>dict['variable_1']</code> versus <code>variable\n",
- "_1</code>). To make this assignment easier for you to implement we have split these functions into three parts: 1) define variables from the dictionary, 2) perform analysis, 3) add results to the dictionary. Note that this is not the most efficient way to write this code; it is done here specifically for clarity and to help you focus on writing the equations properly and understanding the meaning of each term.</p></div>"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "d85b1826",
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "def BLUE(d):\n",
- " \"\"\"Calculate the Best Linear Unbiased Estimator\n",
- " \n",
- " Uses dict as input/output:\n",
- " - inputs defined from existing values in dict\n",
- " - outputs defined as new values in dict\n",
- " \"\"\"\n",
- "\n",
- " y = d['y']\n",
- " A = d['A']\n",
- " Sigma_Y = d['Sigma_Y']\n",
- "\n",
- " # Sigma_X_hat = YOUR_CODE_HERE\n",
- " # x_hat = YOUR_CODE_HERE\n",
- " \n",
- " # y_hat = YOUR_CODE_HERE\n",
- "\n",
- " # e_hat = YOUR_CODE_HERE\n",
- "\n",
- " # Sigma_Y_hat = YOUR_CODE_HERE\n",
- " # std_y = YOUR_CODE_HERE\n",
- "\n",
- " # Sigma_e_hat = YOUR_CODE_HERE\n",
- " # std_e_hat = YOUR_CODE_HERE\n",
- "\n",
- " # SOLUTION:\n",
- " Sigma_X_hat = np.linalg.inv(A.T @ np.linalg.inv(Sigma_Y) @ A)\n",
- " x_hat = Sigma_X_hat @ A.T @ np.linalg.inv(Sigma_Y) @ y\n",
- " \n",
- " y_hat = A @ x_hat\n",
- "\n",
- " e_hat = y - y_hat\n",
- "\n",
- " Sigma_Y_hat = A @ Sigma_X_hat @ A.T\n",
- " std_y = np.sqrt(Sigma_Y_hat.diagonal())\n",
- "\n",
- " Sigma_e_hat = Sigma_Y - Sigma_Y_hat\n",
- " std_e_hat = np.sqrt(Sigma_e_hat.diagonal())\n",
- "\n",
- " d['Sigma_X_hat'] = Sigma_X_hat\n",
- " d['x_hat'] = x_hat\n",
- " d['y_hat'] = y_hat\n",
- " d['e_hat'] = e_hat\n",
- " d['Sigma_Y_hat'] = Sigma_Y_hat\n",
- " d['std_y'] = std_y\n",
- " d['Sigma_e_hat'] = Sigma_e_hat\n",
- " d['std_e_hat'] = std_e_hat\n",
- "\n",
- " return d"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f6c74941",
- "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 4.2:</b> \n",
- " \n",
- "Now that you have completed the function, apply it to our two models and then print values for the estimated parameters.\n",
- "</p>\n",
- "</div>\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "4a592ac1",
- "metadata": {
- "tags": []
- },
- "outputs": [],
- "source": [
- "model_insar = BLUE(model_insar)\n",
- "x_hat_insar = model_insar['x_hat']\n",
- "\n",
- "# YOUR_CODE_HERE\n",
- "\n",
- "# SOLUTION:\n",
- "print ('The InSAR-estimated offset is', np.round(x_hat_insar[0],3), 'mm')\n",
- "print ('The InSAR-estimated velocity is', np.round(x_hat_insar[1],4), 'mm/day')\n",
- "print ('The InSAR-estimated velocity is', np.round(x_hat_insar[1]*365,4), 'mm/year')\n",
- "print ('The InSAR-estimated GW factor is', np.round(x_hat_insar[2],3), '[-]\\n')\n",
- "\n",
- "model_gnss = BLUE(model_gnss)\n",
- "x_hat_gnss = model_gnss['x_hat']\n",
- "\n",
- "# YOUR_CODE_HERE\n",
- "\n",
- "# SOLUTION:\n",
- "print ('The GNSS-estimated offset is', np.round(x_hat_gnss[0],3), 'mm')\n",
- "print ('The GNSS-estimated velocity is', np.round(x_hat_gnss[1],4), 'mm/day')\n",
- "print ('The GNSS-estimated velocity is', np.round(x_hat_gnss[1]*365,4), 'mm/year')\n",
- "print ('The GNSS-estimated GW factor is', np.round(x_hat_gnss[2],3), '[-]')"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "bef2f3be",
- "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 4.3:</b> \n",
- "Do the values that you just estimated make sense? Explain, using quantitative results.\n",
- "\n",
- "<em>Hint: all you need to do is use the figures created above to verify that the parameter values printed above are reasonable (e.g., order of magnitude, units, etc).</em> \n",
- "</p>\n",
- "</div>\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "f6740791",
- "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",
- "As long as the velocity is negative and around -0.02 mm/day or -10 mm/yr it makes sense if you compare with what you see in the plots with observations. Since load is applied on soil layers we expect the road to subside. We also expect to see a positive value for the GW factor.\n",
- " \n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "65e42a43",
- "metadata": {
- "id": "65e42a43"
- },
- "source": [
- "## 5. Evaluate the precision\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "68f79fcb",
- "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 5:</b> \n",
- " \n",
- "What is the precision of the final estimates? \n",
- " \n",
- "Print the full covariance matrix of your estimates, and give an interpretation of the numbers in the covariance matrix.\n",
- "</p>\n",
- "</div>\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "835eefc8",
- "metadata": {
- "colab": {
- "base_uri": "https://localhost:8080/"
- },
- "executionInfo": {
- "elapsed": 16,
- "status": "ok",
- "timestamp": 1664699882186,
- "user": {
- "displayName": "C Yin",
- "userId": "14075875094781565898"
- },
- "user_tz": -120
- },
- "id": "835eefc8",
- "outputId": "ad47d53d-c147-4b25-fb03-b967a3bb6f96",
- "tags": []
- },
- "outputs": [],
- "source": [
- "Sigma_X_hat_insar = model_insar['Sigma_X_hat']\n",
- "\n",
- "# YOUR_CODE_HERE\n",
- "\n",
- "# SOLUTION:\n",
- "print ('Covariance matrix of estimated parameters (InSAR):')\n",
- "print (Sigma_X_hat_insar)\n",
- "print ('\\nThe standard deviation for the InSAR-estimated offset is', \n",
- " np.round(np.sqrt(Sigma_X_hat_insar[0,0]),3), 'mm')\n",
- "print ('The standard deviation for the InSAR-estimated velocity is', \n",
- " np.round(np.sqrt(Sigma_X_hat_insar[1,1]),4), 'mm/day')\n",
- "print ('The standard deviation for the InSAR-estimated GW factor is', \n",
- " np.round(np.sqrt(Sigma_X_hat_insar[2,2]),3), '[-]\\n')\n",
- "\n",
- "Sigma_X_hat_gnss = model_gnss['Sigma_X_hat']\n",
- "\n",
- "# YOUR_CODE_HERE\n",
- "\n",
- "# SOLUTION:\n",
- "print ('Covariance matrix of estimated parameters (GNSS):')\n",
- "print (Sigma_X_hat_gnss)\n",
- "print ('\\nThe standard deviation for the GNSS-estimated offset is', \n",
- " np.round(np.sqrt(Sigma_X_hat_gnss[0,0]),3), 'mm')\n",
- "print ('The standard deviation for the GNSS-estimated velocity is', \n",
- " np.round(np.sqrt(Sigma_X_hat_gnss[1,1]),4), 'mm/day')\n",
- "print ('The standard deviation for the GNSS-estimated GW factor is', \n",
- " np.round(np.sqrt(Sigma_X_hat_gnss[2,2]),3), '[-]')"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "8e62592d",
- "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",
- "As shown above, the standard deviations of the estimated parameters are equal to the square root of the diagonal elements. Compared with the estimated values, the standard deviations seem quite small, except for the estimated offsets. Meaning that the complete estimated model can be shifted up or down. \n",
- " \n",
- "The off-diagonal elements show the covariances between the estimated parameters, which are non-zeros since the estimates are all computed as function of the same vector of observations and the same model. A different value for the estimated velocity would imply a different value for the GW factor and offset. \n",
- " \n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "886efe26",
- "metadata": {
- "id": "886efe26"
- },
- "source": [
- "## 6. Present and reflect on estimation results"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "7ef41ce8",
- "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 6.1:</b> \n",
- " \n",
- "Complete the function below to help us compute the confidence intervals, then apply the function. Use a confidence interval of 96% in your analysis.\n",
- "\n",
- "<em>Hint: it can be used in exactly the same way as the <code>BLUE</code> function above, although it has one extra input.</em>\n",
- "</p>\n",
- "</div>\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "2711da12",
- "metadata": {},
- "outputs": [],
- "source": [
- "def get_CI(d, alpha):\n",
- " \"\"\"Compute the confidence intervals.\n",
- " \n",
- " Uses dict as input/output:\n",
- " - inputs defined from existing values in dict\n",
- " - outputs defined as new values in dict\n",
- " \"\"\"\n",
- "\n",
- " std_e_hat = d['std_e_hat']\n",
- " std_y = d['std_y']\n",
- "\n",
- " # k = YOUR_CODE_HERE\n",
- " # CI_y = YOUR_CODE_HERE\n",
- " # CI_res = YOUR_CODE_HERE\n",
- "\n",
- " # SOLUTION:\n",
- " k = norm.ppf(1 - 0.5*alpha)\n",
- " CI_y = k*std_y\n",
- " CI_res = k*std_e_hat\n",
- " CI_y_hat = k*np.sqrt(d['Sigma_Y_hat'].diagonal())\n",
- "\n",
- " d['alpha'] = alpha\n",
- " d['CI_y'] = CI_y\n",
- " d['CI_res'] = CI_res\n",
- " d['CI_Y_hat'] = CI_y_hat\n",
- "\n",
- " return d"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "d9a41ea5",
- "metadata": {},
- "outputs": [],
- "source": [
- "# model_insar = YOUR_CODE_HERE\n",
- "# model_gnss = YOUR_CODE_HERE\n",
- "\n",
- "# SOLUTION:\n",
- "model_insar = get_CI(model_insar, 0.04)\n",
- "model_gnss = get_CI(model_gnss, 0.04)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "53cf3663",
- "metadata": {},
- "source": [
- "At this point we have all the important results entered in our dictionary and we will be able to use the plots that have been written for you in the next Tasks. In case you would like to easily see all of the key-value pairs that have been added to the dictionary, you can run the cell below:"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "b3bb808e",
- "metadata": {},
- "outputs": [],
- "source": [
- "print(\"Keys and Values (type) for model_insar:\")\n",
- "for key, value in model_insar.items():\n",
- " print(f\"{key:16s} --> {type(value)}\")\n",
- "print(\"\\nKeys and Values (type) for model_gnss:\")\n",
- "for key, value in model_gnss.items():\n",
- " print(f\"{key:16s} --> {type(value)}\")"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "aaf72e41",
- "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 6.2:</b> \n",
- " \n",
- "Read the contents of file <code>functions.py</code> and identify what it is doing: you should be able to recognize that they use our model dictionary as an input and create three different figures. Note also that the function to create the figures have already been imported at the top of this notebook.\n",
- "\n",
- "Use the functions provided to visualize the results of our two models.\n",
- "</p>\n",
- "</div>\n"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "e8da1f6c-23de-4f80-a76c-5f0a9b4020b4",
- "metadata": {
- "id": "0491cc69"
- },
- "source": [
- "<div style=\"background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px\"> <p><strong>Note</strong>: remember that you will have to use the same function to look at <em>both</em> models when writing your interpretation in the Report.</p></div>"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "ec7c8bef",
- "metadata": {},
- "outputs": [],
- "source": [
- "# _, _ = plot_model(YOUR_CODE_HERE)\n",
- "\n",
- "# SOLUTION:\n",
- "_, _ = plot_model(model_insar)\n",
- "_, _ = plot_model(model_gnss)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "104d155d",
- "metadata": {},
- "outputs": [],
- "source": [
- "# _, _ = plot_residual(YOUR_CODE_HERE)\n",
- "\n",
- "# SOLUTION:\n",
- "_, _ = plot_residual(model_insar)\n",
- "_, _ = plot_residual(model_gnss)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "1dc93ce9",
- "metadata": {},
- "outputs": [],
- "source": [
- "# _, _ = plot_residual_histogram(YOUR_CODE_HERE)\n",
- "\n",
- "# SOLUTION:\n",
- "_, _ = plot_residual_histogram(model_insar)\n",
- "_, _ = plot_residual_histogram(model_gnss)"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "1ae74dd4",
- "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: the True Model}$\n",
- " \n",
- "The data used in this exercise was generated using Monte Carlo Simulation. It is added to the plots here to illustrate where and how our models differ (it is your job to interpret \"why\"). \n",
- "</p>\n",
- "</div>"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "113f3809",
- "metadata": {},
- "outputs": [],
- "source": [
- "k_true = 0.15\n",
- "R_true = -22 \n",
- "a_true = 180\n",
- "d0_true = 10\n",
- "\n",
- "disp_insar = (d0_true + R_true*(1 - np.exp(-days_insar/a_true)) +\n",
- " k_true*GW_at_InSAR_times)\n",
- "disp_gnss = (d0_true + R_true*(1 - np.exp(-days_gnss/a_true)) +\n",
- " k_true*GW_at_GNSS_times)\n",
- "\n",
- "plot_model(model_insar, alt_model=('True model', times_insar, disp_insar));\n",
- "plot_model(model_gnss, alt_model=('True model', times_gnss, disp_gnss));"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "aa877dd5",
- "metadata": {},
- "outputs": [],
- "source": [
- "import ipywidgets as widgets\n",
- "from ipywidgets import interact\n",
- "\n",
- "# Function to update the plot based on slider values\n",
- "def update_plot(x0, x1, x2):\n",
- " plt.figure(figsize=(15,5))\n",
- " for m in [model_gnss]: #[model_insar, model_gnss]:\n",
- " plt.plot(m['times'], m['y'], 'o', label=m['data_type'])\n",
- " plt.ylabel('Displacement [mm]')\n",
- " plt.xlabel('Time')\n",
- " \n",
- " y_fit = model_gnss['A'] @ [x0, x1, x2]\n",
- " if (x0 == 0) & (x1 == 0) & (x2 == 1):\n",
- " plt.plot(model_gnss['times'], y_fit, 'r', label='Groundwater data', linewidth=2)\n",
- " else:\n",
- " plt.plot(model_gnss['times'], y_fit, 'r', label='Fit (GNSS)', linewidth=2)\n",
- "\n",
- " W = np.linalg.inv(model_gnss['Sigma_Y'])\n",
- " ss_res = (model_gnss['y'] - y_fit).T @ W @ (model_gnss['y'] - y_fit)\n",
- " plt.title(f'Mean of squared residuals: {ss_res:.0f}')\n",
- " plt.grid()\n",
- " plt.legend()\n",
- " plt.show()\n",
- "\n",
- "# Create sliders for x0, x1, and x2\n",
- "x0_slider = widgets.FloatSlider(value=0, min=-10, max=10, step=0.1, description='x0')\n",
- "x1_slider = widgets.FloatSlider(value=0, min=-0.1, max=0.1, step=0.001, description='x1')\n",
- "x2_slider = widgets.FloatSlider(value=0, min=-1, max=1, step=0.01, description='x2')\n",
- "\n",
- "# Use interact to create the interactive plot\n",
- "interact(update_plot, x0=x0_slider, x1=x1_slider, x2=x2_slider)\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "d8b31540",
- "metadata": {},
- "outputs": [],
- "source": [
- "xhat_slider_plot(model_gnss['A'], model_gnss['y'], model_gnss['times'], model_gnss['Sigma_Y'])"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "3203d779",
- "metadata": {},
- "source": [
- "**End of notebook.**\n",
- "<h2 style=\"height: 60px\">\n",
- "</h2>\n",
- "<h3 style=\"position: absolute; display: flex; flex-grow: 0; flex-shrink: 0; flex-direction: row-reverse; bottom: 60px; right: 50px; margin: 0; border: 0\">\n",
- " <style>\n",
- " .markdown {width:100%; position: relative}\n",
- " article { position: relative }\n",
- " </style>\n",
- " <a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\">\n",
- " <img alt=\"Creative Commons License\" style=\"border-width:; width:88px; height:auto; padding-top:10px\" src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\" />\n",
- " </a>\n",
- " <a rel=\"TU Delft\" href=\"https://www.tudelft.nl/en/ceg\">\n",
- " <img alt=\"TU Delft\" style=\"border-width:0; width:100px; height:auto; padding-bottom:0px\" src=\"https://gitlab.tudelft.nl/mude/public/-/raw/main/tu-logo/TU_P1_full-color.png\" />\n",
- " </a>\n",
- " <a rel=\"MUDE\" href=\"http://mude.citg.tudelft.nl/\">\n",
- " <img alt=\"MUDE\" style=\"border-width:0; width:100px; height:auto; padding-bottom:0px\" src=\"https://gitlab.tudelft.nl/mude/public/-/raw/main/mude-logo/MUDE_Logo-small.png\" />\n",
- " </a>\n",
- " \n",
- "</h3>\n",
- "<span style=\"font-size: 75%\">\n",
- "© Copyright 2024 <a rel=\"MUDE\" href=\"http://mude.citg.tudelft.nl/\">MUDE</a> TU Delft. This work is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\">CC BY 4.0 License</a>."
- ]
- }
- ],
- "metadata": {
- "colab": {
- "collapsed_sections": [],
- "provenance": []
- },
- "kernelspec": {
- "display_name": "Python 3 (ipykernel)",
- "language": "python",
- "name": "python3"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 3
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython3",
- "version": "3.12.4"
- },
- "vscode": {
- "interpreter": {
- "hash": "b932539803d9742d977fbe9ca28a706a58466dc5c9df0c7af4e41c76d82e5a85"
- }
- },
- "widgets": {
- "application/vnd.jupyter.widget-state+json": {
- "state": {},
- "version_major": 2,
- "version_minor": 0
- }
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/synced_files/test.md b/synced_files/test.md
deleted file mode 100644
index 7be68f6516bf10d3df9fadd25c4b6e6c1ad62078..0000000000000000000000000000000000000000
--- a/synced_files/test.md
+++ /dev/null
@@ -1,975 +0,0 @@
----
-jupyter:
- jupytext:
- text_representation:
- extension: .md
- format_name: markdown
- format_version: '1.3'
- jupytext_version: 1.16.4
- kernelspec:
- display_name: Python 3 (ipykernel)
- language: python
- name: python3
----
-
-# GA 1.3: test
-
-<h1 style="position: absolute; display: flex; flex-grow: 0; flex-shrink: 0; flex-direction: row-reverse; top: 60px;right: 30px; margin: 0; border: 0">
- <style>
- .markdown {width:100%; position: relative}
- article { position: relative }
- </style>
- <img src="https://gitlab.tudelft.nl/mude/public/-/raw/main/tu-logo/TU_P1_full-color.png" style="width:100px" />
- <img src="https://gitlab.tudelft.nl/mude/public/-/raw/main/mude-logo/MUDE_Logo-small.png" style="width:100px" />
-</h1>
-<h2 style="height: 10px">
-</h2>
-
-*[CEGM1000 MUDE](http://mude.citg.tudelft.nl/): Week 1.3. Due: Friday, September 20, 2024.*
-
-
-<div style="background-color:#ffa6a6; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"><p><b>Note:</b> don't forget to read the "Assignment Context" section of the README, it contains important information to understand this analysis.</p></div>
-
-```python
-import numpy as np
-from scipy import interpolate
-from scipy.stats import norm
-import pandas as pd
-import matplotlib.pyplot as plt
-
-
-from functions import *
-
-np.set_printoptions(precision=3)
-```
-
-## Part 0: Dictionary Review
-
-As described above, several functions in this assignment require the use of a Python dictionary to make it easier to keep track of important data, variables and results for the various _models_ we will be constructing and validating.
-
-_It may be useful to revisit PA 1.1, where there was a brief infroduction to dictionaires. That PA contains all the dictionary info you need for GA 1.3. A [read-only copy is here](https://mude.citg.tudelft.nl/2024/files/Week_1_1/PA_1_1_Catch_Them_All.html) and [the source code (notebook) is here](https://gitlab.tudelft.nl/mude/2024-week-1-1)._
-
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-
-$\textbf{Task 0.1}$
-
-Read and run the cell below to make sure you remember how to use a dictionary.
-
-Modify the function to print some of the other key-value pairs of the dictionary.
-
-<em>It may also be useful to use the cell below when working on later tasks in this assignment.</em>
-
-</p>
-</div>
-
-```python
-my_dictionary = {'key1': 'value1',
- 'key2': 'value2',
- 'name': 'Dictionary Example',
- 'a_list': [1, 2, 3],
- 'an_array': np.array([1, 2, 3]),
- 'a_string': 'hello'
- }
-
-def function_that_uses_my_dictionary(d):
- print(d['key1'])
-
- # SOLUTION:
- print(d['name'])
- print(d['a_list'])
- print(d['an_array'])
- print(d['a_string'])
-
- if 'new_key' in d:
- print('new_key exists and has value:', d['new_key'])
- return
-
-function_that_uses_my_dictionary(my_dictionary)
-```
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-
-$\textbf{Task 0.2}$
-
-Test your knowledge by adding a new key <code>new_key</code> and then executing the function to print the value.
-
-</p>
-</div>
-
-```python
-# YOUR_CODE_HERE
-# function_that_uses_my_dictionary(my_dictionary)
-
-# SOLUTION:
-my_dictionary['new_key'] = 'new_value'
-function_that_uses_my_dictionary(my_dictionary)
-```
-
-<!-- #region id="160d6250" -->
-## Task 1: Preparing the data
-
-Within this assignment you will work with two types of data: InSAR data and GNSS data. The cell below will load the data and visualize the observed displacements time. In this task we use the package `pandas`, which is really useful for handling time series. We will learn how to use it later in the quarter; for now, you only need to recognize that it imports the data as a `dataframe` object, which we then convert into a numpy array using the code below.
-<!-- #endregion -->
-
-<div style="background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px"> <p>Tip: note that we have converted all observations to millimeters.</p></div>
-
-```python
-gnss = pd.read_csv('./data/gnss_observations.csv')
-times_gnss = pd.to_datetime(gnss['times'])
-y_gnss = (gnss['observations[m]']).to_numpy()*1000
-
-insar = pd.read_csv('./data/insar_observations.csv')
-times_insar = pd.to_datetime(insar['times'])
-y_insar = (insar['observations[m]']).to_numpy()*1000
-
-gw = pd.read_csv('./data/groundwater_levels.csv')
-times_gw = pd.to_datetime(gw['times'])
-y_gw = (gw['observations[mm]']).to_numpy()
-```
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 1.1:</b>
-
-Once you have used the cell above to import the data, investigate the data sets using the code cell below. Then provide some relevant summary information in the Markdown cell.
-
-<em>Hint: at the least, you should be able to tell how many data points are in each data set and get an understanding of the mean and standard deviation of each. Make sure you compare the different datasets and use consistent units.</em>
-
-</p>
-</div>
-
-<!-- #region id="0491cc69" -->
-<div style="background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px"> <p>The code below gives some examples of the quantitative and qualitative ways you could have looked at the data. It is more than you were expected to do; the important thing is that you showed the ability to learn something about the data and describe aspects that are relevant to our problem. We use a dictionary to easily access the different data series using their names, which are entered as the dictionary keys (also not expected of you, but it's hopefully fun to learn useful tricks).</div>
-<!-- #endregion -->
-
-```python
-# YOUR_CODE_HERE
-
-# SOLUTION:
-data_list = ['y_gnss', 'y_insar', 'y_gw']
-data_dict = {data_list[0]: y_gnss,
- data_list[1]: y_insar,
- data_list[2]: y_gw}
-def print_summary(data):
- '''Summarize an array with simple print statements.'''
- print('Minimum = ', data.min())
- print('Maximum = ', data.max())
- print('Mean = ', data.mean())
- print('Std dev = ', data.std())
- print('Shape = ', data.shape)
- print('First value = ', data[0])
- print('Last value = ', data[-1])
- print('\n')
-
-for item in data_list:
- print('Summary for array: ', item)
- print('------------------------------------------------')
- print_summary(data_dict[item])
-```
-
-```python
-# SOLUTION:
-times_dict = {data_list[0]: times_gnss,
- data_list[1]: times_insar,
- data_list[2]: times_gw}
-def plot_data(times, data, label):
- plt.figure(figsize=(15,4))
- plt.plot(times, data, 'co', mec='black')
- plt.title(label)
- plt.xlabel('Times')
- plt.ylabel('Data [mm]')
- plt.show()
-
-plt.figure(figsize=(15,4))
-for i in range(3):
- plot_data(times_dict[data_list[i]],
- data_dict[data_list[i]],
- data_list[i])
-```
-
-<div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-
-$\textbf{Solution}$
-
-There are a lot more GNSS data points than InSAR or groundwater. The GNSS observations also have more noise, and what seem to be outliers. In this case the mean and standard deviation do not mean much, because there is clearly a trend with time. We can at least confirm that the time periods of measurements overlap, although the intervals between measurements is certainly not uniform (note that you don't need to do anything with the times, since they are pandas time series and we have not covered them yet).
-</p>
-</div>
-
-
-You may have noticed that the groundwater data is available for different times than the GNSS and InSAR data. You will therefore have to *interpolate* the data to the same times for a further analysis. You can use the SciPy function ```interpolate.interp1d``` (read its [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html)).
-
-The cells below do the following:
-1. Define a function to convert the time unit
-2. Convert the time stamps for all data
-3. Use `interp1d` to interpolate the groundwater measurements at the time of the satellite measurements
-
-```python
-def to_days_years(times):
- '''Convert the observation times to days and years.'''
-
- times_datetime = pd.to_datetime(times)
- time_diff = (times_datetime - times_datetime[0])
- days_diff = (time_diff / np.timedelta64(1,'D')).astype(int)
-
- days = days_diff.to_numpy()
- years = days/365
-
- return days, years
-```
-
-```python
-days_gnss, years_gnss = to_days_years(times_gnss)
-days_insar, years_insar = to_days_years(times_insar)
-days_gw, years_gw = to_days_years(times_gw)
-
-interp = interpolate.interp1d(days_gw, y_gw)
-
-GW_at_GNSS_times = interp(days_gnss)
-GW_at_InSAR_times = interp(days_insar)
-```
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 1.2:</b>
-
-Answer/complete the code and Markdown cells below:
-<ol>
- <li>What is <code>interp</code>? (what kind of object is it, and how does it work?)</li>
- <li>How did the groundwater observation array change? Be quantitative. </li>
-</ol>
-
-</p>
-</div>
-
-```python
-# YOUR_CODE_HERE
-
-# SOLUTION:
-print('array size of GW_at_GNSS_times', len(GW_at_GNSS_times))
-print('array size of GW_at_InSAR_times', len(GW_at_InSAR_times))
-print('array size of GW before interpolation', len(y_gw))
-
-print('\nFirst values of times_gw:')
-print(times_gw[0:2])
-print('\nFirst values of y_gw:')
-print(y_gw[0:2])
-print('\nFirst values of times_gnss:')
-print(times_gnss[0:2])
-print('\nFirst values of GW_at_GNSS_times:')
-print(GW_at_GNSS_times[0:2])
-```
-
-**Write your answer in this Markdown cell.**
-
-
-<div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-
-$\textbf{Solution}$
-<ol>
- <li><code>interp</code> is a function that will return a value (gw level) for the input(s) (date(s)). The interpolated value is found by linearly interpolating between the two nearest times in the gw observations.</li>
- <li>The observation arrays of <code>GW_at_GNSS_times</code> and <code>GW_at_INSAR_times</code> changed in size to match the size of the GNSS and InSAR observations, respectively.</li>
-</ol>
-</p>
-</div>
-
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 1.3:</b>
-
-Create a single plot to compare observed displacement for the GNSS and InSAR data sets.
-
-</p>
-</div>
-
-```python
-# plt.figure(figsize=(15,5))
-# plt.plot(YOUR_CODE_HERE, YOUR_CODE_HERE,
-# 'o', mec='black', label = 'GNSS')
-# plt.plot(YOUR_CODE_HERE, YOUR_CODE_HERE,
-# 'o', mec='black', label = 'InSAR')
-# plt.legend()
-# plt.ylabel('Displacement [mm]')
-# plt.xlabel('Time')
-# plt.show()
-
-# SOLUTION:
-plt.figure(figsize=(15,5))
-plt.plot(times_gnss, y_gnss, 'o', mec='black', label = 'GNSS')
-plt.plot(times_insar, y_insar, 'o', mec='black', label = 'InSAR')
-plt.legend()
-plt.ylabel('Displacement [mm]')
-plt.xlabel('Time')
-plt.show()
-```
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 1.4:</b>
-Describe the datasets based on the figure above and your observations from the previous tasks. What kind of deformation do you see? And what are the differences between both datasets? Be quantitative.
-</p>
-</div>
-
-
-<div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-
-$\textbf{Solution}$
-
-The points obviously show subsidence, the displacement shows a similar pattern for both datasets. The GNSS data is much noisier than InSAR (range is around 60 mm versus only a few mm), but has a higher sampling rate. Also there seem to be more outliers in the GNSS data compared to InSAR, especially at the start of the observation period. InSAR has only observations every 6 days but is less noisy.
-</p>
-</div>
-
-
-Before we move on, it is time to do a little bit of housekeeping.
-
-Have you found it confusing to keep track of two sets of variables---one for each data type? Let's use a dictionary to store relevant information about each model. We will use this in the plotting functions for this task (and again next week), so make sure you take the time to see what is happening. Review also Part 0 at the top of this notebook if you need a refresher on dictionaries.
-
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 1.5:</b>
-
-Run the cell below to define a dictionary for storing information about the two (future) models.
-
-</p>
-</div>
-
-
-```python
-model_insar = {'data_type': 'InSAR',
- 'y':y_insar,
- 'times':times_insar,
- 'groundwater': GW_at_InSAR_times
- }
-
-model_gnss = {'data_type': 'GNSS',
- 'y':y_gnss,
- 'times':times_gnss,
- 'groundwater': GW_at_GNSS_times
- }
-```
-
-<!-- #region id="76c9115b" -->
-## Task 2: Set-up linear functional model
-
-We want to investigate how we could model the observed displacements of the road. Because the road is built in the Green Heart we expect that the observed displacements are related to the groundwater level. Furthermore, we assume that the displacements can be modeled using a constant velocity. The model is defined as
-$$
-d = d_0 + vt + k \ \textrm{GW},
-$$
-where $d$ is the displacement, $t$ is time and $\textrm{GW}$ is the groundwater level (that we assume to be deterministic).
-
-Therefore, the model has 3 unknowns:
-1. $d_0$, as the initial displacement at $t_0$;
-2. $v$, as the displacement velocity;
-3. $k$, as the 'groundwater factor', which can be seen as the response of the soil to changes in the groundwater level.
-
-
-As a group you will construct the **functional model** that is defined as
-$$
-\mathbb{E}(Y) = \mathrm{A x}.
-$$
-
-
-<!-- #endregion -->
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 2.1:</b>
-
-Construct the design matrix $A$ (for both InSAR and GNSS observations), then show the first 5 observations and confirm the dimensions of $A$.
-</p>
-</div>
-
-```python
-# YOUR_CODE_HERE
-
-# SOLUTION:
-A_insar = np.ones((len(times_insar), 3))
-A_insar[:,1] = days_insar
-A_insar[:,2] = GW_at_InSAR_times
-
-print ('The first 5 rows of the A matrix (InSAR) are:')
-print (A_insar[0:5, :])
-
-print ('The first 5 observations [mm] of y_insar are:')
-print (y_insar[0:5])
-
-m_insar = np.shape(A_insar)[0]
-n_insar = np.shape(A_insar)[1]
-print(f'm = {m_insar} and n = {n_insar}')
-```
-
-```python
-# YOUR_CODE_HERE
-
-# SOLUTION:
-A_gnss = np.ones((len(times_gnss), 3))
-A_gnss[:,1] = days_gnss
-A_gnss[:,2] = GW_at_GNSS_times
-
-print ('The first 5 rows of the A matrix (GNSS) are:')
-print (A_gnss[0:5, :])
-
-print ('\nThe first 5 observations [mm] of y_gnss are:')
-print (y_gnss[0:5])
-
-m_gnss = np.shape(A_gnss)[0]
-n_gnss = np.shape(A_gnss)[1]
-print(f'm = {m_gnss} and n = {n_gnss}')
-```
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-
-$\textbf{Task 2.2}$
-
-Answer the following questions:
-
-- What is the dimension of the observables' vector $Y$?
-- What are the unknowns of the functional model?
-- What is the redundancy for this model?
-
-</p>
-</div>
-
-
-
-<div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-
-$\textbf{Solution}$
-
-For InSAR:
-<ol>
- <li>The number of observations is 61.</li>
- <li>The number of unknowns is 3.</li>
- <li>The redundancy is 58.</li>
-</ol>
-
-For GNSS:
-<ol>
- <li>The number of observations is 730.</li>
- <li>The number of unknowns is 3.</li>
- <li>The redundancy is 727.</li>
-</ol>
-</p>
-</div>
-
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 2.3:</b>
-
-Add the A matrix to the dictionaries for each model. This will be used to plot results later in the notebook.
-
-</p>
-</div>
-
-
-```python
-# model_insar['A'] = YOUR_CODE_HERE
-# model_gnss['A'] = YOUR_CODE_HERE
-
-# SOLUTION:
-model_insar['A'] = A_insar
-model_gnss['A'] = A_gnss
-
-print("Keys and Values (type) for model_insar:")
-for key, value in model_insar.items():
- print(f"{key:16s} --> {type(value)}")
-print("\nKeys and Values (type) for model_gnss:")
-for key, value in model_gnss.items():
- print(f"{key:16s} --> {type(value)}")
-```
-
-<!-- #region id="9325d32b" -->
-## 3. Set-up stochastic model
-
-We will use the Best Linear Unbiased Estimator (BLUE) to solve for the unknown parameters. Therefore we also need a stochastic model, which is defined as
-$$
-\mathbb{D}(Y) = \Sigma_{Y}.
-$$
-where $\Sigma_{Y}$ is the covariance matrix of the observables' vector.
-
-
-
-<!-- #endregion -->
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 3.1:</b>
-
-Construct the covariance matrix for each type of data and assume that
-
-- the observables are independent
-
-- the observables are normally distributed
-
-- the observables' standard deviation is
-
- - $\sigma_\textrm{InSAR} = 2$ mm
- - $\sigma_\textrm{GNSS} = 15$ mm
-
-</p>
-</div>
-
-
-```python colab={"base_uri": "https://localhost:8080/"} executionInfo={"elapsed": 40, "status": "ok", "timestamp": 1664699881875, "user": {"displayName": "C Yin", "userId": "14075875094781565898"}, "user_tz": -120} id="163acdb3" outputId="8bc99da3-a61e-4a25-8a54-c90f33299a57"
-# YOUR_CODE_HERE
-
-# SOLUTION:
-std_insar = 2 #mm
-
-Sigma_Y_insar = np.identity(len(times_insar))*std_insar**2
-
-print ('Sigma_Y (InSAR) is defined as:')
-print (Sigma_Y_insar)
-```
-
-```python
-# YOUR_CODE_HERE
-
-# SOLUTION:
-std_gnss = 15 #mm (corrected from original value of 5 mm)
-
-Sigma_Y_gnss = np.identity(len(times_gnss))*std_gnss**2
-
-print ('\nSigma_Y (GNSS) is defined as:')
-print (Sigma_Y_gnss)
-```
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-
-$\textbf{Task 3.2}$
-
-Answer the following questions:
-
-- What information is contained in the covariance matrix?
-- How do you implement the assumption that all observations are independent?
-- What is the dimension of $\Sigma_{Y}$?
-- How do you create $\Sigma_{Y}$?
-
-</p>
-</div>
-
-
-_Write your answer in this cell._
-
-
-<div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-
-$\textbf{Solution}$
-
-- The covariance matrix contains information on the quality of the observations, where an entry on the diagonal represents the variance of one observation at a particular epoch. If there is an indication that for instance the quality for a particular time interval differs, different $\sigma$ values can be put in the stochastic model for these epochs.
-- The off-diagonal terms in the matrix are related to the correlation between observations at different epochs, where a zero value on the off-diagonal indicates zero correlation.
-- The dimension of the matrix is 61x61 for InSAR and 730x730 for GNSS.
-- See code.
-</p>
-</div>
-
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 3.3:</b>
-
-Add <code>Sigma_Y</code> to the dictionaries for each model.
-
-</p>
-</div>
-
-
-```python
-# model_insar['Sigma_Y] = YOUR_CODE_HERE
-# model_gnss['Sigma_Y'] = YOUR_CODE_HERE
-
-# SOLUTION:
-model_insar['Sigma_Y'] = Sigma_Y_insar
-model_gnss['Sigma_Y'] = Sigma_Y_gnss
-```
-
-<!-- #region id="09e965bf" -->
-## 4. Apply best linear unbiased estimation
-
-
-<!-- #endregion -->
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 4.1:</b>
-
-Write a function to apply BLUE in the cell below and use the function to estimate the unknowns for the model using the data.
-
-Compute the modeled displacements ($\hat{\mathrm{y}}$), and corresponding residuals ($\hat{\mathrm{\epsilon}}$), as well as associated values (as requested by the blank code lines).
-</p>
-</div>
-
-
-
-
-<div style="background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px"> <p><strong>Note on code implementation</strong>: you'll see that the functions in this assignment use a dictionary; this greatly reduces the number of input/output variables needed in a function. However, it can make the code inside the function more difficult to read due to the key syntax (e.g., <code>dict['variable_1']</code> versus <code>variable
-_1</code>). To make this assignment easier for you to implement we have split these functions into three parts: 1) define variables from the dictionary, 2) perform analysis, 3) add results to the dictionary. Note that this is not the most efficient way to write this code; it is done here specifically for clarity and to help you focus on writing the equations properly and understanding the meaning of each term.</p></div>
-
-```python
-def BLUE(d):
- """Calculate the Best Linear Unbiased Estimator
-
- Uses dict as input/output:
- - inputs defined from existing values in dict
- - outputs defined as new values in dict
- """
-
- y = d['y']
- A = d['A']
- Sigma_Y = d['Sigma_Y']
-
- # Sigma_X_hat = YOUR_CODE_HERE
- # x_hat = YOUR_CODE_HERE
-
- # y_hat = YOUR_CODE_HERE
-
- # e_hat = YOUR_CODE_HERE
-
- # Sigma_Y_hat = YOUR_CODE_HERE
- # std_y = YOUR_CODE_HERE
-
- # Sigma_e_hat = YOUR_CODE_HERE
- # std_e_hat = YOUR_CODE_HERE
-
- # SOLUTION:
- Sigma_X_hat = np.linalg.inv(A.T @ np.linalg.inv(Sigma_Y) @ A)
- x_hat = Sigma_X_hat @ A.T @ np.linalg.inv(Sigma_Y) @ y
-
- y_hat = A @ x_hat
-
- e_hat = y - y_hat
-
- Sigma_Y_hat = A @ Sigma_X_hat @ A.T
- std_y = np.sqrt(Sigma_Y_hat.diagonal())
-
- Sigma_e_hat = Sigma_Y - Sigma_Y_hat
- std_e_hat = np.sqrt(Sigma_e_hat.diagonal())
-
- d['Sigma_X_hat'] = Sigma_X_hat
- d['x_hat'] = x_hat
- d['y_hat'] = y_hat
- d['e_hat'] = e_hat
- d['Sigma_Y_hat'] = Sigma_Y_hat
- d['std_y'] = std_y
- d['Sigma_e_hat'] = Sigma_e_hat
- d['std_e_hat'] = std_e_hat
-
- return d
-```
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 4.2:</b>
-
-Now that you have completed the function, apply it to our two models and then print values for the estimated parameters.
-</p>
-</div>
-
-
-
-```python
-model_insar = BLUE(model_insar)
-x_hat_insar = model_insar['x_hat']
-
-# YOUR_CODE_HERE
-
-# SOLUTION:
-print ('The InSAR-estimated offset is', np.round(x_hat_insar[0],3), 'mm')
-print ('The InSAR-estimated velocity is', np.round(x_hat_insar[1],4), 'mm/day')
-print ('The InSAR-estimated velocity is', np.round(x_hat_insar[1]*365,4), 'mm/year')
-print ('The InSAR-estimated GW factor is', np.round(x_hat_insar[2],3), '[-]\n')
-
-model_gnss = BLUE(model_gnss)
-x_hat_gnss = model_gnss['x_hat']
-
-# YOUR_CODE_HERE
-
-# SOLUTION:
-print ('The GNSS-estimated offset is', np.round(x_hat_gnss[0],3), 'mm')
-print ('The GNSS-estimated velocity is', np.round(x_hat_gnss[1],4), 'mm/day')
-print ('The GNSS-estimated velocity is', np.round(x_hat_gnss[1]*365,4), 'mm/year')
-print ('The GNSS-estimated GW factor is', np.round(x_hat_gnss[2],3), '[-]')
-```
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 4.3:</b>
-Do the values that you just estimated make sense? Explain, using quantitative results.
-
-<em>Hint: all you need to do is use the figures created above to verify that the parameter values printed above are reasonable (e.g., order of magnitude, units, etc).</em>
-</p>
-</div>
-
-
-
-
-<div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-
-$\textbf{Solution}$
-
-As long as the velocity is negative and around -0.02 mm/day or -10 mm/yr it makes sense if you compare with what you see in the plots with observations. Since load is applied on soil layers we expect the road to subside. We also expect to see a positive value for the GW factor.
-
-</p>
-</div>
-
-<!-- #region id="65e42a43" -->
-## 5. Evaluate the precision
-
-
-<!-- #endregion -->
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 5:</b>
-
-What is the precision of the final estimates?
-
-Print the full covariance matrix of your estimates, and give an interpretation of the numbers in the covariance matrix.
-</p>
-</div>
-
-
-
-```python colab={"base_uri": "https://localhost:8080/"} executionInfo={"elapsed": 16, "status": "ok", "timestamp": 1664699882186, "user": {"displayName": "C Yin", "userId": "14075875094781565898"}, "user_tz": -120} id="835eefc8" outputId="ad47d53d-c147-4b25-fb03-b967a3bb6f96"
-Sigma_X_hat_insar = model_insar['Sigma_X_hat']
-
-# YOUR_CODE_HERE
-
-# SOLUTION:
-print ('Covariance matrix of estimated parameters (InSAR):')
-print (Sigma_X_hat_insar)
-print ('\nThe standard deviation for the InSAR-estimated offset is',
- np.round(np.sqrt(Sigma_X_hat_insar[0,0]),3), 'mm')
-print ('The standard deviation for the InSAR-estimated velocity is',
- np.round(np.sqrt(Sigma_X_hat_insar[1,1]),4), 'mm/day')
-print ('The standard deviation for the InSAR-estimated GW factor is',
- np.round(np.sqrt(Sigma_X_hat_insar[2,2]),3), '[-]\n')
-
-Sigma_X_hat_gnss = model_gnss['Sigma_X_hat']
-
-# YOUR_CODE_HERE
-
-# SOLUTION:
-print ('Covariance matrix of estimated parameters (GNSS):')
-print (Sigma_X_hat_gnss)
-print ('\nThe standard deviation for the GNSS-estimated offset is',
- np.round(np.sqrt(Sigma_X_hat_gnss[0,0]),3), 'mm')
-print ('The standard deviation for the GNSS-estimated velocity is',
- np.round(np.sqrt(Sigma_X_hat_gnss[1,1]),4), 'mm/day')
-print ('The standard deviation for the GNSS-estimated GW factor is',
- np.round(np.sqrt(Sigma_X_hat_gnss[2,2]),3), '[-]')
-```
-
-<div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-
-$\textbf{Solution}$
-
-As shown above, the standard deviations of the estimated parameters are equal to the square root of the diagonal elements. Compared with the estimated values, the standard deviations seem quite small, except for the estimated offsets. Meaning that the complete estimated model can be shifted up or down.
-
-The off-diagonal elements show the covariances between the estimated parameters, which are non-zeros since the estimates are all computed as function of the same vector of observations and the same model. A different value for the estimated velocity would imply a different value for the GW factor and offset.
-
-</p>
-</div>
-
-<!-- #region id="886efe26" -->
-## 6. Present and reflect on estimation results
-<!-- #endregion -->
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 6.1:</b>
-
-Complete the function below to help us compute the confidence intervals, then apply the function. Use a confidence interval of 96% in your analysis.
-
-<em>Hint: it can be used in exactly the same way as the <code>BLUE</code> function above, although it has one extra input.</em>
-</p>
-</div>
-
-
-```python
-def get_CI(d, alpha):
- """Compute the confidence intervals.
-
- Uses dict as input/output:
- - inputs defined from existing values in dict
- - outputs defined as new values in dict
- """
-
- std_e_hat = d['std_e_hat']
- std_y = d['std_y']
-
- # k = YOUR_CODE_HERE
- # CI_y = YOUR_CODE_HERE
- # CI_res = YOUR_CODE_HERE
-
- # SOLUTION:
- k = norm.ppf(1 - 0.5*alpha)
- CI_y = k*std_y
- CI_res = k*std_e_hat
- CI_y_hat = k*np.sqrt(d['Sigma_Y_hat'].diagonal())
-
- d['alpha'] = alpha
- d['CI_y'] = CI_y
- d['CI_res'] = CI_res
- d['CI_Y_hat'] = CI_y_hat
-
- return d
-```
-
-```python
-# model_insar = YOUR_CODE_HERE
-# model_gnss = YOUR_CODE_HERE
-
-# SOLUTION:
-model_insar = get_CI(model_insar, 0.04)
-model_gnss = get_CI(model_gnss, 0.04)
-```
-
-At this point we have all the important results entered in our dictionary and we will be able to use the plots that have been written for you in the next Tasks. In case you would like to easily see all of the key-value pairs that have been added to the dictionary, you can run the cell below:
-
-```python
-print("Keys and Values (type) for model_insar:")
-for key, value in model_insar.items():
- print(f"{key:16s} --> {type(value)}")
-print("\nKeys and Values (type) for model_gnss:")
-for key, value in model_gnss.items():
- print(f"{key:16s} --> {type(value)}")
-```
-
-<div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-<b>Task 6.2:</b>
-
-Read the contents of file <code>functions.py</code> and identify what it is doing: you should be able to recognize that they use our model dictionary as an input and create three different figures. Note also that the function to create the figures have already been imported at the top of this notebook.
-
-Use the functions provided to visualize the results of our two models.
-</p>
-</div>
-
-
-<!-- #region id="0491cc69" -->
-<div style="background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px"> <p><strong>Note</strong>: remember that you will have to use the same function to look at <em>both</em> models when writing your interpretation in the Report.</p></div>
-<!-- #endregion -->
-
-```python
-# _, _ = plot_model(YOUR_CODE_HERE)
-
-# SOLUTION:
-_, _ = plot_model(model_insar)
-_, _ = plot_model(model_gnss)
-```
-
-```python
-# _, _ = plot_residual(YOUR_CODE_HERE)
-
-# SOLUTION:
-_, _ = plot_residual(model_insar)
-_, _ = plot_residual(model_gnss)
-```
-
-```python
-# _, _ = plot_residual_histogram(YOUR_CODE_HERE)
-
-# SOLUTION:
-_, _ = plot_residual_histogram(model_insar)
-_, _ = plot_residual_histogram(model_gnss)
-```
-
-<div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-<p>
-
-$\textbf{Solution: the True Model}$
-
-The data used in this exercise was generated using Monte Carlo Simulation. It is added to the plots here to illustrate where and how our models differ (it is your job to interpret "why").
-</p>
-</div>
-
-```python
-k_true = 0.15
-R_true = -22
-a_true = 180
-d0_true = 10
-
-disp_insar = (d0_true + R_true*(1 - np.exp(-days_insar/a_true)) +
- k_true*GW_at_InSAR_times)
-disp_gnss = (d0_true + R_true*(1 - np.exp(-days_gnss/a_true)) +
- k_true*GW_at_GNSS_times)
-
-plot_model(model_insar, alt_model=('True model', times_insar, disp_insar));
-plot_model(model_gnss, alt_model=('True model', times_gnss, disp_gnss));
-```
-
-```python
-import ipywidgets as widgets
-from ipywidgets import interact
-
-# Function to update the plot based on slider values
-def update_plot(x0, x1, x2):
- plt.figure(figsize=(15,5))
- for m in [model_gnss]: #[model_insar, model_gnss]:
- plt.plot(m['times'], m['y'], 'o', label=m['data_type'])
- plt.ylabel('Displacement [mm]')
- plt.xlabel('Time')
-
- y_fit = model_gnss['A'] @ [x0, x1, x2]
- if (x0 == 0) & (x1 == 0) & (x2 == 1):
- plt.plot(model_gnss['times'], y_fit, 'r', label='Groundwater data', linewidth=2)
- else:
- plt.plot(model_gnss['times'], y_fit, 'r', label='Fit (GNSS)', linewidth=2)
-
- W = np.linalg.inv(model_gnss['Sigma_Y'])
- ss_res = (model_gnss['y'] - y_fit).T @ W @ (model_gnss['y'] - y_fit)
- plt.title(f'Mean of squared residuals: {ss_res:.0f}')
- plt.grid()
- plt.legend()
- plt.show()
-
-# Create sliders for x0, x1, and x2
-x0_slider = widgets.FloatSlider(value=0, min=-10, max=10, step=0.1, description='x0')
-x1_slider = widgets.FloatSlider(value=0, min=-0.1, max=0.1, step=0.001, description='x1')
-x2_slider = widgets.FloatSlider(value=0, min=-1, max=1, step=0.01, description='x2')
-
-# Use interact to create the interactive plot
-interact(update_plot, x0=x0_slider, x1=x1_slider, x2=x2_slider)
-
-```
-
-```python
-xhat_slider_plot(model_gnss['A'], model_gnss['y'], model_gnss['times'], model_gnss['Sigma_Y'])
-```
-
-**End of notebook.**
-<h2 style="height: 60px">
-</h2>
-<h3 style="position: absolute; display: flex; flex-grow: 0; flex-shrink: 0; flex-direction: row-reverse; bottom: 60px; right: 50px; margin: 0; border: 0">
- <style>
- .markdown {width:100%; position: relative}
- article { position: relative }
- </style>
- <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">
- <img alt="Creative Commons License" style="border-width:; width:88px; height:auto; padding-top:10px" src="https://i.creativecommons.org/l/by/4.0/88x31.png" />
- </a>
- <a rel="TU Delft" href="https://www.tudelft.nl/en/ceg">
- <img alt="TU Delft" style="border-width:0; width:100px; height:auto; padding-bottom:0px" src="https://gitlab.tudelft.nl/mude/public/-/raw/main/tu-logo/TU_P1_full-color.png" />
- </a>
- <a rel="MUDE" href="http://mude.citg.tudelft.nl/">
- <img alt="MUDE" style="border-width:0; width:100px; height:auto; padding-bottom:0px" src="https://gitlab.tudelft.nl/mude/public/-/raw/main/mude-logo/MUDE_Logo-small.png" />
- </a>
-
-</h3>
-<span style="font-size: 75%">
-© Copyright 2024 <a rel="MUDE" href="http://mude.citg.tudelft.nl/">MUDE</a> TU Delft. This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">CC BY 4.0 License</a>.
diff --git a/synced_files/test.py b/synced_files/test.py
deleted file mode 100644
index 5a16119f85c6d1952fb4fadf06b9a1d9b1105bb3..0000000000000000000000000000000000000000
--- a/synced_files/test.py
+++ /dev/null
@@ -1,967 +0,0 @@
-# ---
-# jupyter:
-# jupytext:
-# text_representation:
-# extension: .py
-# format_name: percent
-# format_version: '1.3'
-# jupytext_version: 1.16.4
-# kernelspec:
-# display_name: Python 3 (ipykernel)
-# language: python
-# name: python3
-# ---
-
-# %% [markdown]
-# # GA 1.3: test
-#
-# <h1 style="position: absolute; display: flex; flex-grow: 0; flex-shrink: 0; flex-direction: row-reverse; top: 60px;right: 30px; margin: 0; border: 0">
-# <style>
-# .markdown {width:100%; position: relative}
-# article { position: relative }
-# </style>
-# <img src="https://gitlab.tudelft.nl/mude/public/-/raw/main/tu-logo/TU_P1_full-color.png" style="width:100px" />
-# <img src="https://gitlab.tudelft.nl/mude/public/-/raw/main/mude-logo/MUDE_Logo-small.png" style="width:100px" />
-# </h1>
-# <h2 style="height: 10px">
-# </h2>
-#
-# *[CEGM1000 MUDE](http://mude.citg.tudelft.nl/): Week 1.3. Due: Friday, September 20, 2024.*
-
-# %% [markdown]
-# <div style="background-color:#ffa6a6; color: black; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px; width: 95%"><p><b>Note:</b> don't forget to read the "Assignment Context" section of the README, it contains important information to understand this analysis.</p></div>
-
-# %%
-import numpy as np
-from scipy import interpolate
-from scipy.stats import norm
-import pandas as pd
-import matplotlib.pyplot as plt
-
-
-from functions import *
-
-np.set_printoptions(precision=3)
-
-# %% [markdown]
-# ## Part 0: Dictionary Review
-#
-# As described above, several functions in this assignment require the use of a Python dictionary to make it easier to keep track of important data, variables and results for the various _models_ we will be constructing and validating.
-#
-# _It may be useful to revisit PA 1.1, where there was a brief infroduction to dictionaires. That PA contains all the dictionary info you need for GA 1.3. A [read-only copy is here](https://mude.citg.tudelft.nl/2024/files/Week_1_1/PA_1_1_Catch_Them_All.html) and [the source code (notebook) is here](https://gitlab.tudelft.nl/mude/2024-week-1-1)._
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-#
-# $\textbf{Task 0.1}$
-#
-# Read and run the cell below to make sure you remember how to use a dictionary.
-#
-# Modify the function to print some of the other key-value pairs of the dictionary.
-#
-# <em>It may also be useful to use the cell below when working on later tasks in this assignment.</em>
-#
-# </p>
-# </div>
-
-# %%
-my_dictionary = {'key1': 'value1',
- 'key2': 'value2',
- 'name': 'Dictionary Example',
- 'a_list': [1, 2, 3],
- 'an_array': np.array([1, 2, 3]),
- 'a_string': 'hello'
- }
-
-def function_that_uses_my_dictionary(d):
- print(d['key1'])
-
- # SOLUTION:
- print(d['name'])
- print(d['a_list'])
- print(d['an_array'])
- print(d['a_string'])
-
- if 'new_key' in d:
- print('new_key exists and has value:', d['new_key'])
- return
-
-function_that_uses_my_dictionary(my_dictionary)
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-#
-# $\textbf{Task 0.2}$
-#
-# Test your knowledge by adding a new key <code>new_key</code> and then executing the function to print the value.
-#
-# </p>
-# </div>
-
-# %%
-# YOUR_CODE_HERE
-# function_that_uses_my_dictionary(my_dictionary)
-
-# SOLUTION:
-my_dictionary['new_key'] = 'new_value'
-function_that_uses_my_dictionary(my_dictionary)
-
-# %% [markdown] id="160d6250"
-# ## Task 1: Preparing the data
-#
-# Within this assignment you will work with two types of data: InSAR data and GNSS data. The cell below will load the data and visualize the observed displacements time. In this task we use the package `pandas`, which is really useful for handling time series. We will learn how to use it later in the quarter; for now, you only need to recognize that it imports the data as a `dataframe` object, which we then convert into a numpy array using the code below.
-
-# %% [markdown]
-# <div style="background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px"> <p>Tip: note that we have converted all observations to millimeters.</p></div>
-
-# %%
-gnss = pd.read_csv('./data/gnss_observations.csv')
-times_gnss = pd.to_datetime(gnss['times'])
-y_gnss = (gnss['observations[m]']).to_numpy()*1000
-
-insar = pd.read_csv('./data/insar_observations.csv')
-times_insar = pd.to_datetime(insar['times'])
-y_insar = (insar['observations[m]']).to_numpy()*1000
-
-gw = pd.read_csv('./data/groundwater_levels.csv')
-times_gw = pd.to_datetime(gw['times'])
-y_gw = (gw['observations[mm]']).to_numpy()
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 1.1:</b>
-#
-# Once you have used the cell above to import the data, investigate the data sets using the code cell below. Then provide some relevant summary information in the Markdown cell.
-#
-# <em>Hint: at the least, you should be able to tell how many data points are in each data set and get an understanding of the mean and standard deviation of each. Make sure you compare the different datasets and use consistent units.</em>
-#
-# </p>
-# </div>
-
-# %% [markdown] id="0491cc69"
-# <div style="background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px"> <p>The code below gives some examples of the quantitative and qualitative ways you could have looked at the data. It is more than you were expected to do; the important thing is that you showed the ability to learn something about the data and describe aspects that are relevant to our problem. We use a dictionary to easily access the different data series using their names, which are entered as the dictionary keys (also not expected of you, but it's hopefully fun to learn useful tricks).</div>
-
-# %%
-# YOUR_CODE_HERE
-
-# SOLUTION:
-data_list = ['y_gnss', 'y_insar', 'y_gw']
-data_dict = {data_list[0]: y_gnss,
- data_list[1]: y_insar,
- data_list[2]: y_gw}
-def print_summary(data):
- '''Summarize an array with simple print statements.'''
- print('Minimum = ', data.min())
- print('Maximum = ', data.max())
- print('Mean = ', data.mean())
- print('Std dev = ', data.std())
- print('Shape = ', data.shape)
- print('First value = ', data[0])
- print('Last value = ', data[-1])
- print('\n')
-
-for item in data_list:
- print('Summary for array: ', item)
- print('------------------------------------------------')
- print_summary(data_dict[item])
-
-# %%
-# SOLUTION:
-times_dict = {data_list[0]: times_gnss,
- data_list[1]: times_insar,
- data_list[2]: times_gw}
-def plot_data(times, data, label):
- plt.figure(figsize=(15,4))
- plt.plot(times, data, 'co', mec='black')
- plt.title(label)
- plt.xlabel('Times')
- plt.ylabel('Data [mm]')
- plt.show()
-
-plt.figure(figsize=(15,4))
-for i in range(3):
- plot_data(times_dict[data_list[i]],
- data_dict[data_list[i]],
- data_list[i])
-
-
-# %% [markdown]
-# <div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-#
-# $\textbf{Solution}$
-#
-# There are a lot more GNSS data points than InSAR or groundwater. The GNSS observations also have more noise, and what seem to be outliers. In this case the mean and standard deviation do not mean much, because there is clearly a trend with time. We can at least confirm that the time periods of measurements overlap, although the intervals between measurements is certainly not uniform (note that you don't need to do anything with the times, since they are pandas time series and we have not covered them yet).
-# </p>
-# </div>
-
-# %% [markdown]
-# You may have noticed that the groundwater data is available for different times than the GNSS and InSAR data. You will therefore have to *interpolate* the data to the same times for a further analysis. You can use the SciPy function ```interpolate.interp1d``` (read its [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html)).
-#
-# The cells below do the following:
-# 1. Define a function to convert the time unit
-# 2. Convert the time stamps for all data
-# 3. Use `interp1d` to interpolate the groundwater measurements at the time of the satellite measurements
-
-# %%
-def to_days_years(times):
- '''Convert the observation times to days and years.'''
-
- times_datetime = pd.to_datetime(times)
- time_diff = (times_datetime - times_datetime[0])
- days_diff = (time_diff / np.timedelta64(1,'D')).astype(int)
-
- days = days_diff.to_numpy()
- years = days/365
-
- return days, years
-
-
-# %%
-days_gnss, years_gnss = to_days_years(times_gnss)
-days_insar, years_insar = to_days_years(times_insar)
-days_gw, years_gw = to_days_years(times_gw)
-
-interp = interpolate.interp1d(days_gw, y_gw)
-
-GW_at_GNSS_times = interp(days_gnss)
-GW_at_InSAR_times = interp(days_insar)
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 1.2:</b>
-#
-# Answer/complete the code and Markdown cells below:
-# <ol>
-# <li>What is <code>interp</code>? (what kind of object is it, and how does it work?)</li>
-# <li>How did the groundwater observation array change? Be quantitative. </li>
-# </ol>
-#
-# </p>
-# </div>
-
-# %%
-# YOUR_CODE_HERE
-
-# SOLUTION:
-print('array size of GW_at_GNSS_times', len(GW_at_GNSS_times))
-print('array size of GW_at_InSAR_times', len(GW_at_InSAR_times))
-print('array size of GW before interpolation', len(y_gw))
-
-print('\nFirst values of times_gw:')
-print(times_gw[0:2])
-print('\nFirst values of y_gw:')
-print(y_gw[0:2])
-print('\nFirst values of times_gnss:')
-print(times_gnss[0:2])
-print('\nFirst values of GW_at_GNSS_times:')
-print(GW_at_GNSS_times[0:2])
-
-# %% [markdown]
-# **Write your answer in this Markdown cell.**
-
-# %% [markdown]
-# <div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-#
-# $\textbf{Solution}$
-# <ol>
-# <li><code>interp</code> is a function that will return a value (gw level) for the input(s) (date(s)). The interpolated value is found by linearly interpolating between the two nearest times in the gw observations.</li>
-# <li>The observation arrays of <code>GW_at_GNSS_times</code> and <code>GW_at_INSAR_times</code> changed in size to match the size of the GNSS and InSAR observations, respectively.</li>
-# </ol>
-# </p>
-# </div>
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 1.3:</b>
-#
-# Create a single plot to compare observed displacement for the GNSS and InSAR data sets.
-#
-# </p>
-# </div>
-
-# %%
-# plt.figure(figsize=(15,5))
-# plt.plot(YOUR_CODE_HERE, YOUR_CODE_HERE,
-# 'o', mec='black', label = 'GNSS')
-# plt.plot(YOUR_CODE_HERE, YOUR_CODE_HERE,
-# 'o', mec='black', label = 'InSAR')
-# plt.legend()
-# plt.ylabel('Displacement [mm]')
-# plt.xlabel('Time')
-# plt.show()
-
-# SOLUTION:
-plt.figure(figsize=(15,5))
-plt.plot(times_gnss, y_gnss, 'o', mec='black', label = 'GNSS')
-plt.plot(times_insar, y_insar, 'o', mec='black', label = 'InSAR')
-plt.legend()
-plt.ylabel('Displacement [mm]')
-plt.xlabel('Time')
-plt.show()
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 1.4:</b>
-# Describe the datasets based on the figure above and your observations from the previous tasks. What kind of deformation do you see? And what are the differences between both datasets? Be quantitative.
-# </p>
-# </div>
-
-# %% [markdown]
-# <div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-#
-# $\textbf{Solution}$
-#
-# The points obviously show subsidence, the displacement shows a similar pattern for both datasets. The GNSS data is much noisier than InSAR (range is around 60 mm versus only a few mm), but has a higher sampling rate. Also there seem to be more outliers in the GNSS data compared to InSAR, especially at the start of the observation period. InSAR has only observations every 6 days but is less noisy.
-# </p>
-# </div>
-
-# %% [markdown]
-# Before we move on, it is time to do a little bit of housekeeping.
-#
-# Have you found it confusing to keep track of two sets of variables---one for each data type? Let's use a dictionary to store relevant information about each model. We will use this in the plotting functions for this task (and again next week), so make sure you take the time to see what is happening. Review also Part 0 at the top of this notebook if you need a refresher on dictionaries.
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 1.5:</b>
-#
-# Run the cell below to define a dictionary for storing information about the two (future) models.
-#
-# </p>
-# </div>
-#
-
-# %%
-model_insar = {'data_type': 'InSAR',
- 'y':y_insar,
- 'times':times_insar,
- 'groundwater': GW_at_InSAR_times
- }
-
-model_gnss = {'data_type': 'GNSS',
- 'y':y_gnss,
- 'times':times_gnss,
- 'groundwater': GW_at_GNSS_times
- }
-
-# %% [markdown] id="76c9115b"
-# ## Task 2: Set-up linear functional model
-#
-# We want to investigate how we could model the observed displacements of the road. Because the road is built in the Green Heart we expect that the observed displacements are related to the groundwater level. Furthermore, we assume that the displacements can be modeled using a constant velocity. The model is defined as
-# $$
-# d = d_0 + vt + k \ \textrm{GW},
-# $$
-# where $d$ is the displacement, $t$ is time and $\textrm{GW}$ is the groundwater level (that we assume to be deterministic).
-#
-# Therefore, the model has 3 unknowns:
-# 1. $d_0$, as the initial displacement at $t_0$;
-# 2. $v$, as the displacement velocity;
-# 3. $k$, as the 'groundwater factor', which can be seen as the response of the soil to changes in the groundwater level.
-#
-#
-# As a group you will construct the **functional model** that is defined as
-# $$
-# \mathbb{E}(Y) = \mathrm{A x}.
-# $$
-#
-#
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 2.1:</b>
-#
-# Construct the design matrix $A$ (for both InSAR and GNSS observations), then show the first 5 observations and confirm the dimensions of $A$.
-# </p>
-# </div>
-
-# %%
-# YOUR_CODE_HERE
-
-# SOLUTION:
-A_insar = np.ones((len(times_insar), 3))
-A_insar[:,1] = days_insar
-A_insar[:,2] = GW_at_InSAR_times
-
-print ('The first 5 rows of the A matrix (InSAR) are:')
-print (A_insar[0:5, :])
-
-print ('The first 5 observations [mm] of y_insar are:')
-print (y_insar[0:5])
-
-m_insar = np.shape(A_insar)[0]
-n_insar = np.shape(A_insar)[1]
-print(f'm = {m_insar} and n = {n_insar}')
-
-# %%
-# YOUR_CODE_HERE
-
-# SOLUTION:
-A_gnss = np.ones((len(times_gnss), 3))
-A_gnss[:,1] = days_gnss
-A_gnss[:,2] = GW_at_GNSS_times
-
-print ('The first 5 rows of the A matrix (GNSS) are:')
-print (A_gnss[0:5, :])
-
-print ('\nThe first 5 observations [mm] of y_gnss are:')
-print (y_gnss[0:5])
-
-m_gnss = np.shape(A_gnss)[0]
-n_gnss = np.shape(A_gnss)[1]
-print(f'm = {m_gnss} and n = {n_gnss}')
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-#
-# $\textbf{Task 2.2}$
-#
-# Answer the following questions:
-#
-# - What is the dimension of the observables' vector $Y$?
-# - What are the unknowns of the functional model?
-# - What is the redundancy for this model?
-#
-# </p>
-# </div>
-#
-
-# %% [markdown]
-# <div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-#
-# $\textbf{Solution}$
-#
-# For InSAR:
-# <ol>
-# <li>The number of observations is 61.</li>
-# <li>The number of unknowns is 3.</li>
-# <li>The redundancy is 58.</li>
-# </ol>
-#
-# For GNSS:
-# <ol>
-# <li>The number of observations is 730.</li>
-# <li>The number of unknowns is 3.</li>
-# <li>The redundancy is 727.</li>
-# </ol>
-# </p>
-# </div>
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 2.3:</b>
-#
-# Add the A matrix to the dictionaries for each model. This will be used to plot results later in the notebook.
-#
-# </p>
-# </div>
-#
-
-# %%
-# model_insar['A'] = YOUR_CODE_HERE
-# model_gnss['A'] = YOUR_CODE_HERE
-
-# SOLUTION:
-model_insar['A'] = A_insar
-model_gnss['A'] = A_gnss
-
-print("Keys and Values (type) for model_insar:")
-for key, value in model_insar.items():
- print(f"{key:16s} --> {type(value)}")
-print("\nKeys and Values (type) for model_gnss:")
-for key, value in model_gnss.items():
- print(f"{key:16s} --> {type(value)}")
-
-# %% [markdown] id="9325d32b"
-# ## 3. Set-up stochastic model
-#
-# We will use the Best Linear Unbiased Estimator (BLUE) to solve for the unknown parameters. Therefore we also need a stochastic model, which is defined as
-# $$
-# \mathbb{D}(Y) = \Sigma_{Y}.
-# $$
-# where $\Sigma_{Y}$ is the covariance matrix of the observables' vector.
-#
-#
-#
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 3.1:</b>
-#
-# Construct the covariance matrix for each type of data and assume that
-#
-# - the observables are independent
-#
-# - the observables are normally distributed
-#
-# - the observables' standard deviation is
-#
-# - $\sigma_\textrm{InSAR} = 2$ mm
-# - $\sigma_\textrm{GNSS} = 15$ mm
-#
-# </p>
-# </div>
-#
-
-# %% colab={"base_uri": "https://localhost:8080/"} executionInfo={"elapsed": 40, "status": "ok", "timestamp": 1664699881875, "user": {"displayName": "C Yin", "userId": "14075875094781565898"}, "user_tz": -120} id="163acdb3" outputId="8bc99da3-a61e-4a25-8a54-c90f33299a57"
-# YOUR_CODE_HERE
-
-# SOLUTION:
-std_insar = 2 #mm
-
-Sigma_Y_insar = np.identity(len(times_insar))*std_insar**2
-
-print ('Sigma_Y (InSAR) is defined as:')
-print (Sigma_Y_insar)
-
-# %%
-# YOUR_CODE_HERE
-
-# SOLUTION:
-std_gnss = 15 #mm (corrected from original value of 5 mm)
-
-Sigma_Y_gnss = np.identity(len(times_gnss))*std_gnss**2
-
-print ('\nSigma_Y (GNSS) is defined as:')
-print (Sigma_Y_gnss)
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-#
-# $\textbf{Task 3.2}$
-#
-# Answer the following questions:
-#
-# - What information is contained in the covariance matrix?
-# - How do you implement the assumption that all observations are independent?
-# - What is the dimension of $\Sigma_{Y}$?
-# - How do you create $\Sigma_{Y}$?
-#
-# </p>
-# </div>
-
-# %% [markdown]
-# _Write your answer in this cell._
-
-# %% [markdown]
-# <div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-#
-# $\textbf{Solution}$
-#
-# - The covariance matrix contains information on the quality of the observations, where an entry on the diagonal represents the variance of one observation at a particular epoch. If there is an indication that for instance the quality for a particular time interval differs, different $\sigma$ values can be put in the stochastic model for these epochs.
-# - The off-diagonal terms in the matrix are related to the correlation between observations at different epochs, where a zero value on the off-diagonal indicates zero correlation.
-# - The dimension of the matrix is 61x61 for InSAR and 730x730 for GNSS.
-# - See code.
-# </p>
-# </div>
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 3.3:</b>
-#
-# Add <code>Sigma_Y</code> to the dictionaries for each model.
-#
-# </p>
-# </div>
-#
-
-# %%
-# model_insar['Sigma_Y] = YOUR_CODE_HERE
-# model_gnss['Sigma_Y'] = YOUR_CODE_HERE
-
-# SOLUTION:
-model_insar['Sigma_Y'] = Sigma_Y_insar
-model_gnss['Sigma_Y'] = Sigma_Y_gnss
-
-
-# %% [markdown] id="09e965bf"
-# ## 4. Apply best linear unbiased estimation
-#
-#
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 4.1:</b>
-#
-# Write a function to apply BLUE in the cell below and use the function to estimate the unknowns for the model using the data.
-#
-# Compute the modeled displacements ($\hat{\mathrm{y}}$), and corresponding residuals ($\hat{\mathrm{\epsilon}}$), as well as associated values (as requested by the blank code lines).
-# </p>
-# </div>
-#
-#
-
-# %% [markdown]
-# <div style="background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px"> <p><strong>Note on code implementation</strong>: you'll see that the functions in this assignment use a dictionary; this greatly reduces the number of input/output variables needed in a function. However, it can make the code inside the function more difficult to read due to the key syntax (e.g., <code>dict['variable_1']</code> versus <code>variable
-# _1</code>). To make this assignment easier for you to implement we have split these functions into three parts: 1) define variables from the dictionary, 2) perform analysis, 3) add results to the dictionary. Note that this is not the most efficient way to write this code; it is done here specifically for clarity and to help you focus on writing the equations properly and understanding the meaning of each term.</p></div>
-
-# %%
-def BLUE(d):
- """Calculate the Best Linear Unbiased Estimator
-
- Uses dict as input/output:
- - inputs defined from existing values in dict
- - outputs defined as new values in dict
- """
-
- y = d['y']
- A = d['A']
- Sigma_Y = d['Sigma_Y']
-
- # Sigma_X_hat = YOUR_CODE_HERE
- # x_hat = YOUR_CODE_HERE
-
- # y_hat = YOUR_CODE_HERE
-
- # e_hat = YOUR_CODE_HERE
-
- # Sigma_Y_hat = YOUR_CODE_HERE
- # std_y = YOUR_CODE_HERE
-
- # Sigma_e_hat = YOUR_CODE_HERE
- # std_e_hat = YOUR_CODE_HERE
-
- # SOLUTION:
- Sigma_X_hat = np.linalg.inv(A.T @ np.linalg.inv(Sigma_Y) @ A)
- x_hat = Sigma_X_hat @ A.T @ np.linalg.inv(Sigma_Y) @ y
-
- y_hat = A @ x_hat
-
- e_hat = y - y_hat
-
- Sigma_Y_hat = A @ Sigma_X_hat @ A.T
- std_y = np.sqrt(Sigma_Y_hat.diagonal())
-
- Sigma_e_hat = Sigma_Y - Sigma_Y_hat
- std_e_hat = np.sqrt(Sigma_e_hat.diagonal())
-
- d['Sigma_X_hat'] = Sigma_X_hat
- d['x_hat'] = x_hat
- d['y_hat'] = y_hat
- d['e_hat'] = e_hat
- d['Sigma_Y_hat'] = Sigma_Y_hat
- d['std_y'] = std_y
- d['Sigma_e_hat'] = Sigma_e_hat
- d['std_e_hat'] = std_e_hat
-
- return d
-
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 4.2:</b>
-#
-# Now that you have completed the function, apply it to our two models and then print values for the estimated parameters.
-# </p>
-# </div>
-#
-#
-
-# %%
-model_insar = BLUE(model_insar)
-x_hat_insar = model_insar['x_hat']
-
-# YOUR_CODE_HERE
-
-# SOLUTION:
-print ('The InSAR-estimated offset is', np.round(x_hat_insar[0],3), 'mm')
-print ('The InSAR-estimated velocity is', np.round(x_hat_insar[1],4), 'mm/day')
-print ('The InSAR-estimated velocity is', np.round(x_hat_insar[1]*365,4), 'mm/year')
-print ('The InSAR-estimated GW factor is', np.round(x_hat_insar[2],3), '[-]\n')
-
-model_gnss = BLUE(model_gnss)
-x_hat_gnss = model_gnss['x_hat']
-
-# YOUR_CODE_HERE
-
-# SOLUTION:
-print ('The GNSS-estimated offset is', np.round(x_hat_gnss[0],3), 'mm')
-print ('The GNSS-estimated velocity is', np.round(x_hat_gnss[1],4), 'mm/day')
-print ('The GNSS-estimated velocity is', np.round(x_hat_gnss[1]*365,4), 'mm/year')
-print ('The GNSS-estimated GW factor is', np.round(x_hat_gnss[2],3), '[-]')
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 4.3:</b>
-# Do the values that you just estimated make sense? Explain, using quantitative results.
-#
-# <em>Hint: all you need to do is use the figures created above to verify that the parameter values printed above are reasonable (e.g., order of magnitude, units, etc).</em>
-# </p>
-# </div>
-#
-#
-
-# %% [markdown]
-# <div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-#
-# $\textbf{Solution}$
-#
-# As long as the velocity is negative and around -0.02 mm/day or -10 mm/yr it makes sense if you compare with what you see in the plots with observations. Since load is applied on soil layers we expect the road to subside. We also expect to see a positive value for the GW factor.
-#
-# </p>
-# </div>
-
-# %% [markdown] id="65e42a43"
-# ## 5. Evaluate the precision
-#
-#
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 5:</b>
-#
-# What is the precision of the final estimates?
-#
-# Print the full covariance matrix of your estimates, and give an interpretation of the numbers in the covariance matrix.
-# </p>
-# </div>
-#
-#
-
-# %% colab={"base_uri": "https://localhost:8080/"} executionInfo={"elapsed": 16, "status": "ok", "timestamp": 1664699882186, "user": {"displayName": "C Yin", "userId": "14075875094781565898"}, "user_tz": -120} id="835eefc8" outputId="ad47d53d-c147-4b25-fb03-b967a3bb6f96"
-Sigma_X_hat_insar = model_insar['Sigma_X_hat']
-
-# YOUR_CODE_HERE
-
-# SOLUTION:
-print ('Covariance matrix of estimated parameters (InSAR):')
-print (Sigma_X_hat_insar)
-print ('\nThe standard deviation for the InSAR-estimated offset is',
- np.round(np.sqrt(Sigma_X_hat_insar[0,0]),3), 'mm')
-print ('The standard deviation for the InSAR-estimated velocity is',
- np.round(np.sqrt(Sigma_X_hat_insar[1,1]),4), 'mm/day')
-print ('The standard deviation for the InSAR-estimated GW factor is',
- np.round(np.sqrt(Sigma_X_hat_insar[2,2]),3), '[-]\n')
-
-Sigma_X_hat_gnss = model_gnss['Sigma_X_hat']
-
-# YOUR_CODE_HERE
-
-# SOLUTION:
-print ('Covariance matrix of estimated parameters (GNSS):')
-print (Sigma_X_hat_gnss)
-print ('\nThe standard deviation for the GNSS-estimated offset is',
- np.round(np.sqrt(Sigma_X_hat_gnss[0,0]),3), 'mm')
-print ('The standard deviation for the GNSS-estimated velocity is',
- np.round(np.sqrt(Sigma_X_hat_gnss[1,1]),4), 'mm/day')
-print ('The standard deviation for the GNSS-estimated GW factor is',
- np.round(np.sqrt(Sigma_X_hat_gnss[2,2]),3), '[-]')
-
-
-# %% [markdown]
-# <div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-#
-# $\textbf{Solution}$
-#
-# As shown above, the standard deviations of the estimated parameters are equal to the square root of the diagonal elements. Compared with the estimated values, the standard deviations seem quite small, except for the estimated offsets. Meaning that the complete estimated model can be shifted up or down.
-#
-# The off-diagonal elements show the covariances between the estimated parameters, which are non-zeros since the estimates are all computed as function of the same vector of observations and the same model. A different value for the estimated velocity would imply a different value for the GW factor and offset.
-#
-# </p>
-# </div>
-
-# %% [markdown] id="886efe26"
-# ## 6. Present and reflect on estimation results
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 6.1:</b>
-#
-# Complete the function below to help us compute the confidence intervals, then apply the function. Use a confidence interval of 96% in your analysis.
-#
-# <em>Hint: it can be used in exactly the same way as the <code>BLUE</code> function above, although it has one extra input.</em>
-# </p>
-# </div>
-#
-
-# %%
-def get_CI(d, alpha):
- """Compute the confidence intervals.
-
- Uses dict as input/output:
- - inputs defined from existing values in dict
- - outputs defined as new values in dict
- """
-
- std_e_hat = d['std_e_hat']
- std_y = d['std_y']
-
- # k = YOUR_CODE_HERE
- # CI_y = YOUR_CODE_HERE
- # CI_res = YOUR_CODE_HERE
-
- # SOLUTION:
- k = norm.ppf(1 - 0.5*alpha)
- CI_y = k*std_y
- CI_res = k*std_e_hat
- CI_y_hat = k*np.sqrt(d['Sigma_Y_hat'].diagonal())
-
- d['alpha'] = alpha
- d['CI_y'] = CI_y
- d['CI_res'] = CI_res
- d['CI_Y_hat'] = CI_y_hat
-
- return d
-
-
-# %%
-# model_insar = YOUR_CODE_HERE
-# model_gnss = YOUR_CODE_HERE
-
-# SOLUTION:
-model_insar = get_CI(model_insar, 0.04)
-model_gnss = get_CI(model_gnss, 0.04)
-
-# %% [markdown]
-# At this point we have all the important results entered in our dictionary and we will be able to use the plots that have been written for you in the next Tasks. In case you would like to easily see all of the key-value pairs that have been added to the dictionary, you can run the cell below:
-
-# %%
-print("Keys and Values (type) for model_insar:")
-for key, value in model_insar.items():
- print(f"{key:16s} --> {type(value)}")
-print("\nKeys and Values (type) for model_gnss:")
-for key, value in model_gnss.items():
- print(f"{key:16s} --> {type(value)}")
-
-# %% [markdown]
-# <div style="background-color:#AABAB2; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-# <b>Task 6.2:</b>
-#
-# Read the contents of file <code>functions.py</code> and identify what it is doing: you should be able to recognize that they use our model dictionary as an input and create three different figures. Note also that the function to create the figures have already been imported at the top of this notebook.
-#
-# Use the functions provided to visualize the results of our two models.
-# </p>
-# </div>
-#
-
-# %% [markdown] id="0491cc69"
-# <div style="background-color:#facb8e; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px"> <p><strong>Note</strong>: remember that you will have to use the same function to look at <em>both</em> models when writing your interpretation in the Report.</p></div>
-
-# %%
-# _, _ = plot_model(YOUR_CODE_HERE)
-
-# SOLUTION:
-_, _ = plot_model(model_insar)
-_, _ = plot_model(model_gnss)
-
-# %%
-# _, _ = plot_residual(YOUR_CODE_HERE)
-
-# SOLUTION:
-_, _ = plot_residual(model_insar)
-_, _ = plot_residual(model_gnss)
-
-# %%
-# _, _ = plot_residual_histogram(YOUR_CODE_HERE)
-
-# SOLUTION:
-_, _ = plot_residual_histogram(model_insar)
-_, _ = plot_residual_histogram(model_gnss)
-
-# %% [markdown]
-# <div style="background-color:#FAE99E; color: black; width:95%; vertical-align: middle; padding:15px; margin: 10px; border-radius: 10px">
-# <p>
-#
-# $\textbf{Solution: the True Model}$
-#
-# The data used in this exercise was generated using Monte Carlo Simulation. It is added to the plots here to illustrate where and how our models differ (it is your job to interpret "why").
-# </p>
-# </div>
-
-# %%
-k_true = 0.15
-R_true = -22
-a_true = 180
-d0_true = 10
-
-disp_insar = (d0_true + R_true*(1 - np.exp(-days_insar/a_true)) +
- k_true*GW_at_InSAR_times)
-disp_gnss = (d0_true + R_true*(1 - np.exp(-days_gnss/a_true)) +
- k_true*GW_at_GNSS_times)
-
-plot_model(model_insar, alt_model=('True model', times_insar, disp_insar));
-plot_model(model_gnss, alt_model=('True model', times_gnss, disp_gnss));
-
-# %%
-import ipywidgets as widgets
-from ipywidgets import interact
-
-# Function to update the plot based on slider values
-def update_plot(x0, x1, x2):
- plt.figure(figsize=(15,5))
- for m in [model_gnss]: #[model_insar, model_gnss]:
- plt.plot(m['times'], m['y'], 'o', label=m['data_type'])
- plt.ylabel('Displacement [mm]')
- plt.xlabel('Time')
-
- y_fit = model_gnss['A'] @ [x0, x1, x2]
- if (x0 == 0) & (x1 == 0) & (x2 == 1):
- plt.plot(model_gnss['times'], y_fit, 'r', label='Groundwater data', linewidth=2)
- else:
- plt.plot(model_gnss['times'], y_fit, 'r', label='Fit (GNSS)', linewidth=2)
-
- W = np.linalg.inv(model_gnss['Sigma_Y'])
- ss_res = (model_gnss['y'] - y_fit).T @ W @ (model_gnss['y'] - y_fit)
- plt.title(f'Mean of squared residuals: {ss_res:.0f}')
- plt.grid()
- plt.legend()
- plt.show()
-
-# Create sliders for x0, x1, and x2
-x0_slider = widgets.FloatSlider(value=0, min=-10, max=10, step=0.1, description='x0')
-x1_slider = widgets.FloatSlider(value=0, min=-0.1, max=0.1, step=0.001, description='x1')
-x2_slider = widgets.FloatSlider(value=0, min=-1, max=1, step=0.01, description='x2')
-
-# Use interact to create the interactive plot
-interact(update_plot, x0=x0_slider, x1=x1_slider, x2=x2_slider)
-
-
-# %%
-xhat_slider_plot(model_gnss['A'], model_gnss['y'], model_gnss['times'], model_gnss['Sigma_Y'])
-
-# %% [markdown]
-# **End of notebook.**
-# <h2 style="height: 60px">
-# </h2>
-# <h3 style="position: absolute; display: flex; flex-grow: 0; flex-shrink: 0; flex-direction: row-reverse; bottom: 60px; right: 50px; margin: 0; border: 0">
-# <style>
-# .markdown {width:100%; position: relative}
-# article { position: relative }
-# </style>
-# <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">
-# <img alt="Creative Commons License" style="border-width:; width:88px; height:auto; padding-top:10px" src="https://i.creativecommons.org/l/by/4.0/88x31.png" />
-# </a>
-# <a rel="TU Delft" href="https://www.tudelft.nl/en/ceg">
-# <img alt="TU Delft" style="border-width:0; width:100px; height:auto; padding-bottom:0px" src="https://gitlab.tudelft.nl/mude/public/-/raw/main/tu-logo/TU_P1_full-color.png" />
-# </a>
-# <a rel="MUDE" href="http://mude.citg.tudelft.nl/">
-# <img alt="MUDE" style="border-width:0; width:100px; height:auto; padding-bottom:0px" src="https://gitlab.tudelft.nl/mude/public/-/raw/main/mude-logo/MUDE_Logo-small.png" />
-# </a>
-#
-# </h3>
-# <span style="font-size: 75%">
-# © Copyright 2024 <a rel="MUDE" href="http://mude.citg.tudelft.nl/">MUDE</a> TU Delft. This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">CC BY 4.0 License</a>.