Step3.AnalyseData.scipy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Video Actor Synchroncy and Causality (VASC)\n",
    "## RAEng: Measuring Responsive Caregiving Project\n",
    "### Caspar Addyman, 2020\n",
    "### https://github.com/infantlab/VASC\n",
    "\n",
    "# Step 3: Analyse the data using scipy statsmodels\n",
    "\n",
    "This script correlates and compares the timeseries of wireframes for the two figures in the video `[\"parent\", \"infant\"]`\n",
    "\n",
    "We start by reloading the saved parquet file containing the multi-index numpy array of all [OpenPose](https://github.com/CMU-Perceptual-Computing-Lab/openpose) data from all pairs of individuals.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Automatic pdb calling has been turned ON\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import os\n",
    "import json\n",
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pyarrow.parquet as pq\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import logging\n",
    "import ipywidgets as widgets  #let's us add buttons and sliders to this page.\n",
    "from ipycanvas import Canvas\n",
    "\n",
    "import vasc #a module of our own functions (found in vasc.py in this folder)\n",
    "\n",
    "#turn on debugging\n",
    "logger = logging.getLogger()\n",
    "logger.setLevel(logging.INFO)\n",
    "%pdb on"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Settings?\n",
    "\n",
    "Load a json file that tells us where to find our videos and where to save the data. You should create a different settings file for each project. Then you don't need to change any other values in the script for Step 1 or Step 2.\n",
    "\n",
    "TODO - write a helper to create a settings file\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Existing settings.json found..\n"
     ]
    }
   ],
   "source": [
    "settingsjson = \".\\\\Drum.Tutorial.settings.json\"\n",
    "\n",
    "try:\n",
    "    with open(settingsjson) as json_file:\n",
    "        settings = json.load(json_file)\n",
    "        print(\"Existing settings.json found..\")\n",
    "except json.JSONDecodeError:\n",
    "    logging.exception(\"Settings file was not valid JSON.\")\n",
    "except Exception as e:\n",
    "        emsg = str(e)\n",
    "        #show the error\n",
    "        print(\"Error: \",emsg)\n",
    "        print(\"No setting.json file found!\\nPlease see Step 0 for instructions\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "includeHands = settings[\"flags\"][\"includeHands\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "anon = settings[\"flags\"][\"anon\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "videos_in:  .\\DrumTutorial\\videos\n",
      "videos_out:  .\\DrumTutorial\\\n",
      "videos_out_openpose:  .\\DrumTutorial\\openpose\n",
      "videos_out_timeseries:  .\\DrumTutorial\\timeseries\n",
      "videos_out_analyses:  .\\DrumTutorial\\analyses\n"
     ]
    }
   ],
   "source": [
    "# where's the project data folder? (with trailing slash)\n",
    "projectpath = settings[\"paths\"][\"project\"]\n",
    "#where are your video files?\n",
    "videos_in = settings[\"paths\"][\"videos_in\"]\n",
    "\n",
    "# locations of videos and output\n",
    "videos_out = settings[\"paths\"][\"videos_out\"]\n",
    "videos_out_openpose   = settings[\"paths\"][\"videos_out_openpose\"]\n",
    "videos_out_timeseries = settings[\"paths\"][\"videos_out_timeseries\"]\n",
    "videos_out_analyses   = settings[\"paths\"][\"videos_out_analyses\"]\n",
    "\n",
    "\n",
    "print(\"videos_in: \", videos_in)\n",
    "print(\"videos_out: \", videos_out)\n",
    "print(\"videos_out_openpose: \", videos_out_openpose)\n",
    "print(\"videos_out_timeseries: \", videos_out_timeseries)\n",
    "print(\"videos_out_analyses: \", videos_out_analyses)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Load the clean data as a DataFrame\n",
    "\n",
    "Reload the clean data file created in step 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Existing clean.json found..\n"
     ]
    }
   ],
   "source": [
    "#retrieve the list of base names of processed videos.\n",
    "videosjson = settings[\"paths\"][\"videos_out\"] + '\\\\' + settings[\"filenames\"][\"clean_json\"]\n",
    "try:\n",
    "    with open(videosjson) as json_file:\n",
    "        videos = json.load(json_file)\n",
    "        print(\"Existing clean.json found..\")\n",
    "except:\n",
    "    videos = {}\n",
    "    print(\"No clean.json file found, please locate the file or complete Step 2 first.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reading parquet file:\n",
      "video  1fa339b_04-test-trials                                                 \\\n",
      "person                 infant                                                  \n",
      "coord                      0        1         2        3        4         5    \n",
      "0                     603.999  402.094  0.812414  619.638  488.323  0.591055   \n",
      "1                     605.907  400.237  0.817730  619.558  490.241  0.580483   \n",
      "2                     607.815  400.268  0.798972  621.568  490.254  0.570408   \n",
      "3                     605.932  402.122  0.807359  617.672  484.432  0.616228   \n",
      "4                     605.928  402.124  0.807981  617.688  482.488  0.613706   \n",
      "\n",
      "video                                        ... b22644a_14-test-trials       \\\n",
      "person                                       ...                 parent        \n",
      "coord        6        7         8        9   ...                     65   66   \n",
      "0       478.519  484.430  0.513021  374.709  ...                    0.0  0.0   \n",
      "1       476.577  488.301  0.516335  372.761  ...                    0.0  0.0   \n",
      "2       478.491  486.372  0.500429  376.642  ...                    0.0  0.0   \n",
      "3       476.522  484.420  0.517500  380.506  ...                    0.0  0.0   \n",
      "4       476.541  482.474  0.516069  380.512  ...                    0.0  0.0   \n",
      "\n",
      "video                                           \n",
      "person                                          \n",
      "coord    67   68   69   70   71   72   73   74  \n",
      "0       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  \n",
      "1       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  \n",
      "2       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  \n",
      "3       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  \n",
      "4       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  \n",
      "\n",
      "[5 rows x 4500 columns]\n"
     ]
    }
   ],
   "source": [
    "print('reading parquet file:')\n",
    "df = pq.read_table(videos_out_timeseries + '\\\\' + settings[\"filenames\"][\"cleandataparquet\"]).to_pandas()\n",
    "if includeHands:\n",
    "    lh = pq.read_table(videos_out_timeseries + '\\\\' + settings[\"filenames\"][\"lefthandparquet\"]).to_pandas()\n",
    "    rh = pq.read_table(videos_out_timeseries + '\\\\' + settings[\"filenames\"][\"righthandparquet\"]).to_pandas()\n",
    "\n",
    "#sort the column names as this helps with indexing\n",
    "df = df.sort_index(axis = 1)\n",
    "print(df.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 Process the data\n",
    "\n",
    "Next we set all 0 values to as missing value `np.nan` to enable interpolation.\n",
    "Then use numpy's built in `interpolate` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.replace(0.0, np.nan)\n",
    "if includeHands:\n",
    "    rh = rh.replace(0.0, np.nan)\n",
    "    lh = lh.replace(0.0, np.nan)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(956, 4500)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#are we going to use all the data or a subset?\n",
    "first = 0\n",
    "last = 4500\n",
    "\n",
    "df = df.truncate(before  = first, after = last)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#linear interpolate missing values\n",
    "df = df.interpolate()\n",
    "#may still have some NaNs at start so back fill these with first non-zero value\n",
    "#eg  [NaN, NaN, 3.1, 3.2, ...] -> [3.1, 3.1, 3.1, 3.2, ...]\n",
    "df = df.fillna(method = 'backfill')\n",
    "\n",
    "if includeHands:\n",
    "    rh = rh.interpolate()\n",
    "    rh = rh.fillna(method = 'backfill')\n",
    "    lh = lh.interpolate()\n",
    "    lh = lh.fillna(method = 'backfill')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "video  1fa339b_04-test-trials                                                 \\\n",
      "person                 infant                                                  \n",
      "coord                      0        1         2        3        4         5    \n",
      "0                     603.999  402.094  0.812414  619.638  488.323  0.591055   \n",
      "1                     605.907  400.237  0.817730  619.558  490.241  0.580483   \n",
      "2                     607.815  400.268  0.798972  621.568  490.254  0.570408   \n",
      "3                     605.932  402.122  0.807359  617.672  484.432  0.616228   \n",
      "4                     605.928  402.124  0.807981  617.688  482.488  0.613706   \n",
      "\n",
      "video                                        ... b22644a_14-test-trials      \\\n",
      "person                                       ...                 parent       \n",
      "coord        6        7         8        9   ...                     65  66   \n",
      "0       478.519  484.430  0.513021  374.709  ...                    NaN NaN   \n",
      "1       476.577  488.301  0.516335  372.761  ...                    NaN NaN   \n",
      "2       478.491  486.372  0.500429  376.642  ...                    NaN NaN   \n",
      "3       476.522  484.420  0.517500  380.506  ...                    NaN NaN   \n",
      "4       476.541  482.474  0.516069  380.512  ...                    NaN NaN   \n",
      "\n",
      "video                                   \n",
      "person                                  \n",
      "coord   67  68  69  70  71  72  73  74  \n",
      "0      NaN NaN NaN NaN NaN NaN NaN NaN  \n",
      "1      NaN NaN NaN NaN NaN NaN NaN NaN  \n",
      "2      NaN NaN NaN NaN NaN NaN NaN NaN  \n",
      "3      NaN NaN NaN NaN NaN NaN NaN NaN  \n",
      "4      NaN NaN NaN NaN NaN NaN NaN NaN  \n",
      "\n",
      "[5 rows x 4500 columns]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(956, 4500)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#take a quick look\n",
    "print(df.head())\n",
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(lh.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2.1 Mean movements\n",
    "We create a dictionary of the subsets of OpenPose coordinates we want to average and then call `mean` on the Pandas dataframe. e.g.\n",
    "\n",
    "```\n",
    "meanpoints = {\n",
    "               \"headx\" : [0, 3, 45, 48, 51, 54],\n",
    "               \"heady\" : [1, 4, 46, 49, 52, 55],\n",
    "               \"allx\" :  [0, 3, 6, 9, ...],\n",
    "               \"ally\" :  [1, 4, 7, 10, ...]\n",
    "             }\n",
    "```\n",
    "\n",
    "Then we call the `vasc.averageCoordinateTimeSeries` function to average across sets of coordinates. For a given set of videos and people. For example\n",
    "\n",
    "In:\n",
    "```\n",
    "videos = \"All\"\n",
    "people = \"Both\"\n",
    "df2 = vasc.averageCoordinateTimeSeries(df,meanpoints,videos,people)\n",
    "df2.head\n",
    "```\n",
    "\n",
    "Out:\n",
    "```\n",
    "person      infant                                          parent\n",
    "avgs         headx       heady          xs          ys       headx\n",
    "501     565.996600  369.840600  534.895615  398.482538  471.686200\n",
    "502     567.231800  369.887600  534.354198  398.706552  471.849400\n",
    "503     567.228600  370.159600  534.444328  398.678133  471.711600\n",
    "504     566.912600  369.857000  535.369536  398.551636  472.309400\n",
    "...            ...         ...         ...         ...         ...\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "meanpoints = {\"head\" : vasc.headxys,\n",
    "              \"headx\": vasc.headx,\n",
    "              \"heady\": vasc.heady,\n",
    "              \"arms\" : vasc.armsxys,\n",
    "              \"armsx\": vasc.armsx,\n",
    "              \"armsy\": vasc.armsy,\n",
    "              \"leftarm\" : vasc.leftarmxys,\n",
    "              \"leftarmx\": vasc.leftarmx,\n",
    "              \"leftarmy\": vasc.leftarmy,\n",
    "              \"rightarm\" : vasc.rightarmxys,\n",
    "              \"rightarmx\": vasc.rightarmx,\n",
    "              \"rightarmy\": vasc.rightarmy,\n",
    "              \"all\"  : vasc.xys,\n",
    "              \"allx\" : vasc.xs,\n",
    "              \"ally\" : vasc.ys\n",
    "             }\n",
    "\n",
    "vids = \"All\"\n",
    "people = [\"infant\",\"parent\"]\n",
    "\n",
    "#average across the points in each group (all points of head etc. )\n",
    "avgdf = vasc.averageCoordinateTimeSeries(df,meanpoints,vids,people)\n",
    "\n",
    "hps = {\"hand\" : vasc.hxys,\n",
    "       \"handx\" : vasc.hxs,\n",
    "       \"handy\" : vasc.hys}\n",
    "\n",
    "avglh = vasc.averageCoordinateTimeSeries(lh,hps,vids,people)\n",
    "avgrh = vasc.averageCoordinateTimeSeries(rh,hps,vids,people)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "avgdf.head"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "%pdb on"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2.2 Rolling window of movements\n",
    "\n",
    "One thing we'd like to know is if mothers move in response to infants. The raw time series are probably too noisy to tell us this so instead we can look at few alternatives\n",
    "\n",
    "1. **Smoothed** - if we average the signal over a short rolling window we smooth out any high-frequency jitter.\n",
    "2. **Variance** - the variance of movement over a short rolling window. First we apply short (10 frame) rolling window to each coordinate of the body and use the stddev or variance function `std()` or `var()` . Then we take averages as in the step above. However, this time we combine x and y coordinates as this is now a movement index.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "win = 10 #10 frames better for rhythm detcion\n",
    "halfwin = math.floor(win/2)\n",
    "\n",
    "smoothdf = df.rolling(window = 5).mean()\n",
    "smoothdf = smoothdf.truncate(before  = first, after = last)\n",
    "\n",
    "vardf = df.rolling(window = win, min_periods = halfwin).var()\n",
    "vardf = vardf.truncate(before  = first , after = last) # cut out the empty bits at the start\n",
    "\n",
    "smoothdf = vasc.averageCoordinateTimeSeries(smoothdf,meanpoints,vids,people)\n",
    "vardf = vasc.averageCoordinateTimeSeries(vardf,meanpoints,vids,people)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualising the data\n",
    "\n",
    "Let's create a widget to plot some graphs of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "458d5d780a7148d9b892e8d1456b364f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "vidlist = [] #used to fill dropdown options\n",
    "for vid in videos:\n",
    "    vidlist.append(vid)\n",
    "\n",
    "pickvid = widgets.Dropdown(\n",
    "    options= vidlist,\n",
    "    value= vidlist[0],\n",
    "    description='Subject:'\n",
    ")\n",
    "\n",
    "features = []\n",
    "for f in meanpoints:\n",
    "    features.append(f)\n",
    "\n",
    "pickfeature = widgets.Dropdown(\n",
    "    options= features,\n",
    "    value= features[0],\n",
    "    description='Feature:'\n",
    ")\n",
    "\n",
    "linetypes = [\"Mean point\", \"Smoothed Mean (5 frames)\",\"Variance over 2 secs\"]\n",
    "picktype = widgets.Dropdown(\n",
    "    options= linetypes,\n",
    "    value= linetypes[0],\n",
    "    description='Line type:'\n",
    ")\n",
    "\n",
    "def pickvid_change(change):\n",
    "    if change['name'] == 'value' and (change['new'] != change['old']):\n",
    "        updateAll(True)\n",
    "\n",
    "def pickfeature_change(change):\n",
    "    if change['name'] == 'value' and (change['new'] != change['old']):\n",
    "        updateAll(True)\n",
    "\n",
    "def picktype_change(change):\n",
    "    if change['name'] == 'value' and (change['new'] != change['old']):\n",
    "        updateAll(True)\n",
    "\n",
    "pickvid.observe(pickvid_change, 'value')\n",
    "pickfeature.observe(pickfeature_change, 'value')\n",
    "picktype.observe(picktype_change, 'value')\n",
    "button_update = widgets.Button(description=\"Redraw\")\n",
    "output = widgets.Output()\n",
    "\n",
    "\n",
    "def drawGraphs(vid, feature, linetype):\n",
    "    \"\"\"Plot input signals\"\"\"\n",
    "    plt.ion()\n",
    "\n",
    "    f,ax=plt.subplots(4,1,figsize=(14,10),sharex=True)\n",
    "    ax[0].set_title('Infant')\n",
    "    ax[1].set_title('Parent')\n",
    "    ax[1].set_xlabel('Frames')\n",
    "\n",
    "    who = [\"infant\",\"parent\"]\n",
    "\n",
    "    if linetype == linetypes[0]:\n",
    "        usedf = avgdf\n",
    "    elif linetype == linetypes[1]:\n",
    "        usedf = smoothdf\n",
    "    else:\n",
    "        usedf = vardf\n",
    "\n",
    "    #to select a single column..\n",
    "    infant = usedf[(vid, people[0], feature)].to_frame()\n",
    "    parent = usedf[(vid, people[1], feature)].to_frame()\n",
    "    n  = np.arange(usedf.shape[0])\n",
    "\n",
    "    #selecting multiple columns slightly messier\n",
    "    #infant = df3.loc[50:,(vid, part[0], ('head','arms', 'all'))]\n",
    "    #parent = df3.loc[50:,(vid, part[1], ('head','arms', 'all'))]\n",
    "\n",
    "    ax[0].plot(n,infant)\n",
    "    ax[1].plot(n,parent, color='b')\n",
    "\n",
    "    #calculate the correlations in a shorter rolling window\n",
    "    r_window_size = 120\n",
    "    rolling_r = usedf[(vid, who[0], feature)].rolling(window=r_window_size, center=True).corr(vardf[(vid, who[1], feature)])\n",
    "\n",
    "\n",
    "    usedf.loc[:,(vid, slice(None), feature)].plot(ax=ax[2])\n",
    "    ax[2].set(xlabel='Frame',ylabel='Movement index for parent and infant')\n",
    "\n",
    "    rolling_r.plot(ax=ax[3])\n",
    "    ax[3].set(xlabel='Frame',ylabel='Pearson r')\n",
    "    ax[3].set_title(\"Local correlation with rolling window size \" + str(r_window_size))\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "def updateAll(forceUpdate = False):\n",
    "    output.clear_output(wait = True)\n",
    "    if forceUpdate:\n",
    "        logging.debug('forceUpdate')\n",
    "        #slider.value = 0\n",
    "        #slider.max = videos[pickvid.value][pickcam.value][\"end\"]\n",
    "    with output:\n",
    "        display(pickvid,pickfeature,picktype,button_update)\n",
    "        drawGraphs(pickvid.value,pickfeature.value,picktype.value)\n",
    "\n",
    "#draw everything for first time\n",
    "updateAll(True)\n",
    "output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Fitting the best periodic function.\n",
    "\n",
    "For each infant and each trial we try to find the best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "who = [\"infant\",\"parent\"]\n",
    "parts = [\"head\",\"arms\",\"all\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'results' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_180/147969753.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     10\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     11\u001b[0m     \u001b[0mthisrow\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m#don't have DyadSynScore yet\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 12\u001b[1;33m     \u001b[0mresults\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mvid\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mthisrow\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'results' is not defined"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "> \u001b[1;32mc:\\users\\cas\\appdata\\local\\temp\\ipykernel_180\\147969753.py\u001b[0m(12)\u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n",
      "\n",
      "ipdb> q\n"
     ]
    }
   ],
   "source": [
    "#loop through colculate for each pair\n",
    "for vid in videos:\n",
    "    thisrow = []\n",
    "    for part in parts:\n",
    "        #to select a single column..\n",
    "        pearson = vardf[(vid, people[0], part)].corr(vardf[(vid, people[1], part)])\n",
    "\n",
    "        thisrow.append(pearson) #this is for correlation\n",
    "        thisrow.append(None) #this is for maximum lag\n",
    "\n",
    "    thisrow.append(None) #don't have DyadSynScore yet\n",
    "    results.loc[vid] = thisrow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting the best periodic function\n",
    "\n",
    "We use the scipy [least squares optimiser](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html#scipy.optimize.curve_fit) to find the best fitting sine wave, estimating phase, frequency and amplitude.\n",
    "\n",
    "Because each video could have a different frame rate (fps) we have to make a note of that each time.\n",
    "\n",
    "This page provides a quick example\n",
    "https://scipy-lectures.org/intro/scipy/auto_examples/plot_curve_fit.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import scipy\n",
    "import numpy as np\n",
    "from scipy import optimize\n",
    "import pylab as plt\n",
    "\n",
    "\n",
    "#the function that we want to fit to our data.\n",
    "#basically a sine wave where we optimize everything - especially frequency.\n",
    "def optimize_func(x_time, amp, freq, phase, mean):\n",
    "    return amp * np.sin(2 * math.pi * freq * x_time + phase) + mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "plotgraphs = True\n",
    "showfirstguess = False\n",
    "savedata = True\n",
    "GuessSMT  = 450  #we expect infants to have a spontaneous ISI of 450 ms\n",
    "GuessFreq = 1000 / GuessSMT\n",
    "\n",
    "bodypart = 'leftarmy'\n",
    "\n",
    "failarray = [math.nan, math.nan, math.nan, math.nan,math.nan]  #empty array if we fail to fit function\n",
    "\n",
    "\n",
    "results = pd.DataFrame(columns = (\"fit measure\", \"best_amp\",\"best_freq\",\"best_phase\",\"best_mean\",\"fps\",\"fitted_isi\",\"message\"), index = videos)\n",
    "\n",
    "for vid in videos:\n",
    "    #retrieve the raw data for this video\n",
    "    frames = videos[vid]['camera1']['frames']    #how many frames?\n",
    "    fps = videos[vid]['camera1']['fps']          #how many frames per second?\n",
    "    x_data = np.linspace(0,frames,frames+1)      #x axis\n",
    "    x_time = x_data / fps                        #x axis in units of time (seconds)\n",
    "    armmov = avgdf[(vid, people[0], bodypart)]  #get the averaged data\n",
    "    y_data = armmov.truncate(before = 0, after = frames).to_numpy()  #convert to numpy array\n",
    "\n",
    "    try:\n",
    "        #have a first guess of the parameters\n",
    "        guess_mean = float(np.mean(armmov))\n",
    "        guess_amp =  2* float(np.std(armmov))\n",
    "        guess_phase = 0.0\n",
    "        guess_freq = GuessFreq\n",
    "\n",
    "        # we'll use this to plot our first estimate. This might already be good enough for you\n",
    "        data_first_guess = optimize_func(x_time, guess_amp, guess_freq, guess_phase, guess_mean)\n",
    "\n",
    "        #now optimize\n",
    "        params, params_covariance = optimize.curve_fit(optimize_func, x_time, y_data,\n",
    "                                                   p0=[guess_amp, guess_freq,guess_phase,guess_mean])\n",
    "\n",
    "        #what is the best fit?\n",
    "        fitted = optimize_func(x_time, params[0],params[1],params[2],params[3])\n",
    "\n",
    "        #how far is this from raw data?\n",
    "        best_fit = scipy.spatial.distance.euclidean(y_data,fitted )\n",
    "\n",
    "        if savedata:\n",
    "            #save all values to dataframe\n",
    "            vals = [best_fit]\n",
    "            vals.extend(params)\n",
    "            vals.append(fps)\n",
    "            vals.append(1000/(params[1]))  #fitted isi\n",
    "            vals.append(\"Success\")\n",
    "            print(vals)\n",
    "            results.loc[vid] = vals\n",
    "\n",
    "        if plotgraphs:\n",
    "            #plot the functions\n",
    "            plt.figure(figsize=(6, 4),)\n",
    "            plt.plot(x_time, y_data, label='Data')\n",
    "            if showfirstguess:\n",
    "                plt.plot(x_time, data_first_guess, label = 'First guess')\n",
    "            plt.plot(x_time, fitted,label='Fitted function')\n",
    "            plt.title(vid + \" \" + bodypart)\n",
    "            plt.legend(loc='best')\n",
    "\n",
    "            plt.show()\n",
    "\n",
    "    except Exception as e:\n",
    "        #emsg = \"\".join(str(e)) #error message string hack :(\n",
    "        emsg = str(e)\n",
    "        #show the error\n",
    "        print(vid, \" Error: \",emsg)\n",
    "        if savedata:\n",
    "            #record error in our results array\n",
    "            vals = failarray.copy()\n",
    "            vals.append(fps)\n",
    "            vals.append(math.nan)\n",
    "            vals.append(emsg)\n",
    "            print(vals)\n",
    "            results.loc[vid] = vals\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save the fitted parameters.\n",
    "results.to_excel(\"Test.LeftArm.FreeFit.xlsx\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting sine waves with fixed frequencies\n",
    "\n",
    "The next thing we try is to try fitting a sine wave with the expected frequency for that trial. We have two conditions A & B.\n",
    "\n",
    "```\n",
    "                      Condition A      Condition B\n",
    "xxxx_04-test-trials      Spontaneous Motor Tempo\n",
    "xxxx_06-test-trials      700ms ISI     400ms ISI\n",
    "xxxx_08-test-trials      500ms ISI     600ms ISI\n",
    "xxxx_10-test-trials      600ms ISI     500ms ISI\n",
    "xxxx_12-test-trials      400ms ISI     700ms ISI\n",
    "xxxx_14-test-trials      Spontaneous Motor Tempo\n",
    "```\n",
    "\n",
    "Therefore, for a given trial type we test the two possible fixed frequency/ISI values to see which fits best.\n",
    "We do not include the Spontaneous Motor Tempo trials in this analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "\n",
    "plotgraphs = True\n",
    "showfirstguess = False\n",
    "savedata = True\n",
    "GuessSMT  = 450  #we expect infants to have a spontaneous ISI of 450 ms\n",
    "GuessFreq = 1000/GuessSMT\n",
    "\n",
    "bodypart = 'rightarmy'\n",
    "\n",
    "failarray = [math.nan, math.nan, math.nan, math.nan,math.nan]  #empty array if we fail to fit function\n",
    "\n",
    "\n",
    "def optimize_ISI(fixeddata, amp, phase, mean):\n",
    "#we need function that we want to fit to our data.\n",
    "#basically a sine wave where we optimize everything - especially frequency.\n",
    "#this is slightly complex because each video could have differnt fps.\n",
    "#so the fixed data is now x_data, fps, ISI\n",
    "#we optimise on phase, amp & mean\n",
    "    isi = fixeddata[\"isi\"]\n",
    "    freq = 1000 / (isi)\n",
    "    x_time = fixeddata[\"x_time\"]\n",
    "    return amp * np.sin(twopi * freq * x_time + phase) + mean\n",
    "\n",
    "\n",
    "results = pd.DataFrame(columns = (\"fit value\",\"best_freq\",\"best_amp\",\"best_phase\",\"best_mean\",\"isi\",\"fps\", \"message\"))\n",
    "\n",
    "for vid in videos:\n",
    "    #retrieve the raw data for this video\n",
    "    frames = videos[vid]['camera1']['frames']    #how many frames?\n",
    "    fps = videos[vid]['camera1']['fps']          #how many frames per second?\n",
    "    x_data = np.linspace(0,frames,frames+1)      #x axis\n",
    "    x_time = x_data / fps                        #x axis in units of time (seconds)\n",
    "    fixeddata = {}\n",
    "    fixeddata[\"fps\"] = fps\n",
    "    fixeddata[\"x_time\"] = x_time\n",
    "    armmov = avgdf[(vid, people[0], bodypart)]  #get the averaged data\n",
    "    y_data = armmov.truncate(before = 0, after = frames).to_numpy()  #convert to numpy array\n",
    "\n",
    "    #WHat ISI\n",
    "    isi = [450, 450]\n",
    "    if '_06-' in vid:\n",
    "        isi[0] = 700\n",
    "        isi[1] = 400\n",
    "    elif '_08-' in vid:\n",
    "        isi[0] = 500\n",
    "        isi[1] = 600\n",
    "    elif '_10-' in vid:\n",
    "        isi[0] = 600\n",
    "        isi[1] = 500\n",
    "    elif '_12-' in vid:\n",
    "        isi[0] = 400\n",
    "        isi[1] = 700\n",
    "    else:\n",
    "        isi[0] = 450\n",
    "        isi[1] = 450\n",
    "\n",
    "    for iisi in isi:\n",
    "        try:\n",
    "            #have a first guess of the parameters\n",
    "            guess_mean = float(np.mean(armmov))\n",
    "            guess_amp =  2* float(np.std(armmov))\n",
    "            guess_phase = 0.0\n",
    "            guess_freq = 1000 / iisi  # freq in Hz\n",
    "\n",
    "            fixeddata[\"isi\"] = iisi\n",
    "\n",
    "            # we'll use this to plot our first estimate. This might already be good enough for you\n",
    "            data_first_guess = optimize_ISI(fixeddata, guess_amp, guess_phase, guess_mean)\n",
    "\n",
    "            #now optimize\n",
    "            params, params_covariance = optimize.curve_fit(optimize_ISI, fixeddata, y_data,\n",
    "                                                       p0=[guess_amp,guess_phase,guess_mean])\n",
    "\n",
    "            #what is the best fit?\n",
    "            fitted = optimize_ISI(fixeddata, params[0],params[1],params[2])\n",
    "\n",
    "            #how far is this from raw data?\n",
    "            best_fit = scipy.spatial.distance.euclidean(y_data,fitted )\n",
    "\n",
    "            if savedata:\n",
    "                #save all values to dataframe\n",
    "                vals = [best_fit]\n",
    "                vals.append(guess_freq)\n",
    "                vals.extend(params)\n",
    "                vals.append(iisi)\n",
    "                vals.append(fps)\n",
    "                vals.append(\"Success\")\n",
    "                print(vals)\n",
    "                results.loc[vid+ \"-\" + str(iisi)] = vals\n",
    "\n",
    "            if plotgraphs:\n",
    "                #plot the functions\n",
    "                plt.figure(figsize=(6, 4),)\n",
    "                plt.plot(x_time, y_data, label='Data')\n",
    "                plt.plot(x_time, fitted,label='Fitted function')\n",
    "                if showfirstguess:\n",
    "                    plt.plot(x_time, data_first_guess, label = 'First guess')\n",
    "                plt.title(vid + \" \" + str(iisi) + \"ms \" + bodypart )\n",
    "                plt.legend(loc='best')\n",
    "\n",
    "                plt.show()\n",
    "\n",
    "        except Exception as e:\n",
    "            #emsg = \"\".join(str(e)) #error message string hack :(\n",
    "            emsg = str(e)\n",
    "            #show the error\n",
    "            print(vid, \" Error: \",emsg)\n",
    "            if savedata:\n",
    "                #record error in our results array\n",
    "                vals = failarray.copy()\n",
    "                vals.append(iisi)\n",
    "                vals.append(fps)\n",
    "                vals.append(emsg)\n",
    "                print(vals)\n",
    "                results.loc[vid+ \"-\" + str(iisi)] = vals\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save the fitted parameters.\n",
    "results.to_excel(\"RightArm.FixedFit.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4 Finding fundamental frequency with FFT\n",
    "\n",
    "\n",
    "A good guide can be found here https://realpython.com/python-scipy-fft/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_end_of_cell_marker": 2,
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "import scipy.fft\n",
    "\n",
    "plotgraphs = True\n",
    "showfirstguess = False\n",
    "savedata = True\n",
    "GuessSMT  = 450  #we expect infants to have a spontaneous ISI of 450 ms\n",
    "GuessFreq = 1000/GuessSMT\n",
    "\n",
    "bodypart = 'rightarmy'\n",
    "\n",
    "failarray = [math.nan, math.nan, math.nan, math.nan,math.nan]  #empty array if we fail to fit function\n",
    "\n",
    "results = pd.DataFrame(columns = (\"peak\", \"400\",\"500\",\"600\",\"700\",\"fps\",\"message\"), index = videos)\n",
    "\n",
    "halfbinwidth = 4\n",
    "targetFreqs = [1000/400, 1000/500, 1000/600, 1000/700]\n",
    "freqResults =pd.DataFrame(columns = (\"ISI400\",\"ISI500\",\"ISI600\",\"ISI700\"))\n",
    "\n",
    "\n",
    "for vid in videos:\n",
    "    #retrieve the raw data for this video\n",
    "    frames = videos[vid]['camera1']['frames']    #how many frames?\n",
    "    fps = videos[vid]['camera1']['fps']          #how many frames per second?\n",
    "    x_data = np.linspace(0,frames,frames+1)      #x axis\n",
    "    x_time = x_data / fps                        #x axis in units of time (seconds)\n",
    "    armmov = avgdf[(vid, people[0], bodypart)]  #get the averaged data\n",
    "    y_data = armmov.truncate(before = 0, after = frames).to_numpy()  #convert to numpy array\n",
    "\n",
    "    #we are just interested in the periodic elements (not absolute value above zero) so substract the mean\n",
    "    y_normed = np.subtract(y_data,np.average(y_data))\n",
    "\n",
    "    try:\n",
    "        yf = scipy.fft.rfft(y_normed)\n",
    "        power = np.abs(yf)**2\n",
    "        xf = scipy.fft.rfftfreq(x_data.size, 1 / fps)\n",
    "\n",
    "        #let's find the maximum frequency.\n",
    "        #we hope this is infant's repetitive movement\n",
    "        mY = np.abs(power) # Find magnitude\n",
    "        peakY = np.max(mY) # Find max peak\n",
    "        locY = np.argmax(mY) # Find its location\n",
    "        frqY = xf[locY] # Get the actual frequency value\n",
    "        print(frqY)\n",
    "\n",
    "        #let's find the power in a given frequency bucket corresponding to target freq\n",
    "        powers = [0,0,0,0]\n",
    "        for p in range(4):\n",
    "            #find the nearest value in our list to\n",
    "            bucket = np.argmax(xf>targetFreqs[p])\n",
    "            powers[p] = sum(mY[bucket-halfbinwidth:bucket+halfbinwidth])\n",
    "\n",
    "        print(\"powers\")\n",
    "        print(powers)\n",
    "\n",
    "        if plotgraphs:\n",
    "            plt.figure(figsize=(12, 4),)\n",
    "            plt.subplot(1,2,1)\n",
    "            plt.plot(xf, power)\n",
    "\n",
    "            # New - Plot the max point\n",
    "            plt.plot(frqY, peakY, 'b.', markersize=18)\n",
    "            # Rest of the code is the same\n",
    "            plt.xlabel('Freq (Hz)')\n",
    "            plt.ylabel('Power')\n",
    "            plt.xlim(right=5)\n",
    "            #make title reflecting peak information\n",
    "            plt.title(vid + \" \" + bodypart + 'Peak value: %f, Location: %f Hz' % (peakY, frqY))\n",
    "\n",
    "            plt.subplot(1,2,2)\n",
    "            plt.plot(x_time, y_normed)\n",
    "            plt.xlabel('Time (s)')\n",
    "            plt.ylabel('Vertical movement')\n",
    "            plt.show()\n",
    "\n",
    "\n",
    "        if savedata:\n",
    "            #save all values to dataframe\n",
    "            vals = [locY]\n",
    "            vals.extend(powers)\n",
    "            vals.append(fps)\n",
    "            vals.append(\"Success\")\n",
    "            print(vals)\n",
    "            results.loc[vid] = vals\n",
    "\n",
    "\n",
    "    except Exception as e:\n",
    "        #emsg = \"\".join(str(e)) #error message string hack :(\n",
    "        emsg = str(e)\n",
    "        #show the error\n",
    "        print(vid, \" Error: \",emsg)\n",
    "        if savedata:\n",
    "            #record error in our results array\n",
    "            vals = failarray.copy()\n",
    "            vals.append(fps)\n",
    "            vals.append(emsg)\n",
    "            print(vals)\n",
    "            results.loc[vid] = vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_end_of_cell_marker": 2,
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "import scipy.fft\n",
    "\n",
    "plotgraphs = True\n",
    "showfirstguess = False\n",
    "savedata = True\n",
    "GuessSMT  = 450  #we expect infants to have a spontaneous ISI of 450 ms\n",
    "GuessFreq = 1000/GuessSMT\n",
    "\n",
    "bodypart = 'handy'\n",
    "\n",
    "failarray = [math.nan, math.nan, math.nan, math.nan,math.nan]  #empty array if we fail to fit function\n",
    "\n",
    "results = pd.DataFrame(columns = (\"peak\", \"400\",\"500\",\"600\",\"700\",\"fps\",\"message\"), index = videos)\n",
    "\n",
    "halfbinwidth = 4\n",
    "targetFreqs = [1000/400, 1000/500, 1000/600, 1000/700]\n",
    "freqResults =pd.DataFrame(columns = (\"ISI400\",\"ISI500\",\"ISI600\",\"ISI700\"))\n",
    "\n",
    "\n",
    "for vid in videos:\n",
    "    #retrieve the raw data for this video\n",
    "    frames = videos[vid]['camera1']['frames']    #how many frames?\n",
    "    fps = videos[vid]['camera1']['fps']          #how many frames per second?\n",
    "    x_data = np.linspace(0,frames,frames+1)      #x axis\n",
    "    x_time = x_data / fps                        #x axis in units of time (seconds)\n",
    "    armmov = avgrh[(vid, people[0], bodypart)]  #get the averaged data\n",
    "    y_data = armmov.truncate(before = 0, after = frames).to_numpy()  #convert to numpy array\n",
    "\n",
    "    #we are just interested in the periodic elements (not absolute value above zero) so substract the mean\n",
    "    y_normed = np.subtract(y_data,np.average(y_data))\n",
    "\n",
    "    try:\n",
    "        yf = scipy.fft.rfft(y_normed)\n",
    "        power = np.abs(yf)**2\n",
    "        xf = scipy.fft.rfftfreq(x_data.size, 1 / fps)\n",
    "\n",
    "        #cutoff at the lower end as these frequencies are not relevant\n",
    "        power[0:10] = 0\n",
    "\n",
    "        #let's find the maximum frequency.\n",
    "        #we hope this is infant's repetitive movement\n",
    "        mY = np.abs(power) # Find magnitude\n",
    "        peakY = np.max(mY) # Find max peak\n",
    "        locY = np.argmax(mY) # Find its location\n",
    "        frqY = xf[locY] # Get the actual frequency value\n",
    "        print(frqY)\n",
    "\n",
    "        #let's find the power in a given frequency bucket corresponding to target freq\n",
    "        powers = [0,0,0,0]\n",
    "        for p in range(4):\n",
    "            #find the nearest value in our list to\n",
    "            bucket = np.argmax(xf>targetFreqs[p])\n",
    "            powers[p] = sum(mY[bucket-halfbinwidth:bucket+halfbinwidth])\n",
    "\n",
    "        print(\"powers\")\n",
    "        print(powers)\n",
    "\n",
    "        if plotgraphs:\n",
    "            plt.figure(figsize=(12, 4),)\n",
    "            plt.subplot(1,2,1)\n",
    "            plt.plot(xf, power)\n",
    "\n",
    "            # New - Plot the max point\n",
    "            plt.plot(frqY, peakY, 'b.', markersize=18)\n",
    "            # Rest of the code is the same\n",
    "            plt.xlabel('Freq (Hz)')\n",
    "            plt.ylabel('Power')\n",
    "            plt.xlim(right=5)\n",
    "            #make title reflecting peak information\n",
    "            plt.title(vid + \" \" + bodypart + 'Peak value: %f, Location: %f Hz' % (peakY, frqY))\n",
    "\n",
    "            plt.subplot(1,2,2)\n",
    "            plt.plot(x_time, y_normed)\n",
    "            plt.xlabel('Time (s)')\n",
    "            plt.ylabel('Vertical movement')\n",
    "            plt.show()\n",
    "\n",
    "\n",
    "        if savedata:\n",
    "            #save all values to dataframe\n",
    "            vals = [locY]\n",
    "            vals.extend(powers)\n",
    "            vals.append(fps)\n",
    "            vals.append(\"Success\")\n",
    "            print(vals)\n",
    "            results.loc[vid] = vals\n",
    "\n",
    "\n",
    "    except Exception as e:\n",
    "        #emsg = \"\".join(str(e)) #error message string hack :(\n",
    "        emsg = str(e)\n",
    "        #show the error\n",
    "        print(vid, \" Error: \",emsg)\n",
    "        if savedata:\n",
    "            #record error in our results array\n",
    "            vals = failarray.copy()\n",
    "            vals.append(fps)\n",
    "            vals.append(emsg)\n",
    "            print(vals)\n",
    "            results.loc[vid] = vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save the fitted parameters.\n",
    "results.to_excel(\"RightHand.FixedBinFreq.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "subsetvids = []\n",
    "count = 0\n",
    "for vid in videos:\n",
    "    subsetvids.append(vid)\n",
    "    count += 1\n",
    "    if count > 12:\n",
    "        break\n",
    "\n",
    "print(subsetvids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbdcbab1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# -\n",
    "#save the fitted parameters.\n",
    "results.to_excel(\"Test.LeftArm.FreeFit.xlsx\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcfb4f19",
   "metadata": {},
   "source": [
    "## Fitting sine waves with fixed frequencies\n",
    "\n",
    "The next thing we try is to try fitting a sine wave with the expected frequency for that trial. We have two conditions A & B.\n",
    "\n",
    "```\n",
    "                      Condition A      Condition B\n",
    "xxxx_04-test-trials      Spontaneous Motor Tempo\n",
    "xxxx_06-test-trials      700ms ISI     400ms ISI\n",
    "xxxx_08-test-trials      500ms ISI     600ms ISI\n",
    "xxxx_10-test-trials      600ms ISI     500ms ISI\n",
    "xxxx_12-test-trials      400ms ISI     700ms ISI\n",
    "xxxx_14-test-trials      Spontaneous Motor Tempo\n",
    "```\n",
    "\n",
    "Therefore, for a given trial type we test the two possible fixed frequency/ISI values to see which fits best.\n",
    "We do not include the Spontaneous Motor Tempo trials in this analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd3bf2aa",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "\n",
    "plotgraphs = True\n",
    "showfirstguess = False\n",
    "savedata = True\n",
    "GuessSMT  = 450  #we expect infants to have a spontaneous ISI of 450 ms\n",
    "GuessFreq = 1000/GuessSMT\n",
    "\n",
    "bodypart = 'rightarmy'\n",
    "\n",
    "failarray = [math.nan, math.nan, math.nan, math.nan,math.nan]  #empty array if we fail to fit function\n",
    "\n",
    "\n",
    "def optimize_ISI(fixeddata, amp, phase, mean):\n",
    "#we need function that we want to fit to our data.\n",
    "#basically a sine wave where we optimize everything - especially frequency.\n",
    "#this is slightly complex because each video could have differnt fps.\n",
    "#so the fixed data is now x_data, fps, ISI\n",
    "#we optimise on phase, amp & mean\n",
    "    isi = fixeddata[\"isi\"]\n",
    "    freq = 1000 / (isi)\n",
    "    x_time = fixeddata[\"x_time\"]\n",
    "    return amp * np.sin(twopi * freq * x_time + phase) + mean\n",
    "\n",
    "\n",
    "results = pd.DataFrame(columns = (\"fit value\",\"best_freq\",\"best_amp\",\"best_phase\",\"best_mean\",\"isi\",\"fps\", \"message\"))\n",
    "\n",
    "for vid in videos:\n",
    "    #retrieve the raw data for this video\n",
    "    frames = videos[vid]['camera1']['frames']    #how many frames?\n",
    "    fps = videos[vid]['camera1']['fps']          #how many frames per second?\n",
    "    x_data = np.linspace(0,frames,frames+1)      #x axis\n",
    "    x_time = x_data / fps                        #x axis in units of time (seconds)\n",
    "    fixeddata = {}\n",
    "    fixeddata[\"fps\"] = fps\n",
    "    fixeddata[\"x_time\"] = x_time\n",
    "    armmov = avgdf[(vid, people[0], bodypart)]  #get the averaged data\n",
    "    y_data = armmov.truncate(before = 0, after = frames).to_numpy()  #convert to numpy array\n",
    "\n",
    "    #WHat ISI\n",
    "    isi = [450, 450]\n",
    "    if '_06-' in vid:\n",
    "        isi[0] = 700\n",
    "        isi[1] = 400\n",
    "    elif '_08-' in vid:\n",
    "        isi[0] = 500\n",
    "        isi[1] = 600\n",
    "    elif '_10-' in vid:\n",
    "        isi[0] = 600\n",
    "        isi[1] = 500\n",
    "    elif '_12-' in vid:\n",
    "        isi[0] = 400\n",
    "        isi[1] = 700\n",
    "    else:\n",
    "        isi[0] = 450\n",
    "        isi[1] = 450\n",
    "\n",
    "    for iisi in isi:\n",
    "        try:\n",
    "            #have a first guess of the parameters\n",
    "            guess_mean = float(np.mean(armmov))\n",
    "            guess_amp =  2* float(np.std(armmov))\n",
    "            guess_phase = 0.0\n",
    "            guess_freq = 1000 / iisi  # freq in Hz\n",
    "\n",
    "            fixeddata[\"isi\"] = iisi\n",
    "\n",
    "            # we'll use this to plot our first estimate. This might already be good enough for you\n",
    "            data_first_guess = optimize_ISI(fixeddata, guess_amp, guess_phase, guess_mean)\n",
    "\n",
    "            #now optimize\n",
    "            params, params_covariance = optimize.curve_fit(optimize_ISI, fixeddata, y_data,\n",
    "                                                       p0=[guess_amp,guess_phase,guess_mean])\n",
    "\n",
    "            #what is the best fit?\n",
    "            fitted = optimize_ISI(fixeddata, params[0],params[1],params[2])\n",
    "\n",
    "            #how far is this from raw data?\n",
    "            best_fit = scipy.spatial.distance.euclidean(y_data,fitted )\n",
    "\n",
    "            if savedata:\n",
    "                #save all values to dataframe\n",
    "                vals = [best_fit]\n",
    "                vals.append(guess_freq)\n",
    "                vals.extend(params)\n",
    "                vals.append(iisi)\n",
    "                vals.append(fps)\n",
    "                vals.append(\"Success\")\n",
    "                print(vals)\n",
    "                results.loc[vid+ \"-\" + str(iisi)] = vals\n",
    "\n",
    "            if plotgraphs:\n",
    "                #plot the functions\n",
    "                plt.figure(figsize=(6, 4),)\n",
    "                plt.plot(x_time, y_data, label='Data')\n",
    "                plt.plot(x_time, fitted,label='Fitted function')\n",
    "                if showfirstguess:\n",
    "                    plt.plot(x_time, data_first_guess, label = 'First guess')\n",
    "                plt.title(vid + \" \" + str(iisi) + \"ms \" + bodypart )\n",
    "                plt.legend(loc='best')\n",
    "\n",
    "                plt.show()\n",
    "\n",
    "        except Exception as e:\n",
    "            #emsg = \"\".join(str(e)) #error message string hack :(\n",
    "            emsg = str(e)\n",
    "            #show the error\n",
    "            print(vid, \" Error: \",emsg)\n",
    "            if savedata:\n",
    "                #record error in our results array\n",
    "                vals = failarray.copy()\n",
    "                vals.append(iisi)\n",
    "                vals.append(fps)\n",
    "                vals.append(emsg)\n",
    "                print(vals)\n",
    "                results.loc[vid+ \"-\" + str(iisi)] = vals\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d5b0f80",
   "metadata": {},
   "outputs": [],
   "source": [
    "#save the fitted parameters.\n",
    "results.to_excel(\"RightArm.FixedFit.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "924aa2ed",
   "metadata": {},
   "source": [
    "## 3.4 Finding fundamental frequency with FFT\n",
    "\n",
    "\n",
    "A good guide can be found here https://realpython.com/python-scipy-fft/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "import scipy.fft\n",
    "\n",
    "plotgraphs = True\n",
    "showfirstguess = False\n",
    "savedata = True\n",
    "GuessSMT  = 450  #we expect infants to have a spontaneous ISI of 450 ms\n",
    "GuessFreq = 1000/GuessSMT\n",
    "\n",
    "bodypart = 'rightarmy'\n",
    "\n",
    "failarray = [math.nan, math.nan, math.nan, math.nan,math.nan]  #empty array if we fail to fit function\n",
    "\n",
    "results = pd.DataFrame(columns = (\"peak\", \"400\",\"500\",\"600\",\"700\",\"fps\",\"message\"), index = videos)\n",
    "\n",
    "halfbinwidth = 4\n",
    "targetFreqs = [1000/400, 1000/500, 1000/600, 1000/700]\n",
    "freqResults =pd.DataFrame(columns = (\"ISI400\",\"ISI500\",\"ISI600\",\"ISI700\"))\n",
    "\n",
    "\n",
    "for vid in videos:\n",
    "    #retrieve the raw data for this video\n",
    "    frames = videos[vid]['camera1']['frames']    #how many frames?\n",
    "    fps = videos[vid]['camera1']['fps']          #how many frames per second?\n",
    "    x_data = np.linspace(0,frames,frames+1)      #x axis\n",
    "    x_time = x_data / fps                        #x axis in units of time (seconds)\n",
    "    armmov = avgdf[(vid, people[0], bodypart)]  #get the averaged data\n",
    "    y_data = armmov.truncate(before = 0, after = frames).to_numpy()  #convert to numpy array\n",
    "\n",
    "    #we are just interested in the periodic elements (not absolute value above zero) so substract the mean\n",
    "    y_normed = np.subtract(y_data,np.average(y_data))\n",
    "\n",
    "    try:\n",
    "        yf = scipy.fft.rfft(y_normed)\n",
    "        power = np.abs(yf)**2\n",
    "        xf = scipy.fft.rfftfreq(x_data.size, 1 / fps)\n",
    "\n",
    "        #let's find the maximum frequency.\n",
    "        #we hope this is infant's repetitive movement\n",
    "        mY = np.abs(power) # Find magnitude\n",
    "        peakY = np.max(mY) # Find max peak\n",
    "        locY = np.argmax(mY) # Find its location\n",
    "        frqY = xf[locY] # Get the actual frequency value\n",
    "        print(frqY)\n",
    "\n",
    "        #let's find the power in a given frequency bucket corresponding to target freq\n",
    "        powers = [0,0,0,0]\n",
    "        for p in range(4):\n",
    "            #find the nearest value in our list to\n",
    "            bucket = np.argmax(xf>targetFreqs[p])\n",
    "            powers[p] = sum(mY[bucket-halfbinwidth:bucket+halfbinwidth])\n",
    "\n",
    "        print(\"powers\")\n",
    "        print(powers)\n",
    "\n",
    "        if plotgraphs:\n",
    "            plt.figure(figsize=(12, 4),)\n",
    "            plt.subplot(1,2,1)\n",
    "            plt.plot(xf, power)\n",
    "\n",
    "            # New - Plot the max point\n",
    "            plt.plot(frqY, peakY, 'b.', markersize=18)\n",
    "            # Rest of the code is the same\n",
    "            plt.xlabel('Freq (Hz)')\n",
    "            plt.ylabel('Power')\n",
    "            plt.xlim(right=5)\n",
    "            #make title reflecting peak information\n",
    "            plt.title(vid + \" \" + bodypart + 'Peak value: %f, Location: %f Hz' % (peakY, frqY))\n",
    "\n",
    "            plt.subplot(1,2,2)\n",
    "            plt.plot(x_time, y_normed)\n",
    "            plt.xlabel('Time (s)')\n",
    "            plt.ylabel('Vertical movement')\n",
    "            plt.show()\n",
    "\n",
    "\n",
    "        if savedata:\n",
    "            #save all values to dataframe\n",
    "            vals = [locY]\n",
    "            vals.extend(powers)\n",
    "            vals.append(fps)\n",
    "            vals.append(\"Success\")\n",
    "            print(vals)\n",
    "            results.loc[vid] = vals\n",
    "\n",
    "\n",
    "    except Exception as e:\n",
    "        #emsg = \"\".join(str(e)) #error message string hack :(\n",
    "        emsg = str(e)\n",
    "        #show the error\n",
    "        print(vid, \" Error: \",emsg)\n",
    "        if savedata:\n",
    "            #record error in our results array\n",
    "            vals = failarray.copy()\n",
    "            vals.append(fps)\n",
    "            vals.append(emsg)\n",
    "            print(vals)\n",
    "            results.loc[vid] = vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "965cb689",
   "metadata": {
    "lines_to_end_of_cell_marker": 2,
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "import scipy.fft\n",
    "\n",
    "plotgraphs = True\n",
    "showfirstguess = False\n",
    "savedata = True\n",
    "GuessSMT  = 450  #we expect infants to have a spontaneous ISI of 450 ms\n",
    "GuessFreq = 1000/GuessSMT\n",
    "\n",
    "bodypart = 'handy'\n",
    "\n",
    "failarray = [math.nan, math.nan, math.nan, math.nan,math.nan]  #empty array if we fail to fit function\n",
    "\n",
    "results = pd.DataFrame(columns = (\"peak\", \"400\",\"500\",\"600\",\"700\",\"fps\",\"message\"), index = videos)\n",
    "\n",
    "halfbinwidth = 4\n",
    "targetFreqs = [1000/400, 1000/500, 1000/600, 1000/700]\n",
    "freqResults =pd.DataFrame(columns = (\"ISI400\",\"ISI500\",\"ISI600\",\"ISI700\"))\n",
    "\n",
    "\n",
    "for vid in videos:\n",
    "    #retrieve the raw data for this video\n",
    "    frames = videos[vid]['camera1']['frames']    #how many frames?\n",
    "    fps = videos[vid]['camera1']['fps']          #how many frames per second?\n",
    "    x_data = np.linspace(0,frames,frames+1)      #x axis\n",
    "    x_time = x_data / fps                        #x axis in units of time (seconds)\n",
    "    armmov = avgrh[(vid, people[0], bodypart)]  #get the averaged data\n",
    "    y_data = armmov.truncate(before = 0, after = frames).to_numpy()  #convert to numpy array\n",
    "\n",
    "    #we are just interested in the periodic elements (not absolute value above zero) so substract the mean\n",
    "    y_normed = np.subtract(y_data,np.average(y_data))\n",
    "\n",
    "    try:\n",
    "        yf = scipy.fft.rfft(y_normed)\n",
    "        power = np.abs(yf)**2\n",
    "        xf = scipy.fft.rfftfreq(x_data.size, 1 / fps)\n",
    "\n",
    "        #cutoff at the lower end as these frequencies are not relevant\n",
    "        power[0:10] = 0\n",
    "\n",
    "        #let's find the maximum frequency.\n",
    "        #we hope this is infant's repetitive movement\n",
    "        mY = np.abs(power) # Find magnitude\n",
    "        peakY = np.max(mY) # Find max peak\n",
    "        locY = np.argmax(mY) # Find its location\n",
    "        frqY = xf[locY] # Get the actual frequency value\n",
    "        print(frqY)\n",
    "\n",
    "        #let's find the power in a given frequency bucket corresponding to target freq\n",
    "        powers = [0,0,0,0]\n",
    "        for p in range(4):\n",
    "            #find the nearest value in our list to\n",
    "            bucket = np.argmax(xf>targetFreqs[p])\n",
    "            powers[p] = sum(mY[bucket-halfbinwidth:bucket+halfbinwidth])\n",
    "\n",
    "        print(\"powers\")\n",
    "        print(powers)\n",
    "\n",
    "        if plotgraphs:\n",
    "            plt.figure(figsize=(12, 4),)\n",
    "            plt.subplot(1,2,1)\n",
    "            plt.plot(xf, power)\n",
    "\n",
    "            # New - Plot the max point\n",
    "            plt.plot(frqY, peakY, 'b.', markersize=18)\n",
    "            # Rest of the code is the same\n",
    "            plt.xlabel('Freq (Hz)')\n",
    "            plt.ylabel('Power')\n",
    "            plt.xlim(right=5)\n",
    "            #make title reflecting peak information\n",
    "            plt.title(vid + \" \" + bodypart + 'Peak value: %f, Location: %f Hz' % (peakY, frqY))\n",
    "\n",
    "            plt.subplot(1,2,2)\n",
    "            plt.plot(x_time, y_normed)\n",
    "            plt.xlabel('Time (s)')\n",
    "            plt.ylabel('Vertical movement')\n",
    "            plt.show()\n",
    "\n",
    "\n",
    "        if savedata:\n",
    "            #save all values to dataframe\n",
    "            vals = [locY]\n",
    "            vals.extend(powers)\n",
    "            vals.append(fps)\n",
    "            vals.append(\"Success\")\n",
    "            print(vals)\n",
    "            results.loc[vid] = vals\n",
    "\n",
    "\n",
    "    except Exception as e:\n",
    "        #emsg = \"\".join(str(e)) #error message string hack :(\n",
    "        emsg = str(e)\n",
    "        #show the error\n",
    "        print(vid, \" Error: \",emsg)\n",
    "        if savedata:\n",
    "            #record error in our results array\n",
    "            vals = failarray.copy()\n",
    "            vals.append(fps)\n",
    "            vals.append(emsg)\n",
    "            print(vals)\n",
    "            results.loc[vid] = vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "766bbb9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#save the fitted parameters.\n",
    "results.to_excel(\"RightHand.FixedBinFreq.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "095ea178",
   "metadata": {},
   "outputs": [],
   "source": [
    "subsetvids = []\n",
    "count = 0\n",
    "for vid in videos:\n",
    "    subsetvids.append(vid)\n",
    "    count += 1\n",
    "    if count > 12:\n",
    "        break\n",
    "\n",
    "print(subsetvids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4 Comparing to human coding.\n",
    "\n",
    "We have a spreadsheet of syhnchrony scores for each parent infant dyad. Here we see if we can find a measure that correlates with the human scores.\n",
    "\n",
    "First, load up the spreadsheet.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "excelpath = projectpath + \"\\\\SS_CARE.xlsx\"\n",
    "\n",
    "filename, file_format = os.path.splitext(excelpath)\n",
    "if file_format and file_format == 'xls':\n",
    "    # use default reader\n",
    "    videolist = pd.read_excel(excelpath)\n",
    "else:\n",
    "    #since dec 2020 read_excel no longer supports xlsx (!?) so need to use openpyxl like so..\n",
    "    videolist = pd.read_excel(excelpath, engine = \"openpyxl\")\n",
    "\n",
    "videolist = videolist.set_index(\"subject\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#take a quick look\n",
    "videolist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#copy the dyad syncrhony and maternal sensitivity scores into our data frame.\n",
    "results[\"DyadSynScore\"] = videolist[\"DyadSyn\"]\n",
    "results[\"MatSensScore\"] = videolist[\"MatSens\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#take a quick look\n",
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#scatter plots of these results.\n",
    "plt.scatter(results[\"DyadSynScore\"],results[\"corrArms\"], )\n",
    "plt.title(\"Correlation between expert rated synchrony and time series correlations\")\n",
    "plt.xlabel(\"Dyad Synchroncy Score\")\n",
    "plt.ylabel(\"Dyad Correlation\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rolling_r.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "d1 = vardf[(vid, who[0], parts[0])]\n",
    "d2 = vardf[(vid, who[1], parts[0])]\n",
    "seconds = 5\n",
    "fps = 25\n",
    "wholeads = who[0] + 'leads <> ' + who[1] + ' leads'\n",
    "rs = [vasc.crosscorr(d1,d2, lag) for lag in range(-int(seconds*fps-1),int(seconds*fps))]\n",
    "offset = np.ceil(len(rs)/2)-np.argmax(rs)\n",
    "f,ax=plt.subplots(figsize=(14,3))\n",
    "ax.plot(rs)\n",
    "ax.axvline(np.ceil(len(rs)/2),color='k',linestyle='--',label='Center')\n",
    "ax.axvline(np.argmax(rs),color='r',linestyle='--',label='Peak synchrony')\n",
    "ax.set(title=f'Offset = {offset} frames\\n' + wholeads,ylim=[.0,1],xlim=[0,300], xlabel='Offset',ylabel='Pearson r')\n",
    "ax.set_xticklabels([int(item-150) for item in ax.get_xticks()]);\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4 Granger Causality\n",
    "\n",
    "The next thing to look at is if the movements of the infant predict the movements of the parent. This would suggest parent is responding to the infant.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "\n",
    "https://towardsdatascience.com/granger-causality-and-vector-auto-regressive-model-for-time-series-forecasting-3226a64889a6\n",
    "\n",
    "https://www.machinelearningplus.com/time-series/time-series-analysis-python/\n",
    "\n",
    "https://towardsdatascience.com/four-ways-to-quantify-synchrony-between-time-series-data-b99136c4a9c9\n"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,py:light"
  },
  "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}