notebook_unit1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Sustainable Software Development, block course, March 2021**  \n",
    "*Scientific Software Center, Institute for Scientific Computing, Dr. Inga Ulusoy*\n",
    "\n",
    "# Analysis of the data\n",
    "\n",
    "Imagine you perform a \"measurement\" of some type and obtain \"scientific data\". You know what your data represents, but you have only a vague idea how different features in the data are connected, and what information you can extract from the data.\n",
    "\n",
    "You would start first with going through the data, making sure your data set is complete and that the result is reasonable. Imagine this already happened.\n",
    "\n",
    "In the next step, you would inspect your data more closely and try to identify structures. That is the step that we are focusing on in this unit.\n",
    "\n",
    "In the `data` folder, you will find several data files (`*.t` and `*.dat`). These are data files generated through some \"new approach\" that hasn't been used in your lab before. No previous analysis software exists, and you are going to establish a protocol for this \"new approach\" and \"publish your results\".\n",
    "\n",
    "The data can be grouped into two categories: \n",
    "1. data to be analyzed using statistical methods;\n",
    "2. data to be analyzed using numerical methods.\n",
    "\n",
    "In your hypothetical lab, you are an \"expert\" in one particular \"method\", and your co-worker is an \"expert\" in the other. Combined these two methods will lead to much more impactful results than if only one of you analyzed the data. Now, the task in this course is to be solved collaboratively with your team member working on one of the analysis approaches, and you working on the other. You will both implement functionality into the same piece of \"software\", but do so collaboratively through git.\n",
    "\n",
    "As you do not know yet which analysis is most meaningful for your data, and how to implement it, you will start with a jupyter notebook. You and your team member will work on the same notebook that will be part of a github repository for your project. This is the task for today. Discuss with your team members who will work on the statistical and who on the numerical analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1\n",
    "\n",
    "Generate a github repository with the relevant files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2\n",
    "\n",
    "Clone the repository to your local machine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3\n",
    "\n",
    "Start working on task 1 for your analysis approach. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4\n",
    "\n",
    "Create your own branch of the repository and commit your changes to your branch; push to the remote repository."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5\n",
    "\n",
    "Open a `pull request` so your team member can review your implementation. Likewise, your team member will ask you to review theirs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 6\n",
    "\n",
    "Merge the changes in your branch into `main`. Resolve conflicts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 7\n",
    "\n",
    "Repeat working on task; committing and pushing to your previously generated branch or a new branch; open a pull request; merge with main; until you have finished all the tasks in your analysis approach. Delete obsolete branches."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Start of the analysis notebook\n",
    "\n",
    "**Author : Group 1**  \n",
    "*Date : 03/01/23*  \n",
    "*Affiliation : The entity under whose name you are working on this*  \n",
    "\n",
    "Place the required modules in the top, followed by required constants and global functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# required modules\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# constants and global functions\n",
    "path = \"../../data/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reading of the data files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistical analysis\n",
    "\n",
    "Find correlations in the data sets. Analyse the data statistically and plot your results.  \n",
    "\n",
    "Here we would want to do everything with pandas and leave the data in a dataframe. The files that are relevant to you are `expect.t`, `npop.t` and `table.dat`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 1: Read in expec.t and plot relevant data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read and plot expec.t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can discard the entries norm, \\<x>, and \\<y> as these are mostly constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# eliminate columns based on the variance: if the variance of the values\n",
    "# in a column is below a given threshold, that column is discarded"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 2: Create plots of the relevant data and save as .pdf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 3: Read in file `npop.t` and analyze correlations in the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in npop.t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# discard all columns with variance below a set threshold - we can consider them as constant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the remaining columns. Seaborn prefers \"long format\" (one column for all measurement values, one column to indicate the type) as input, whereas the cvs is in \"wide format\" (one column per measurement type)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot ideally with seaborn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantify the pairwise correlation in the data\n",
    "\n",
    "- negative correlation: y values decrease for increasing x - large values of one feature correspond to small values of the other feature\n",
    "- weak or no correlation: no trend observable, association between two features is hardly observable\n",
    "- positive correlation: y values increase for decreasing x - small values of one feature correspond to small values of the other feature\n",
    "\n",
    "Remember that correlation does not indicate causation - the reason that two features are associated can lie in their dependence on same factors.\n",
    "\n",
    "Correlate the value pairs using Pearson's $r$. Pearson's $r$ is a measure of the linear relationship between features:\n",
    "\n",
    "$r = \\frac{\\sum_i(x_i − \\bar{x})(y_i − \\bar{y})}{\\sqrt{\\sum_i(x_i − \\bar{x})^2 \\sum_i(y_i − \\bar{y})^2}}$\n",
    "\n",
    "Here, $\\bar{x}$ and $\\bar{y}$ indicate mean values. $i$ runs over the whole data set. For a positive correlation, $r$ is positive, and negative for a negative correlation, with minimum and maximum values of -1 and 1, indicating a perfectly linear relationship. Weakly or not correlated features are characterized by $r$-values close to 0.\n",
    "\n",
    "Other measures of correlation that can be used are Spearman's rank (value pairs follow monotonic function) or Kendall's $\\tau$ (measures ordinal association), but they do not apply here. You can also define measures yourself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print the correlation matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The diagonal values tell us that each value is perfectly correlated with itself. We are not interested in the diagonal values and also not in the correlation with time. We also need to get rid of redundant entries. Finally, we need to find the value pairs that exhibit the highest linear correlation. We still want to know if it is positive or negative correlation, so we cannot get rid of the sign."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get rid of time column, lower triangular and diagonal entries of the correlation matrix\n",
    "# sort the remaing values according to their absolute value, but keep the sign"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the entries in the left column are not repeated if they do not change from the row above (so the fourth feature pair is MO3 and MO6)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 4: Print the resulting data to a file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write to file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 5: Calculate the Euclidean distance (L2 norm) for the vectors in `table.dat`\n",
    "\n",
    "\n",
    "The Euclidean distance measures the distance between to objects that are not points:\n",
    "\n",
    "$d(p,q) = \\sqrt{\\left(p-q\\right)^2}$\n",
    "\n",
    "In this case, consider each of the columns in table.dat as a vector in Euclidean space, where column $r(x)$ and column $v(x)$ denote a pair of vectors that should be compared, as well as $r(y)$ and $v(y)$, and r(z) and v(z).\n",
    "\n",
    "(Background: These are dipole moment components in different gauges, the length and velocity gauge.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in table.dat - I suggest reading it as a numpy array\n",
    "# replace the NaNs by zero"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now calculate how different the vectors in column 2 are from column 3, column 4 from column 5, and column 6 from column 7."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the Euclidean distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the result and save to a .pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print the result to a file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical analysis\n",
    "\n",
    "Analyze the data using autocorrelation functions and discrete Fourier transforms. Plot your results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define some global functions\n",
    "def read_data(file: str):\n",
    "    '''Read i the data of task one'''\n",
    "    data = np.loadtxt(file, skiprows=1)\n",
    "    return data\n",
    "\n",
    "def plot_task1(data):\n",
    "    fig, ax = plt.subplots(data.shape[1]-1, 1, sharex=True)\n",
    "    y_labels = ['x', 'y', 'z']\n",
    "    for i in range(data.shape[1]-1):\n",
    "        ax[i].plot(data[:, 0], data[:, i+1])\n",
    "        ax[i].set_ylabel(y_labels[i])\n",
    "    ax[-1].set_xlabel('time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 1: Read in `efield.t` and Fourier-transform relevant columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read and plot efield.t\n",
    "path_to_datafile = './data/efield.t'\n",
    "data = read_data(path_to_datafile)\n",
    "\n",
    "plot_task1(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we are interested in column 2 since the others are constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# discard the columns with variance below threshold - these are considered constant\n",
    "# in a column is below a given threshold, that column is discarded\n",
    "def select_columns(df, threshold):\n",
    "    variances = df.var()\n",
    "    # filter out the columns with variance below the threshold\n",
    "    selected_columns = variances[variances >= threshold].index.tolist()\n",
    "    # select the columns from the original dataframe\n",
    "    df_selected = df[selected_columns]\n",
    "    return df_selected\n",
    "\n",
    "clean_data = select_columns(data, 0.000001)\n",
    "\n",
    "\n",
    "clean_data.plot(x=\"time\", y=\"y\")\n",
    "\n",
    "clean_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# discrete Fourier transform of the remaining column: You only need the real frequencies\n",
    "N = len(clean_data[\"time\"])\n",
    "dt = clean_data[\"time\"].iloc[-1] - clean_data[\"time\"].iloc[0]\n",
    "fs = dt / N  # sampling frequency: dt/N\n",
    "\n",
    "yfft = np.fft.fft(clean_data[\"y\"])\n",
    "\n",
    "freqs = np.fft.fftfreq(N, 1 / fs)\n",
    "freqs = freqs[: N // 2]  # get positive frequencies only"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 2: Generate a plot of your results to be saved as pdf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot your results\n",
    "plt.plot(freqs, np.abs(yfft[: N // 2]))\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.ylabel(\"Magnitude\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 3: Calculate the autocorrelation function from nstate_i.t\n",
    "The autocorrelation function measures how correlated subsequent vectors are with an initial vector; ie. \n",
    "\n",
    "$\\Psi_{corr} = \\langle \\Psi(t=0) | \\Psi(t) \\rangle = \\int_0^{tfin} \\Psi(0)^* \\Psi(t) dt$\n",
    "\n",
    "Since we are in a numerical representation, the integral can be replaced with a sum; and the given vectors are already normalized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in as numpy array\n",
    "np_data_nstate = pd.read_table(path + \"nstate_i.t\", sep=\"\\s+\").to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# store the time column (column 0) in a vector and drop from array\n",
    "t_vec = np_data_nstate[0]\n",
    "np_data_nstate_mod = np.delete(np_data_nstate, [0], axis=0)\n",
    "print(f\"old: {np_data_nstate.shape}\", f\"new: {np_data_nstate_mod.shape}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# correct the data representation: this is in fact a complex matrix\n",
    "# the real part of each matrix column is contained in numpy array column 0, 2, 4, 6, ...\n",
    "# the imaginary part of each matrix column is contained in numpy array column 1, 3, 5, 7, ...\n",
    "# convert the array that was read as dtype=float into a dtype=complex array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for the autocorrelation function, we want the overlap between the first vector at time 0 and all\n",
    "# subsequent vectors at later times - the sum of the product of initial and subsequent vectors for each time step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 4: Generate a plot of your results to be saved as pdf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the autocorrelation function - real, imaginary and absolute part"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 5: Discrete Fourier transform of the autocorrelation function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# discrete Fourier-transform the autocorrelation function - now we need all frequency components,\n",
    "# also the negative ones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 6: Generate a plot of your results to be saved as pdf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the power spectrum (abs**2)"
   ]
  }
 ],
 "metadata": {
  "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.9.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "eb3e99d3805e0234c6935c2c6ba0163562eea1a2f720dd05ea3caf9e337a59f5"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}