data_exploration/figure/klmfigure.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5510d122",
   "metadata": {},
   "outputs": [],
   "source": [
    "#dataframes\n",
    "import pandas as pd\n",
    "import h5py\n",
    "\n",
    "#speedup\n",
    "import multiprocessing as mp\n",
    "\n",
    "#suppress warnings\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "pd.TimeSeries = pd.Series \n",
    "\n",
    "#math\n",
    "import numpy as np\n",
    "from scipy.spatial.distance import cdist\n",
    "from scipy import signal\n",
    "import math\n",
    "import copy\n",
    "\n",
    "\n",
    "#plots\n",
    "import pylab as plt\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import MaxNLocator\n",
    "import matplotlib\n",
    "from figurefirst import FigureLayout,mpl_functions\n",
    "import figurefirst\n",
    "import matplotlib.image as mpimg\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable # for colorbar\n",
    "import seaborn as sns\n",
    "sns.set()\n",
    "sns.set_style(\"whitegrid\")\n",
    "import figurefirst as fifi\n",
    "from IPython.display import display,SVG\n",
    "\n",
    "#misc\n",
    "import cairosvg\n",
    "import itertools\n",
    "import time\n",
    "np.set_printoptions(suppress=True)\n",
    "\n",
    "FS=7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "360ce55d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def __kalman_forward_update__(xhat_fm, P_fm, y, u, A, B, C, R, Q):\n",
    "    \"\"\"\n",
    "    :param xhat_fm:\n",
    "    :param P_fm:\n",
    "    :param y:\n",
    "    :param u:\n",
    "    :param A:\n",
    "    :param B:\n",
    "    :param C:\n",
    "    :param R:\n",
    "    :param Q:\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    \n",
    "    \n",
    "    I = np.array(np.eye(A.shape[0]))\n",
    "    gammaW = np.array(np.eye(A.shape[0]))\n",
    "\n",
    "    K_f = P_fm@C.T@np.linalg.pinv(C@P_fm@C.T + R)\n",
    "\n",
    "    if y is not None:\n",
    "        new_info = (y - C@xhat_fm)\n",
    "        new_info[np.isnan(new_info)] = 0\n",
    "        xhat_fp = xhat_fm + K_f@new_info\n",
    "        \n",
    "        P_fp = (I - K_f@C)@P_fm\n",
    "        xhat_fm = A@xhat_fp + B@u\n",
    "        P_fm = A@P_fp@A.T + gammaW@Q@gammaW.T\n",
    "        \n",
    "    return xhat_fp, xhat_fm, P_fp, P_fm\n",
    "\n",
    "\n",
    "def __kalman_forward_filter__(xhat_fm, P_fm, y, u, A, B, C, R, Q):\n",
    "    \"\"\"\n",
    "    :param xhat_fm:\n",
    "    :param P_fm:\n",
    "    :param y:\n",
    "    :param u:\n",
    "    :param A:\n",
    "    :param B:\n",
    "    :param C:\n",
    "    :param R:\n",
    "    :param Q:\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    if u is None:\n",
    "        u = np.array(np.zeros([B.shape[1], y.shape[1]]))\n",
    "\n",
    "    xhat_fp = None\n",
    "    P_fp = []\n",
    "    P_fm = [P_fm]\n",
    "\n",
    "    for i in range(y.shape[1]):\n",
    "        _xhat_fp, _xhat_fm, _P_fp, _P_fm = __kalman_forward_update__(xhat_fm[:, [-1]], P_fm[-1], y[:, [i]], u[:, [i]],\n",
    "                                                                     A, B, C, R, Q)\n",
    "        \n",
    "        if xhat_fp is None:\n",
    "            xhat_fp = _xhat_fp\n",
    "        else:\n",
    "            xhat_fp = np.hstack((xhat_fp, _xhat_fp))\n",
    "        xhat_fm = np.hstack((xhat_fm, _xhat_fm))\n",
    "\n",
    "        P_fp.append(_P_fp)\n",
    "        P_fm.append(_P_fm)\n",
    "\n",
    "    return xhat_fp, xhat_fm, P_fp, P_fm\n",
    "\n",
    "\n",
    "def __kalman_backward_smooth__(xhat_fp, xhat_fm, P_fp, P_fm, A):\n",
    "    \"\"\"\n",
    "    :param xhat_fp:\n",
    "    :param xhat_fm:\n",
    "    :param P_fp:\n",
    "    :param P_fm:\n",
    "    :param A:\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    N = xhat_fp.shape[1]\n",
    "\n",
    "    xhat_smooth = copy.copy(xhat_fp)\n",
    "    P_smooth = copy.copy(P_fp)\n",
    "    for t in range(N-2, -1, -1):\n",
    "        L = P_fp[t]@A.T@np.linalg.pinv(P_fm[t])\n",
    "        xhat_smooth[:, [t]] = xhat_fp[:, [t]] + L@(xhat_smooth[:, [t+1]] - xhat_fm[:, [t+1]])\n",
    "        P_smooth[t] = P_fp[t] - L@(P_smooth[t+1] - P_fm[t+1])\n",
    "\n",
    "    return xhat_smooth, P_smooth\n",
    "\n",
    "\n",
    "def klm_cal(data1,data2):\n",
    "    wind_vector = np.vstack((data_2.U.values, data_2.V.values)).T\n",
    "    wind_magnitude = np.linalg.norm(wind_vector, axis=1)\n",
    "    unit_wind_vector = wind_vector / np.vstack((wind_magnitude, wind_magnitude)).T\n",
    "\n",
    "\n",
    "    velocity_vector = np.vstack((data_2.ground_speed_x.values, data_2.ground_speed_y.values)).T\n",
    "\n",
    "    # There is probably a more efficient syntax\n",
    "    velocity_ground_relative_to_wind = [np.dot(velocity_vector[i,:], unit_wind_vector[i,:]) for i in range(len(velocity_vector))]\n",
    "\n",
    "    data_for_kf_original = pd.DataFrame({'time': data_2.time, \n",
    "                                    'velocity_ground_relative_to_wind': velocity_ground_relative_to_wind,\n",
    "                                    'distance_true': data_2.distance_from_source,\n",
    "                                    'distance_odor_observed': np.nan*np.ones_like(velocity_ground_relative_to_wind)})\n",
    "    # kalman filter is slow, let's subsample the data\n",
    "    data_for_kf = copy.copy( data_for_kf_original[::20] )\n",
    "    # put distance from odor predictions in the right places\n",
    "    for ix_odor in data_1.index:\n",
    "        ix_kf = np.argmin( np.abs(data_1.mean_whiff_time[ix_odor] - data_for_kf.time))\n",
    "        data_for_kf.distance_odor_observed.iloc[ix_kf] = data_1.predicted_distance.loc[ix_odor]\n",
    "    #     print(data_for_kf.distance_odor_observed[ix_kf])\n",
    "\n",
    "    dt = np.median(np.diff(data_for_kf.time.values))\n",
    "    y = np.array([data_for_kf.distance_odor_observed.values, \n",
    "                  data_for_kf.velocity_ground_relative_to_wind.values])\n",
    "\n",
    "    ### Define dynamics: assume constant velocity\n",
    "\n",
    "    # constant velocity\n",
    "    A = np.array([[1, dt], [0, 1]])\n",
    "\n",
    "    # no explicit control inputs\n",
    "    B = np.array([[0], [0]])\n",
    "\n",
    "    # we have measurements of distance and velocity\n",
    "    C = np.array([[1, 0], [0, 1]])\n",
    "\n",
    "\n",
    "    ### Define covariances\n",
    "    # measurement covariances: the measurements of the observed distance are worse than the velocity measurements\n",
    "    R = np.array([[1e-1, 0], [0, 1e-3]])\n",
    "    # R = np.array([[1e-2, 0], [0, 1e-3]]) #Forest\n",
    "\n",
    "    # model covariances: the model for constant velocity is probably worse than the model for distance = integral of velocity\n",
    "    Q = np.array([[1e-3, 0], [0, 1e-1]])\n",
    "    # Q = np.array([[1e-2, 0], [0, 1]]) #Forest\n",
    "\n",
    "    ### Define initial conditions\n",
    "    x0 = np.array([[0], [0]])\n",
    "    P0 = np.array(100*np.eye(2))\n",
    "    u = None # control input -- technically this could be the acceleration, but not needed\n",
    "\n",
    "    ### Choose a subset of the measurements for testing\n",
    "    ix_start = 0\n",
    "    ix_end = 100000\n",
    "\n",
    "    y_subset = copy.copy( y[:, ix_start:ix_end] )\n",
    "    y_subset[0,0] = 5 # start with a guess for initial distance (can't be nan)\n",
    "\n",
    "    time = data_for_kf.time[ix_start:ix_end].values\n",
    "\n",
    "    ### Run the filter and smoother\n",
    "    xhat_fp, xhat_fm, P_fp, P_fm = __kalman_forward_filter__(x0, P0, y_subset, u, A, B, C, R, Q)\n",
    "    xhat_smooth, _ = __kalman_backward_smooth__(xhat_fp, xhat_fm, P_fp, P_fm, A)\n",
    "    return xhat_smooth, time, data_for_kf\n",
    "\n",
    "def klm_cal_forest(data1,data2):\n",
    "    wind_vector = np.vstack((data_2.U.values, data_2.V.values)).T\n",
    "    wind_magnitude = np.linalg.norm(wind_vector, axis=1)\n",
    "    unit_wind_vector = wind_vector / np.vstack((wind_magnitude, wind_magnitude)).T\n",
    "\n",
    "\n",
    "    velocity_vector = np.vstack((data_2.ground_speed_x.values, data_2.ground_speed_y.values)).T\n",
    "\n",
    "    # There is probably a more efficient syntax\n",
    "    velocity_ground_relative_to_wind = [np.dot(velocity_vector[i,:], unit_wind_vector[i,:]) for i in range(len(velocity_vector))]\n",
    "\n",
    "    data_for_kf_original = pd.DataFrame({'time': data_2.time, \n",
    "                                    'velocity_ground_relative_to_wind': velocity_ground_relative_to_wind,\n",
    "                                    'distance_true': data_2.distance_from_source,\n",
    "                                    'distance_odor_observed': np.nan*np.ones_like(velocity_ground_relative_to_wind)})\n",
    "    # kalman filter is slow, let's subsample the data\n",
    "    data_for_kf = copy.copy( data_for_kf_original[::20] )\n",
    "    # put distance from odor predictions in the right places\n",
    "    for ix_odor in data_1.index:\n",
    "        ix_kf = np.argmin( np.abs(data_1.mean_whiff_time[ix_odor] - data_for_kf.time))\n",
    "        data_for_kf.distance_odor_observed.iloc[ix_kf] = data_1.predicted_distance.loc[ix_odor]\n",
    "    #     print(data_for_kf.distance_odor_observed[ix_kf])\n",
    "\n",
    "    dt = np.median(np.diff(data_for_kf.time.values))\n",
    "    y = np.array([data_for_kf.distance_odor_observed.values, \n",
    "                  data_for_kf.velocity_ground_relative_to_wind.values])\n",
    "\n",
    "    ### Define dynamics: assume constant velocity\n",
    "\n",
    "    # constant velocity\n",
    "    A = np.array([[1, dt], [0, 1]])\n",
    "\n",
    "    # no explicit control inputs\n",
    "    B = np.array([[0], [0]])\n",
    "\n",
    "    # we have measurements of distance and velocity\n",
    "    C = np.array([[1, 0], [0, 1]])\n",
    "\n",
    "\n",
    "    ### Define covariances\n",
    "    # measurement covariances: the measurements of the observed distance are worse than the velocity measurements\n",
    "    R = np.array([[1e-2, 0], [0, 1e-3]]) #Forest\n",
    "\n",
    "    # model covariances: the model for constant velocity is probably worse than the model for distance = integral of velocity\n",
    "    Q = np.array([[1e-2, 0], [0, 1]]) #Forest\n",
    "\n",
    "    ### Define initial conditions\n",
    "    x0 = np.array([[0], [0]])\n",
    "    P0 = np.array(100*np.eye(2))\n",
    "    u = None # control input -- technically this could be the acceleration, but not needed\n",
    "\n",
    "    ### Choose a subset of the measurements for testing\n",
    "    ix_start = 0\n",
    "    ix_end = 100000\n",
    "\n",
    "    y_subset = copy.copy( y[:, ix_start:ix_end] )\n",
    "    y_subset[0,0] = 5 # start with a guess for initial distance (can't be nan)\n",
    "\n",
    "    time = data_for_kf.time[ix_start:ix_end].values\n",
    "\n",
    "    ### Run the filter and smoother\n",
    "    xhat_fp, xhat_fm, P_fp, P_fm = __kalman_forward_filter__(x0, P0, y_subset, u, A, B, C, R, Q)\n",
    "    xhat_smooth, _ = __kalman_backward_smooth__(xhat_fp, xhat_fm, P_fp, P_fm, A)\n",
    "    return xhat_smooth, time, data_for_kf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f6155ff0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATcAAADlCAYAAAA2o7CiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYaUlEQVR4nO3df0zU9/0H8CechyCRKPbAD9pO1+7YTWVbJZuKmAgIjh+5tlMxTJZIMDE2c3M/WrsQkGk66ZImK4n9o4lbmtrUMCPGgykza7NgskSZWUWIKD9E4YAU1tBSjzuO9/cPv3cVDr3Pfcp97sP7no+/8Pjc3cunH55+jrvP+xMjhBAgIpJMbKQHICIKB5YbEUmJ5UZEUmK5EZGUWG5EJCWWGxFJieVGuqmrq0NOTg7S09PR1dUV6XFIclFdbvxh01dubi7OnDmDVatWRXqUqBDt+3dUlxt/2PSVmZkJRVEiPUbUiPb9e1GkB4ikzMzMSI9AKkxPT2NiYgJmsxkxMTERm0MIAY/Hg8TERMTGGv+4INr376gut2+CP3D6mZiYMNTLKqvViqVLl0Z6jLCRZd9muWnEHzj9mM1mAI/+jnFxcRGbw+12o6uryz+PrGTZt1luGvEHTj++o4e4uDgsXrw4wtMgokczepBl32a5acQfuNCdOHECLS0t+Oyzz7B//34sW7YMTU1NkR6LZpFl347qcuMPm76qqqpQVVUV6TGiRrTv31FdbvxhI5lF+/4t59trRBT1WG5EJCWWGxFJieVGRFJiuRGRlFhuRCQllhsRSYnlRkRSYrkRkZRYbkQkJZYbEUmJ5UZEUmK5EZGUWG5EJCWWGxFJieVGRFJiuRGRlFhuRCSloOVWV1eHnJwcpKenP/FyX16vF7W1tcjLy8OOHTvQ0NAw74OGS29vL0pLS1FQUIDS0lL09fVFdJ5Q8y4qKsLHH3+s85TaGC1rILS8CwsLceTIEZw7d07nKbUxYt66EkFcu3ZNDA4Oiu3bt4vbt2/Puc358+dFRUWF8Hq9YnR0VGRnZ4v79+8He2hDKC8vF42NjUIIIRobG0V5ebmq+7lcLnH9+nXhcrnmdZ5Q8x4cHBSbNm0S3d3d8zpHOBgtayFCy/urr74S//znP8XWrVsXxP5txLxD8U3nCHrklpmZCUVRnrpNc3Mzdu/ejdjYWCQnJyMvLw+XLl2atwIOl9HRUXR0dKC4uBgAUFxcjI6ODoyNjUVsJi15b9y4ES0tLTpNqI0RswZCzzspKQk5OTmG37+Nmree5uXqV06nE2lpaf4/K4qCoaEh1fefnp7GxMQEzGazrtff7O/vR0pKCqampjA1NQUAsFgsuHfvHhITE59638nJSQCPLhwbDkIIuN1u//M8bmBgABaLBZOTk3C73XjmmWfgdDpVP3Yk8h4cHER6evqMrK1WKwYGBlRnPTExEba8fZl88cUXAd978OABli1bhi+//BIAkJqaqnr/jtS+bfS81fA9t9fr1XR/Q1zab2Ji4om/7winnp4euFwutLe3+2+bnJxEd3c3pqenVT1GuOZ2u924e/funOXmm/FxHo9H9WNHKu9f//rXM7L+zW9+A7fbPeO2p+nt7Q3XaPB4PLh3796c/+5utxv9/f3+K5/P9W/yJJHKGjB23qH44osvsGTJkpDvNy/lpigKBgcHkZGRASDwSC4Y305jtVoRFxc3HyOpoigK6urqYLPZYDKZ4PV6MT4+jqysLCQnJ/u3a29vx/r162fcd3p6Gg8fPsSiRYvC8j+y2WzGt771LbzwwgsB33vuuedgNpthtVrhdrvx2WefYd26dSE9NqB/3mrNlXe4xcXF4YUXXsB3vvOdgO+tWbMGCQkJsFqt6Orqwv/+9z8899xzqh7X6FkDkclbDbfbja6uLixdulTT/eel3Hbu3ImGhgbk5+fj888/x5UrV3DmzBnV9/eVQ1xcHBYvXjwfI6mSlpYGm82Gf/zjH7Db7bhw4QJsNtucv4OZa66EhISwzRYbG4vExMQ5/2GLi4vhcDhgt9vx5Zdfoq2tDUeOHFH92JHKOxR6zxUTE/PEPAoLC9HY2Ii8vDyMj4/j448/xocffqj6cQFjZw3on3coTCaTpvsFfUPhxIkT2LZtG4aGhrB//34UFRUBAA4cOICbN28CAOx2O1avXo38/Hzs2bMHr776Kp599llNA+nt2LFj+OCDD1BQUIAPPvgAtbW1EZ0n1Lz37duHl19+GatXr47k2AtWKHkXFxejuroaBw8eXDD7dzSLEUKISA8xOTnpPzQ24v8gbW1t2LhxY6THmJOW7Ji3NjJmDciV9+N4hgIRSYnlRkRSYrkRkZRYbkQkJZYbEUmJ5aaBEVZb8K1mkZGRgfv37+v+/NHMKJ/cjxZa82a5aVBTU4OysjJcvnwZZWVlqK6u1n2G3NxcnDlzJqQzQYiiCcstREZZbUHNahZE0YzlFiKn04nU1FT/KSEmkwkpKSkhrcpBROHHciMiKbHcQqQoCoaHh/1rTHm9XoyMjPAlIpHBsNxCtGLFCthsNjgcDgCAw+GAzWabsUQSEUUey00DI6wk4lvNYnh4GG+++SZefvll3WeIVr/73e/8q4dQ+GnN2xAr8S40zz//fMSv8FVVVYWqqqoZKyeQPv72t78ZdoUPGWnNm0duRCQllhsRSYnlRkRSYrkRkZRYbkQkJZYbEUmJ5UZEUmK5EZGUWG5EJCWWGxFJieVGRFJSdW5pb28vjh49is8//xzLli1DXV0d1qxZM2Ob+vp6fPjhh0hJSQEAvPjii6ipqZn3gaNBKHlbLBa4XC5s2rQJf/jDHyIz8ALGrOWlqtx81wyw2+24cOECqqur8f777wds99JLL+H111+f9yGjTSh5/+pXv+KJ898As5ZX0JelRrlmQLRg3vph1nILWm6hXDOgqakJJSUlqKiowI0bN+Z/2igQat4//elP8cc//hH//e9/9R51wWPWcpu39dz27t2LgwcPwmw24+rVqzh06BCam5uxfPly1Y/R3t4+X+PMu7a2Nl2ep6enBw8fPpzxfC6XC52dnXC5XP7bvve97+FPf/oTFi1ahJs3b+Lw4cP4+9//zrxDwKy/ptf+raeg5fb4NQNMJtMTrxlgsVj8X2dlZUFRFNy5cwc/+tGPVA+zfv16Qy4C2NbWho0bN+ryXGvWrEFdXR1+8IMf+PMeHx9HTk7OnEuZT05OAgBWrlzJvEPErB/Rc/8OhW8hVq2CvixVe82A4eFh/9ednZ0YGBjA2rVrNQ8WrbTk3dfXh8HBQeYdImYtN1UvS48dO4ajR4/i1KlTSEpKQl1dHQDgwIEDOHz4MDZs2IC3334bt27dQmxsLMxmM956660ZR3OkXih5x8TEYGpqCm+++Sbz1oBZyytGCCEiPcTj1wEw4qG7UQ/bAW3ZMW9tZMwakCvvx/EMBSKSEsuNiKTEciMiKbHciEhKLDcikhLLjYikxHIjIimx3IhISiw3IpISy42IpMRyIyIpsdyISEosNyKSEsuNiKTEciMiKbHciEhKLDcikhLLjYikxHIjIimx3IhISiw3IpISy42IpMRyIyIpsdyISEosNyKSkqpy6+3tRWlpKQoKClBaWoq+vr6AbbxeL2pra5GXl4cdO3agoaFhvmeNGqHkXVhYiCNHjuDcuXP6DyoBZi0vVeVWU1ODsrIyXL58GWVlZaiurg7Y5uLFi+jv70dLSwvOnj2L+vp6PHjwYN4Hjgah5O1wOFBbW4t3332XeWvArOUVtNxGR0fR0dGB4uJiAEBxcTE6OjowNjY2Y7vm5mbs3r0bsbGxSE5ORl5eHi5duhSeqSWmJe+kpCTk5OQw7xAxa7ktCraB0+lEamoqTCYTAMBkMiElJQVOpxPJyckztktLS/P/WVEUDA0NqRpCCAEAcLvdIQ2vp8nJSV2ep7+/HykpKZiamsLU1BQAwGKx4N69e0hMTPRvNzAwAIvF4s9s5cqVzDtEzPpreu3fofBl5sswVEHLTQ8ejwcA0NXVFeFJnqy9vV2X5+np6YHL5ZrxfJOTk+ju7sb09HTAbT5er1f1czDvR5j11/Tav7XweDyIj48P+X5By01RFAwPD8Pr9cJkMsHr9WJkZASKogRsNzg4iIyMDACBR3JPk5iYCKvVCrPZjJiYmJD/EjJRFAV1dXWw2Wz+vMfHx5GVlTXjSHnNmjVISEjAunXr4PF44HA4sGrVKlXPwbwfYdbGJoSAx+OZcRQd6gMEtW/fPtHY2CiEEKKxsVHs27cvYJtz586JiooK4fV6xejoqMjOzhb9/f1qHp5mYd76YdbyUlVud+/eFbt27RL5+fli165doru7WwghRGVlpfj000+FEEJMTU2J6upqkZubK3Jzc8VHH30Uvqklx7z1w6zlFSOExt/WEREZGM9QICIpsdyISEosNyKSEsuNiKTEciMiKelabkZeXUTNbPX19di8eTPsdjvsdjtqa2t1ma2urg45OTlIT09/4ifd58qNeWujJe93333XsFkDxs1b676tip6fOykvL5/xgcny8vKAbc6fPx/wgcn79+8bYrZ33nlHnDx5MuyzzHbt2jUxODgotm/fLm7fvj3nNnPltmfPHuatgZa8161bJ06fPi2EMF7WQhg3b637tprcdDtyM/LqImpni5TMzMyA091mm53b1q1b0dnZybw1CDVv8f8fFfWdR2qkrAFj561l31abm27l9rTVRWZvp3V1kXDPBgBNTU0oKSlBRUUFbty4Eda5QjE7t4SEBMTHxzPvMHk8N6fTiaSkJIyMjAAwVta+513IeWvNzRCrgiwUe/fuxcGDB2E2m3H16lUcOnQIzc3NWL58eaRHkxLz1pdseet25Pb46iIAgq4u4uN0OrFy5UpDzGaxWGA2mwEAWVlZUBQFd+7cCetsas3O7eHDh3C5XMw7TB7PTVEUjI+PIyUlBYCxsvY970LOW2tuupXbihUrYLPZ4HA4AAAOhwM2m23G0jIAsHPnTjQ0NGB6ehpjY2O4cuUKCgoKDDHb8PCw/+vOzk4MDAxg7dq1YZ1Nrdm5tba2Ij09nXmHyeO5+ZYy8hWDkbIGFn7emnOb5zc/nsrIKzCome21114TRUVFoqSkRLzyyivik08+0WW248ePi+zsbGGz2cSWLVtEYWFhwGxz5ca8tdGSd319vWGzFsK4eWvdt9XgqiBEJCWeoUBEUmK5EZGUWG5EJCWWGxFJieVGRFJiuRGRlFhuRCQllhsRSYnlRkRSYrkRkZRYbkQkpaDlFtY1zikA89YX85ZX0HLLzc3FmTNnsGrVqiduc/HiRfT396OlpQVnz55FfX09Hjx4MK+DRgvmrS/mLa+gK/FmZmYGfZAnrXFeWVmpaojp6WlMTEzAbDb718aKVhs2bAAALF++HNPT05icnAzY5l//+hdKS0vh8XiwZMkS7N69m3lrFErebrcbZrMZ+fn5qvNm1toJIeDxeJCYmIjY2NB/gzYvy4x/07XhJyYmnviSIFr9/ve/x8TEBNrb2wO+97Of/QwA/N/bsmUL/v3vf6t+bOYdSE3et27dAgB8+9vfRk9Pj6rHZdbfnNVqxdKlS0O+nyGuoeBbwdRqtSIuLi7C0wRqb2/H+vXrIz3GnNxuN7q6uuY84ngS5q2NL2vfVa7UMHrWgPHz9mUYqnkpN98a5xkZGQACj+SC8R2ux8XFYfHixfMx0rwz6lw+zzzzjOptmfc3MzIyovraBwsha8DYeWt9OT8vHwWJ1Nrw9LVt27ZFeoSo0drayv17AQhabidOnMC2bdswNDSE/fv3o6ioCABw4MAB3Lx5EwBgt9uxevVq5OfnY8+ePXj11Vfx7LPPhndymiHYhW1p/vz85z/n/r0ABH1ZWlVVhaqqqoDb33vvPf/XJpMJtbW18zsZkUH5rtxOxsYzFIhISiw3IpISy42IpMRyIyIpsdyISEosNyKSEsuNiKTEciMiKbHciEhKLDcikhLLjYikxHIjIimx3IhISiw3IpISy42IpMRyIyIpsdyISEosNyKSEsuNiKTEciMiKbHciEhKLDcikpKqcuvt7UVpaSkKCgpQWlqKvr6+gG3q6+uxefNm2O122O12XupPZ3/+858jPULUYNYLQ9DrlgJATU0NysrKYLfbceHCBVRXV+P9998P2O6ll17C66+/Pu9DUnC//OUvIz1C1GDWC0PQI7fR0VF0dHT4L0RbXFyMjo4OjI2NhX04IiKtgh65OZ1OpKamwmQyAXh0dfmUlBQ4nU4kJyfP2LapqQmtra2wWCz4xS9+gR/+8IchDdPe3h7S9npqa2uL9AhPdevWLbz44osh3Yd5ayNb1oCx89ZMBHHz5k1RWFg447af/OQnor29fcZtIyMjwu12CyGEaG1tFZs2bRJjY2PBHl4IIYTL5RLXr18XLpdL1fZ6u379eqRHeCJfdoWFhcw7zGTMWgjj5601u6AvSxVFwfDwMLxeLwDA6/ViZGQEiqLM2M5iscBsNgMAsrKyoCgK7ty5E4Y6prmkpKQwb50w64UhaLmtWLECNpsNDocDAOBwOGCz2QJekg4PD/u/7uzsxMDAANauXTvP49KTDA0NMW+dMOuFQdW7pceOHcPRo0dx6tQpJCUloa6uDgBw4MABHD58GBs2bMDbb7+NW7duITY2FmazGW+99RYsFktYh6evvfHGG8xbJ8x6YVBVbs8//zwaGhoCbn/vvff8X/sKjyLjxz/+caRHiBrMemHgGQpEJCWWGxFJieVGRFJiuRGRlFhuRCQllhsRSYnlRkRSYrkRkZRYbkQkJZYbEUmJ5UZEUmK5EZGUWG5EJCWWGxFJieVGRFJiuRGRlFhuRCQllhsRSYnlRkRSYrkRkZRYbkQkJZYbEUlJVbn19vaitLQUBQUFKC0tRV9fX8A2Xq8XtbW1yMvLw44dO+a8FCCFT3Nzc6RHiBrMemFQVW41NTUoKyvD5cuXUVZWhurq6oBtLl68iP7+frS0tODs2bOor6/HgwcP5n1gmttf//pX5q0TZr0wBL0o8+joKDo6OvCXv/wFAFBcXIzjx49jbGwMycnJ/u2am5uxe/duxMbGIjk5GXl5ebh06RIqKyuDDiGEAAC43W6tf4+wm5ycjPQIc/Jllp2dzbzDTNasAWPn7cswVEHLzel0IjU1FSaTCQBgMpmQkpICp9M5o9ycTifS0tL8f1YUBUNDQ6qG8Hg8AICurq6QhtdTe3t7pEd4qrS0NPT09Kjalnl/M7JlDRg7b4/Hg/j4+JDvF7Tc9JCYmAir1Qqz2YyYmJhIjxNRXV1dOHnyJE6fPu2/bf/+/XjjjTdgtVr9t1VWVuK3v/0t0tPT4fF48J///Ef1czDvR5i1sQkh4PF4kJiYqOn+QctNURQMDw/D6/XCZDLB6/ViZGQEiqIEbDc4OIiMjAwAgUdyTxMbG4ulS5dqGF8+aWlpuH37NhYtWuTPu6urC6tWrcLixYv92yUkJMDpdOL73/8+4uPjmbcGzNr4tByx+QR9Q2HFihWw2WxwOBwAAIfDAZvNNuMlKQDs3LkTDQ0NmJ6extjYGK5cuYKCggLNg0Ur5q0fZi05ocLdu3fFrl27RH5+vti1a5fo7u4WQghRWVkpPv30UyGEEFNTU6K6ulrk5uaK3Nxc8dFHH6l5aJoD89YPs5ZXjBAa34ogIjIwnqFARFJiuRGRlFhuRCQllhsRSUnXcjPyCfhqZquvr8fmzZtht9tht9tRW1ury2x1dXXIyclBenr6Ez/pPlduzFsbLXmfOnXKsFkDxs1b676tip5vzZaXl4vGxkYhhBCNjY2ivLw8YJvz58+LiooK4fV6xejoqMjOzhb37983xGzvvPOOOHnyZNhnme3atWticHBQbN++Xdy+fXvObebKbc+ePcxbAy15r1u3Tpw+fVoIYbyshTBu3lr3bTW56Xbk5jsBv7i4GMCjE/A7OjowNjY2Y7snnYBvhNkiJTMzM+CMkNlm57Z161Z0dnYybw1CzVv8/6epfOeRGilrwNh5a9m31eamW7k97QT82dtpPQE/3LMBQFNTE0pKSlBRUYEbN26Eda5QzM5tyZIliI+PZ95h8nhuTqcTSUlJGBkZAWCsrH3Pu5Dz1pqbIU6cXyj27t2LgwcPwmw24+rVqzh06BCam5uxfPnySI8mJeatL9ny1u3I7fET8AEEPQHfx+l0YuXKlYaYzWKxwGw2AwCysrKgKAru3LkT1tnUmp3bV199BZfLxbzD5PHcFEXB+Pg4UlJSABgra9/zLuS8teamW7kZ+SRltbMNDw/7v+7s7MTAwADWrl0b1tnUmp1ba2srvvvd7zLvMHk8N99SRr5iMFLWwMLPW3Nu8/zmx1MZ+SRlNbO99tproqioSJSUlIhXXnlFfPLJJ7rMdvz4cZGdnS1sNpvYsmWLKCwsDJhtrtyYtzZa8q6vrzds1kIYN2+t+7YaPHGeiKTEMxSISEosNyKSEsuNiKTEciMiKbHciEhKLDcikhLLjYikxHIjIin9H8pYrMvznxX3AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 595.276x841.89 with 9 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "layout = fifi.svg_to_axes.FigureLayout('../../../Figure/fig7.svg',\n",
    "                                           autogenlayers=True, make_mplfigures=True, hide_layers=[])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb4c7992",
   "metadata": {},
   "source": [
    "## HWS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "27cdf39e",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir='~/DataAnalysis/data/Sprints/HighRes/KF/'\n",
    "data_1=pd.read_hdf(dir+'AllpredictionFiltered.h5')\n",
    "data_2=pd.read_hdf(dir+'WindyMASigned.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "fedde384",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['mc_min', 'mc_max', 'mc_mean', 'mc_std_dev', 'mc_k', 'wf_min', 'wf_max',\n",
       "       'wf_mean', 'wf_std_dev', 'wf_k', 'wd_min', 'wd_max', 'wd_mean',\n",
       "       'wd_std_dev', 'wd_k', 'ma_min', 'ma_max', 'ma_mean', 'ma_std_dev',\n",
       "       'ma_k', 'st_min', 'st_max', 'st_mean', 'st_std_dev', 'st_k', 'distance',\n",
       "       'mean_whiff_time', 'nwhiffs', 'predicted_distance'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_1.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ea0ea0cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "xhat_smooth,time,data_for_kf = klm_cal(data_1,data_2)\n",
    "\n",
    "ax=layout.axes[('klm','hws')]\n",
    "ax.grid(False)\n",
    "line_width=0.5\n",
    "marker_size=0.5\n",
    "ax.plot(data_1.mean_whiff_time, data_1.predicted_distance, '.', \n",
    "        color='black', label='Predicted Distance',linewidth=line_width,\n",
    "        markersize=marker_size*2,\n",
    "        rasterized= True)\n",
    "ax.plot(time, xhat_smooth[0,:], color='black', label='Smoothed Distance', \n",
    "        linewidth=line_width,markersize=marker_size,\n",
    "        rasterized= True)\n",
    "ax.plot(data_for_kf.time, data_for_kf.distance_true, '--', color='red',\n",
    "        label='True Distance',linewidth=line_width, markersize=marker_size,rasterized= True)\n",
    "\n",
    "ax.set_rasterization_zorder(-20)\n",
    "ax.set_xlim(1000, 2000)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left','bottom'],\n",
    "                            spine_locations={'left':4,'bottom':4},\n",
    "                            smart_bounds=True,\n",
    "                            yticks=[-10,70],\n",
    "                            xticks=[1000,2000],\n",
    "                            linewidth=0.5)\n",
    "\n",
    "ax.set_xticklabels([])\n",
    "# ax.set_xticklabels(['1000','2000'], rotation=90)\n",
    "ax.axvspan(xmin=1000, xmax=1300, alpha=0.4, color='grey')\n",
    "ax.set_xlabel('1000 secs')\n",
    "ax.xaxis.set_label_coords(0.5, -0.1)\n",
    "\n",
    "ax.set_ylabel('Distance, m')\n",
    "ax.yaxis.set_label_coords(-0.08, 0.5)\n",
    "# ax.set_title('WS>3.5m/s', x=0.5, y=0.95)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9e488610",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('klm','hws1')]\n",
    "ax.grid(False)\n",
    "line_width=0.5\n",
    "marker_size=0.5\n",
    "ax.plot(data_1.mean_whiff_time, data_1.predicted_distance, '.', \n",
    "        color='black', label='Predicted Distance',linewidth=line_width,\n",
    "        markersize=marker_size*2,\n",
    "        rasterized= True)\n",
    "ax.plot(time, xhat_smooth[0,:], color='black', label='Smoothed Distance', \n",
    "        linewidth=line_width,markersize=marker_size,\n",
    "        rasterized= True)\n",
    "ax.plot(data_for_kf.time, data_for_kf.distance_true, '--', color='red',\n",
    "        label='True Distance',linewidth=line_width, markersize=marker_size,rasterized= True)\n",
    "\n",
    "ax.set_rasterization_zorder(-20)\n",
    "ax.set_xlim(1000,1300)\n",
    "ax.set_xlabel('300 secs')\n",
    "ax.xaxis.set_label_coords(0.5, -0.1)\n",
    "\n",
    "ax.set_ylabel('Distance, m')\n",
    "ax.yaxis.set_label_coords(-0.08, 0.5)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left','bottom'],\n",
    "                            spine_locations={'bottom':4, 'left':4},\n",
    "                            smart_bounds=True,\n",
    "                            yticks=[-10,60],\n",
    "                            xticks=[1000,1300],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([], rotation=90)\n",
    "\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ed4725ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('klm','hws2')]\n",
    "# f,ax=plt.subplots(1,1)\n",
    "ax.grid(False)\n",
    "error = xhat_smooth[0,:] - data_for_kf.distance_true.values\n",
    "\n",
    "h = ax.hist( np.abs(error), bins=100, color='gray',histtype='stepfilled', rasterized=True)\n",
    "\n",
    "median_error = np.median(np.abs(error))\n",
    "\n",
    "ax.axvline(median_error, 0, ymax = 0.5,color='blue', linewidth=line_width)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "\n",
    "ax.set_ylim(0,1000)\n",
    "ax.set_xlim(0,20)\n",
    "# ax.text(14,4000,'Median\\nError=2.92',color='black')\n",
    "mpl_functions.adjust_spines(ax,['bottom'],\n",
    "                            spine_locations={'left':4,'bottom':4},\n",
    "                            smart_bounds=True,\n",
    "                            yticks=[0,3000],\n",
    "                            xticks=[0,20],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "ax.set_xlabel('20m')\n",
    "ax.xaxis.set_label_coords(0.5, -0.15)\n",
    "\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "9a0e46e1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.901826958741591"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error = xhat_smooth[0,:] - data_for_kf.distance_true.values\n",
    "\n",
    "h = ax.hist( np.abs(error), bins=100, color='gray',histtype='stepfilled', rasterized=True)\n",
    "\n",
    "median_error = np.median(np.abs(error))\n",
    "median_error"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2d5cbee",
   "metadata": {},
   "source": [
    "### LWS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6f3b3946",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir='~/DataAnalysis/data/Sprints/HighRes/KF/'\n",
    "data_1=pd.read_hdf(dir+'AllpredictionFiltered.h5')\n",
    "data_2=pd.read_hdf(dir+'NotWindyMASigned.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "bbba3432",
   "metadata": {},
   "outputs": [],
   "source": [
    "xhat_smooth,time,data_for_kf = klm_cal(data_1,data_2)\n",
    "\n",
    "ax=layout.axes[('klm','lws')]\n",
    "ax.grid(False)\n",
    "line_width=0.5\n",
    "marker_size=0.5\n",
    "ax.plot(data_1.mean_whiff_time, data_1.predicted_distance, '.', \n",
    "        color='black', label='Predicted Distance',linewidth=line_width,\n",
    "        markersize=marker_size*2,\n",
    "        rasterized= True)\n",
    "ax.plot(time, xhat_smooth[0,:], color='black', label='Smoothed Distance', \n",
    "        linewidth=line_width,markersize=marker_size,\n",
    "        rasterized= True)\n",
    "ax.plot(data_for_kf.time, data_for_kf.distance_true, '--', color='red',\n",
    "        label='True Distance',linewidth=line_width, markersize=marker_size,rasterized= True)\n",
    "\n",
    "ax.set_rasterization_zorder(-20)\n",
    "ax.set_xlim(2000,3000)\n",
    "ax.set_ylim(-10,70)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['bottom'],\n",
    "                            spine_locations={'bottom':4},\n",
    "                            smart_bounds=True,\n",
    "#                             yticks=[-10,70],\n",
    "                            xticks=[2000,3000],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([], rotation=90)\n",
    "ax.axvspan(xmin=2200, xmax=2500,alpha=0.4, color='grey')\n",
    "# ax.set_xlabel('1000, s')\n",
    "# ax.xaxis.set_label_coords(0.5, -0.1)\n",
    "\n",
    "# ax.set_title('WS<3.5m/s', x=0.5, y=0.95)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7021d4af",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('klm','lws1')]\n",
    "ax.grid(False)\n",
    "line_width=0.5\n",
    "marker_size=0.5\n",
    "ax.plot(data_1.mean_whiff_time, data_1.predicted_distance, '.', \n",
    "        color='black', label='Predicted Distance',linewidth=line_width,\n",
    "        markersize=marker_size*2,\n",
    "        rasterized= True)\n",
    "ax.plot(time, xhat_smooth[0,:], color='black', label='Smoothed Distance', \n",
    "        linewidth=line_width,markersize=marker_size,\n",
    "        rasterized= True)\n",
    "ax.plot(data_for_kf.time, data_for_kf.distance_true, '--', color='red',\n",
    "        label='True Distance',linewidth=line_width, markersize=marker_size,rasterized= True)\n",
    "\n",
    "ax.set_rasterization_zorder(-20)\n",
    "ax.set_xlim(2200,2500)\n",
    "ax.set_ylim(-10,60)\n",
    "mpl_functions.adjust_spines(ax,['bottom'],\n",
    "                            spine_locations={'bottom':4},\n",
    "                            smart_bounds=True,\n",
    "#                             yticks=[-10,60],\n",
    "                            xticks=[2200,2500],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([], rotation=90)\n",
    "\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f089557e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('klm','lws2')]\n",
    "# f,ax=plt.subplots(1,1)\n",
    "ax.grid(False)\n",
    "error = xhat_smooth[0,:] - data_for_kf.distance_true.values\n",
    "h = ax.hist( np.abs(error), bins=100, color='gray',histtype='stepfilled', rasterized=True)\n",
    "median_error = np.median(np.abs(error))\n",
    "\n",
    "ax.axvline(median_error, 0, ymax=0.5,color='blue', linewidth=line_width)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "ax.set_ylim(0,3000)\n",
    "\n",
    "ax.set_xlim(0,20)\n",
    "# ax.text(14,4000,'Median\\nError=2.92',color='black')\n",
    "mpl_functions.adjust_spines(ax,['bottom'],\n",
    "                            spine_locations={'left':4,'bottom':4},\n",
    "                            smart_bounds=True,\n",
    "                            yticks=[0,3000],\n",
    "                            xticks=[0,20],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "ax.set_xlabel('20m')\n",
    "ax.xaxis.set_label_coords(0.5, -0.15)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "33b8d364",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.261460282768304"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "median_error"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32893b81",
   "metadata": {},
   "source": [
    "### Forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "0188e2d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir='~/DataAnalysis/data/Sprints/HighRes/KF/'\n",
    "data_1=pd.read_hdf(dir+'AllpredictionFiltered.h5')\n",
    "data_2=pd.read_hdf(dir+'ForestMASigned.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "2c4b2931",
   "metadata": {},
   "outputs": [],
   "source": [
    "xhat_smooth,time,data_for_kf = klm_cal_forest(data_1,data_2)\n",
    "\n",
    "ax=layout.axes[('klm','forest')]\n",
    "ax.grid(False)\n",
    "line_width=0.5\n",
    "marker_size=0.5\n",
    "ax.plot(data_1.mean_whiff_time, data_1.predicted_distance, '.', \n",
    "        color='black', label='Predicted Distance',linewidth=line_width,\n",
    "        markersize=marker_size*2,\n",
    "        rasterized= True)\n",
    "ax.plot(time, xhat_smooth[0,:], color='black', label='Smoothed Distance', \n",
    "        linewidth=line_width,markersize=marker_size,\n",
    "        rasterized= True)\n",
    "ax.plot(data_for_kf.time, data_for_kf.distance_true, '--', color='red',\n",
    "        label='True Distance',linewidth=line_width, markersize=marker_size,rasterized= True)\n",
    "\n",
    "ax.set_rasterization_zorder(-20)\n",
    "ax.set_ylim(-10,70)\n",
    "ax.set_xlim(0,1000)\n",
    "mpl_functions.adjust_spines(ax,['bottom'],\n",
    "                            spine_locations={'bottom':4},\n",
    "                            smart_bounds=True,\n",
    "                            yticks=[-10,70],\n",
    "                            xticks=[0,1000],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([], rotation=90)\n",
    "ax.axvspan(xmin=0, xmax=300,  alpha=0.4, color='grey')\n",
    "# ax.set_xlabel('Time, s')\n",
    "# ax.xaxis.set_label_coords(0.5, -0.1)\n",
    "\n",
    "# ax.set_title('Forest', x=0.5, y=0.95)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "57d1af18",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('klm','forest1')]\n",
    "ax.grid(False)\n",
    "line_width=0.5\n",
    "marker_size=0.5\n",
    "ax.plot(data_1.mean_whiff_time, data_1.predicted_distance, '.', \n",
    "        color='black', label='Predicted Distance',linewidth=line_width,\n",
    "        markersize=marker_size*2,\n",
    "        rasterized= True)\n",
    "ax.plot(time, xhat_smooth[0,:], color='black', label='Smoothed Distance', \n",
    "        linewidth=line_width,markersize=marker_size,\n",
    "        rasterized= True)\n",
    "ax.plot(data_for_kf.time, data_for_kf.distance_true, '--', color='red',\n",
    "        label='True Distance',linewidth=line_width, markersize=marker_size,rasterized= True)\n",
    "\n",
    "ax.set_rasterization_zorder(-20)\n",
    "ax.set_xlim(0,300)\n",
    "# ax.set_xlabel('Time, s')\n",
    "# ax.xaxis.set_label_coords(0.5, -0.1)\n",
    "\n",
    "# ax.set_ylabel('Distance, m')\n",
    "# ax.yaxis.set_label_coords(-0.08, 0.5)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['bottom'],\n",
    "                            spine_locations={'bottom':4},\n",
    "                            smart_bounds=True,\n",
    "#                             yticks=[-10,60],\n",
    "                            xticks=[0,300],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([], rotation=90)\n",
    "# ax.set_xlabel('Time, s')\n",
    "# ax.set_ylabel('Distance, m')\n",
    "# ax.set_title('WS>3.5m/s')\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "1b3dab3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('klm','forest2')]\n",
    "# f,ax=plt.subplots(1,1)\n",
    "ax.grid(False)\n",
    "error = xhat_smooth[0,:] - data_for_kf.distance_true.values\n",
    "\n",
    "h = ax.hist( np.abs(error), bins=100, color='gray',histtype='stepfilled', rasterized=True)\n",
    "\n",
    "median_error = np.median(np.abs(error))\n",
    "ax.set_ylim(0,3000)\n",
    "ax.axvline(median_error, 0, ymax=0.5,color='blue', linewidth=line_width)\n",
    "# ax.text(14,4000,'Median\\nError=13.03',color='black')\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "# ax.set_xlabel('Error, m')\n",
    "# ax.xaxis.set_label_coords(0.5, -0.1)\n",
    "\n",
    "# ax.set_ylabel('Density')\n",
    "# ax.yaxis.set_label_coords(-0.08, 0.5)\n",
    "ax.set_xlim(0,20)\n",
    "# ax.text(14,4000,'Median\\nError=2.92',color='black')\n",
    "mpl_functions.adjust_spines(ax,['bottom'],\n",
    "                            spine_locations={'left':4,'bottom':4},\n",
    "                            smart_bounds=True,\n",
    "                            yticks=[0,3000],\n",
    "                            xticks=[0,20],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "ax.set_xlabel('20m')\n",
    "ax.xaxis.set_label_coords(0.5, -0.15)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "183a1ead",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.767390691648034"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "median_error"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1358bfd",
   "metadata": {},
   "source": [
    "### Environment Specific Line"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fcb981c",
   "metadata": {},
   "source": [
    "### HWS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "e401fcb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir='~/DataAnalysis/data/Sprints/HighRes/KF/'\n",
    "data_1=pd.read_hdf(dir+'WindyPredictions.h5')\n",
    "data_2=pd.read_hdf(dir+'WindyMASigned.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "7612f22b",
   "metadata": {},
   "outputs": [],
   "source": [
    "xhat_smooth,time,data_for_kf = klm_cal(data_1,data_2)\n",
    "\n",
    "ax=layout.axes[('klm','hws')]\n",
    "ax.grid(False)\n",
    "line_width=0.5\n",
    "marker_size=0.5\n",
    "ax.plot(time, xhat_smooth[0,:], color='g', label='Smoothed Distance', \n",
    "        linewidth=line_width,markersize=marker_size,\n",
    "        rasterized= True)\n",
    "\n",
    "\n",
    "ax.set_rasterization_zorder(-20)\n",
    "ax.set_xlim(1000, 2000)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left','bottom'],\n",
    "                            spine_locations={'left':4,'bottom':4},\n",
    "                            smart_bounds=True,\n",
    "                            yticks=[-10,70],\n",
    "                            xticks=[1000,2000],\n",
    "                            linewidth=0.5)\n",
    "\n",
    "ax.set_xticklabels([])\n",
    "ax.axvspan(xmin=1000, xmax=1300, alpha=0.4, color='grey')\n",
    "ax.set_xlabel('1000 secs')\n",
    "ax.xaxis.set_label_coords(0.5, -0.1)\n",
    "\n",
    "ax.set_ylabel('Distance, m')\n",
    "ax.yaxis.set_label_coords(-0.08, 0.5)\n",
    "# ax.set_title('WS>3.5m/s', x=0.5, y=0.95)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "a56db51d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('klm','hws1')]\n",
    "ax.grid(False)\n",
    "line_width=0.8\n",
    "marker_size=0.5\n",
    "\n",
    "\n",
    "ax.plot(time, xhat_smooth[0,:], color='g', label='Smoothed Distance', \n",
    "        linewidth=line_width,markersize=marker_size,\n",
    "        rasterized= True)\n",
    "\n",
    "ax.set_rasterization_zorder(-20)\n",
    "ax.set_xlim(1000,1300)\n",
    "ax.set_xlabel('300 secs')\n",
    "ax.xaxis.set_label_coords(0.5, -0.1)\n",
    "\n",
    "ax.set_ylabel('Distance, m')\n",
    "ax.yaxis.set_label_coords(-0.08, 0.5)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left','bottom'],\n",
    "                            spine_locations={'bottom':4, 'left':4},\n",
    "                            smart_bounds=True,\n",
    "                            yticks=[-10,60],\n",
    "                            xticks=[1000,1300],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([], rotation=90)\n",
    "\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0bb59d0",
   "metadata": {},
   "source": [
    "### LWS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "890bbccd",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir='~/DataAnalysis/data/Sprints/HighRes/KF/'\n",
    "data_1=pd.read_hdf(dir+'NotWindyPredictions.h5')\n",
    "data_2=pd.read_hdf(dir+'NotWindyMASigned.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "e74ce285",
   "metadata": {},
   "outputs": [],
   "source": [
    "xhat_smooth,time,data_for_kf = klm_cal(data_1,data_2)\n",
    "\n",
    "ax=layout.axes[('klm','lws')]\n",
    "ax.grid(False)\n",
    "line_width=0.8\n",
    "marker_size=0.5\n",
    "\n",
    "ax.plot(time, xhat_smooth[0,:], color='g', label='Smoothed Distance', \n",
    "        linewidth=line_width,markersize=marker_size,\n",
    "        rasterized= True)\n",
    "\n",
    "ax.set_rasterization_zorder(-20)\n",
    "ax.set_xlim(2000,3000)\n",
    "ax.set_ylim(-10,70)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['bottom'],\n",
    "                            spine_locations={'bottom':4},\n",
    "                            smart_bounds=True,\n",
    "#                             yticks=[-10,70],\n",
    "                            xticks=[2000,3000],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([], rotation=90)\n",
    "ax.axvspan(xmin=2200, xmax=2500,alpha=0.4, color='grey')\n",
    "# ax.set_xlabel('1000, s')\n",
    "# ax.xaxis.set_label_coords(0.5, -0.1)\n",
    "\n",
    "# ax.set_title('WS<3.5m/s', x=0.5, y=0.95)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "55d735fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('klm','lws1')]\n",
    "ax.grid(False)\n",
    "line_width=0.8\n",
    "marker_size=0.5\n",
    "\n",
    "ax.plot(time, xhat_smooth[0,:], color='g', label='Smoothed Distance', \n",
    "        linewidth=line_width,markersize=marker_size,\n",
    "        rasterized= True)\n",
    "\n",
    "ax.set_rasterization_zorder(-20)\n",
    "ax.set_xlim(2200,2500)\n",
    "ax.set_ylim(-10,60)\n",
    "mpl_functions.adjust_spines(ax,['bottom'],\n",
    "                            spine_locations={'bottom':4},\n",
    "                            smart_bounds=True,\n",
    "#                             yticks=[-10,60],\n",
    "                            xticks=[2200,2500],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([], rotation=90)\n",
    "\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd16f1f5",
   "metadata": {},
   "source": [
    "### Forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "9e959014",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir='~/DataAnalysis/data/Sprints/HighRes/KF/'\n",
    "data_1=pd.read_hdf(dir+'ForestPredictions.h5')\n",
    "data_2=pd.read_hdf(dir+'ForestMASigned.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "5f8b6fe1",
   "metadata": {},
   "outputs": [],
   "source": [
    "xhat_smooth,time,data_for_kf = klm_cal_forest(data_1,data_2)\n",
    "\n",
    "ax=layout.axes[('klm','forest')]\n",
    "ax.grid(False)\n",
    "line_width=0.8\n",
    "marker_size=0.5\n",
    "\n",
    "ax.plot(time, xhat_smooth[0,:], color='g', label='Smoothed Distance', \n",
    "        linewidth=line_width,markersize=marker_size,\n",
    "        rasterized= True)\n",
    "\n",
    "ax.set_rasterization_zorder(-20)\n",
    "ax.set_ylim(-10,70)\n",
    "ax.set_xlim(0,1000)\n",
    "mpl_functions.adjust_spines(ax,['bottom'],\n",
    "                            spine_locations={'bottom':4},\n",
    "                            smart_bounds=True,\n",
    "                            yticks=[-10,70],\n",
    "                            xticks=[0,1000],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([], rotation=90)\n",
    "ax.axvspan(xmin=0, xmax=300,  alpha=0.4, color='grey')\n",
    "\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "dc07b46b",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('klm','forest1')]\n",
    "ax.grid(False)\n",
    "line_width=0.8\n",
    "marker_size=0.5\n",
    "\n",
    "ax.plot(time, xhat_smooth[0,:], color='g', label='Smoothed Distance', \n",
    "        linewidth=line_width,markersize=marker_size,\n",
    "        rasterized= True)\n",
    "\n",
    "\n",
    "ax.set_rasterization_zorder(-20)\n",
    "ax.set_xlim(0,300)\n",
    "\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['bottom'],\n",
    "                            spine_locations={'bottom':4},\n",
    "                            smart_bounds=True,\n",
    "#                             yticks=[-10,60],\n",
    "                            xticks=[0,300],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([], rotation=90)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['klm'], 'klm', cleartarget=True)\n",
    "layout.write_svg('../../../Figure/fig7.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3a87f9f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "vscode38kernel",
   "language": "python",
   "name": "pywork38"
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}