Step3b.AnalyseData.scipy.synchrony.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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "settingsjson = \"C:\\\\Users\\\\cas\\\\OneDrive - Goldsmiths College\\\\Projects\\\\Little Drummers\\\\VASC\\\\settings.json\"\n",
    "settingsjson = \"E:\\\\little.drummer.20220111\\\\LD.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. (You might be missing a comma or a bracket.)\")\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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "includeHands = settings[\"flags\"][\"includeHands\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anon = settings[\"flags\"][\"anon\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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",
    "print(videos_in)\n",
    "print(videos_out)\n",
    "print(videos_out_openpose)\n",
    "print(videos_out_timeseries)\n",
    "print(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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('reading parquet file(s):')\n",
    "df = pq.read_table(videos_out_timeseries + '\\\\' + settings[\"filenames\"][\"cleandataparquet\"]).to_pandas()\n",
    "\n",
    "#sort the column names as this helps with indexing\n",
    "df = df.sort_index(axis = 1)\n",
    "\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",
    "    lh = lh.sort_index(axis = 1)\n",
    "    rh = rh.sort_index(axis = 1)\n",
    "    \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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#are we going to use all the data or a subset?\n",
    "first = 0\n",
    "last = df.shape[0]\n",
    "\n",
    "df = df.truncate(before  = first, after = last)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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 (per frame) 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": null,
   "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",
    "\n",
    "handpoints = {\"hand\" : vasc.hxys,\n",
    "       \"handx\" : vasc.hxs,\n",
    "       \"handy\" : vasc.hys}\n",
    "\n",
    "avglh = vasc.averageCoordinateTimeSeries(lh,handpoints,vids,people)\n",
    "avgrh = vasc.averageCoordinateTimeSeries(rh,handpoints,vids,people)\n"
   ]
  },
  {
   "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 Combining hand and wrist data\n",
    "\n",
    "If we have hand data from OpenPose then we can combine this with the body data to get more accurate movements for the left and right hands. Similar to what we did on the step before but combining across multiple dataframes.\n",
    "\n",
    "TODO: We may also use a weighted sum of the points from elbow, wrist and hand.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#to calculate the mean, we need to let routine know what points to use. \n",
    "#From the arm this the points we want to combine\n",
    "leftarmpoints = {\"hand\" : vasc.leftarmxys,\n",
    "               \"handx\": vasc.leftarmx,\n",
    "               \"handy\": vasc.leftarmy}\n",
    "\n",
    "rightarmpoints = {\"hand\" : vasc.rightarmxys,\n",
    "               \"handx\": vasc.rightarmx,\n",
    "               \"handy\": vasc.rightarmy}\n",
    "\n",
    "leftwristpoints = {\"hand\" : vasc.leftwristxys,\n",
    "               \"handx\": vasc.leftwristx,\n",
    "               \"handy\": vasc.leftwristy}\n",
    "\n",
    "rightwristpoints = {\"hand\" : vasc.rightwristxys,\n",
    "               \"handx\": vasc.rightwristx,\n",
    "               \"handy\": vasc.rightwristy}\n",
    "\n",
    "handpoints = {\"hand\" : vasc.hxys,\n",
    "       \"handx\" : vasc.hxs,\n",
    "       \"handy\" : vasc.hys}\n",
    "\n",
    "\n",
    "rightarmhand = vasc.averageArmHandTimeSeries(df,rh,rightarmpoints,handpoints,vids,people)\n",
    "leftarmhand  = vasc.averageArmHandTimeSeries(df,lh,leftarmpoints,handpoints,vids,people)\n",
    "\n",
    "#combine hand and wrist data weighting in ratio\n",
    "wristtohandweightratio = 21\n",
    "\n",
    "rightwristhand = vasc.averageArmHandTimeSeries(df,rh,rightwristpoints,handpoints,vids,people,wristtohandweightratio)\n",
    "leftwristhand  = vasc.averageArmHandTimeSeries(df,lh,leftwristpoints,handpoints,vids,people, wristtohandweightratio)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rightarmhand"
   ]
  },
  {
   "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": null,
   "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": null,
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "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",
    "series = [\"Body\"]\n",
    "if includeHands:\n",
    "    series = [\"Body\",\"Left Hand\",\"Right Hand\",\"Left Arm & Hand\",\"Right Arm & Hand\",\"Left Wrist & Hand (Weighted)\",\"Right Wrist & Hand (Weighted)\"]\n",
    "\n",
    "pickseries = widgets.Dropdown(\n",
    "    options= series,\n",
    "    value= series[0],\n",
    "    description='Series:'\n",
    ")\n",
    "\n",
    "    \n",
    "mainfeatures = []\n",
    "for f in meanpoints:\n",
    "    mainfeatures.append(f)\n",
    "\n",
    "handfeatures = []\n",
    "for f in handpoints:\n",
    "    handfeatures.append(f)\n",
    "    \n",
    "\n",
    "pickfeature = widgets.Dropdown(\n",
    "    options= mainfeatures,\n",
    "    value= mainfeatures[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 pickseries_change(change):\n",
    "    if change['name'] == 'value' and (change['new'] != change['old']):\n",
    "        if pickseries.value == series[0]:\n",
    "            pickfeature.options = mainfeatures\n",
    "        else:\n",
    "            pickfeature.options = handfeatures\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",
    "pickseries.observe(pickseries_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",
    "txt = widgets.Label(value=\"Which arm and start/end points do we use?\")\n",
    "leftright = widgets.RadioButtons(\n",
    "    options=['left', 'right'],\n",
    "    value='right', # Defaults to 'pineapple'\n",
    "    description='Arm',\n",
    "    disabled=False\n",
    ")\n",
    "def leftright_change(change):\n",
    "    if change['name'] == 'value' and (change['new'] != change['old']):\n",
    "        videos[pickvid.value][\"camera1\"][\"side\"] = leftright.value\n",
    "\n",
    "leftright.observe(leftright_change,'value')\n",
    "\n",
    "## a couple of sliders to set the start and end points.\n",
    "startslider = widgets.IntSlider(\n",
    "    value=0,\n",
    "    min=0,\n",
    "    max=videos[pickvid.value][\"camera1\"][\"end\"] ,\n",
    "    step=1,\n",
    "    description='Start Frame:',\n",
    "    continuous_update=False,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='d',\n",
    "    layout=widgets.Layout(width='800px')\n",
    ")\n",
    "endslider = widgets.IntSlider(\n",
    "    value=161,\n",
    "    min=0,\n",
    "    max=161,\n",
    "    step=1,\n",
    "    description='End Frame:',\n",
    "    continuous_update=False,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='d',\n",
    "    layout=widgets.Layout(width='800px')\n",
    ")\n",
    "\n",
    "def slider_change(slider):\n",
    "    videos[pickvid.value][\"camera1\"][\"start\"] = startslider.value\n",
    "    videos[pickvid.value][\"camera1\"][\"end\"] = max(startslider.value,endslider.value) #can't go below start.\n",
    "    #logging.info (videos[vid][\"camera1\"])\n",
    "    updateAll(False)\n",
    "\n",
    "startslider.observe(slider_change, 'value')\n",
    "endslider.observe(slider_change, 'value')\n",
    "    \n",
    "def drawGraphs(vid, seriestype, feature, linetype):\n",
    "    \"\"\"Plot input signals\"\"\"\n",
    "    plt.ion()\n",
    "\n",
    "#    f,ax=plt.subplots(4,1,figsize=(14,10),sharex=True)\n",
    "    f,ax=plt.subplots(2,1,figsize=(14,6),sharex=False)\n",
    "    ax[0].set_title('Infant')\n",
    "    ax[0].set_xlabel('Frames')\n",
    "#    ax[1].set_title('Parent')\n",
    "    ax[1].set_xlabel('Seconds')\n",
    "\n",
    "    who = [\"infant\",\"parent\"]\n",
    "\n",
    "    #what time series (smoothed, etc) are we plotting?\n",
    "    if linetype == linetypes[0]:\n",
    "        usedf = avgdf\n",
    "    elif linetype == linetypes[1]:\n",
    "        usedf = smoothdf\n",
    "    else:\n",
    "        usedf = vardf\n",
    "    \n",
    "    #now pick the correct dataseries\n",
    "    if seriestype == series[0]:\n",
    "        #nothing to change\n",
    "        usedf = usedf\n",
    "    elif seriestype == series[1]:\n",
    "        usedf = avglh\n",
    "    elif seriestype == series[2]:\n",
    "        usedf = avgrh\n",
    "    elif seriestype == series[3]:\n",
    "        usedf = lefttarmhand\n",
    "    elif seriestype == series[4]:\n",
    "        usedf = rightarmhand\n",
    "    elif seriestype == series[5]:\n",
    "        usedf = leftwristhand\n",
    "    elif seriestype == series[6]:\n",
    "        usedf = rightwristhand    \n",
    "    else:\n",
    "        raise ValueError(\"Unknown series type selected\")\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",
    "    fps = videos[vid]['camera1']['fps']  \n",
    "    x_time = n / fps #scale by frame per second to get real time\n",
    "    ax[0].plot(infant)\n",
    "    #add vertical lines for the start and end points for analysis\n",
    "    #todo make vertical lines work as times\n",
    "    starttime = videos[vid][\"camera1\"][\"start\"] \n",
    "    endtime = videos[vid][\"camera1\"][\"end\"]  \n",
    "    #logging.info(videos[vid][\"camera1\"])\n",
    "    ax[0].axvline(x=starttime,c='tab:green')\n",
    "    ax[0].axvline(x=endtime,c='tab:red')\n",
    "    ax[1].plot(x_time, infant /fps)\n",
    "#    ax[1].plot(parent, color='b')\n",
    "    \n",
    "#    ax[2].plot(usedf.loc[:,(vid, slice(None), feature)])\n",
    "#    ax[2].set(xlabel='Time',ylabel='Movement index for parent and infant')\n",
    "\n",
    "#    if seriestype == series[0]:\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",
    "#        rolling_r.plot(ax=ax[3])\n",
    "#        ax[3].set(xlabel='Time (seconds)',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",
    "        fps = videos[pickvid.value]['camera1']['fps']  \n",
    "        startslider.max =  videos[pickvid.value][\"camera1\"][\"frames\"]\n",
    "        endslider.max =  videos[pickvid.value][\"camera1\"][\"frames\"]  \n",
    "        leftright.value = videos[pickvid.value][\"camera1\"][\"side\"] \n",
    "        logging.debug('forceUpdate')\n",
    "    with output:\n",
    "        display(pickvid,pickseries,pickfeature,picktype,txt,leftright,startslider,endslider,button_update)\n",
    "        drawGraphs(pickvid.value,pickseries.value,pickfeature.value,picktype.value)\n",
    "\n",
    "#draw everything for first time\n",
    "updateAll(True)\n",
    "output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting at start time across all videos\n",
    "\n",
    "The next cell loops through all entries in `clean.json` and sets the `start` time the same for all cases. \n",
    "Remember that video data is frame by frame and each video could have a different frame rate (`fps`). So we\n",
    "need to convert time into frames.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "starttime = 0 \n",
    "starttime = 5  #(in seconds)\n",
    "\n",
    "for vid in videos:\n",
    "    for cam in videos[vid]:\n",
    "        videos[vid][cam][\"side\"] = \"right\"\n",
    "        videos[vid][cam][\"start\"] = int(starttime * videos[vid][cam][\"fps\"]) #convert time to number of frames\n",
    "        videos[vid][cam][\"end\"] = videos[vid][cam][\"frames\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "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": [
    "who = [\"infant\",\"parent\"]\n",
    "parts = [\"head\",\"arms\",\"all\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "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": null,
   "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(videos_out + \"//Test.LeftArm.FreeFit.xlsx\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "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(videos_out + \"//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/\n",
    "\n",
    "Filtering explained here https://scipy-lectures.org/intro/scipy/auto_examples/plot_fftpack.html\n",
    "\n",
    "\n",
    "In this first block of code we want to see the freequency power at each target ISI/frequency and the absolute peak "
   ]
  },
  {
   "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",
    "#whichdf = avgdf\n",
    "#bodypart = 'rightarmy'\n",
    "\n",
    "whichdf = leftwristhand\n",
    "bodypart = 'handy'\n",
    "\n",
    "failarray = [math.nan, 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\",\"450\", \"400\",\"500\",\"600\",\"700\",\"fps\",\"message\"), index = videos)\n",
    "\n",
    "halfbinwidth = 4\n",
    "targetFreqs = [1000/450, 1000/400, 1000/500, 1000/600, 1000/700]\n",
    "#freqResults =pd.DataFrame(columns = (\"SMT450\",\"ISI400\",\"ISI500\",\"ISI600\",\"ISI700\"))\n",
    "\n",
    "#what is the lower end cutoff in Hertz? \n",
    "#note - since each video has a different speed (diff fps) we need to convert/this value for each vid.\n",
    "lowendcutoffinHz = .5\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 = whichdf[(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",
    "        yfft = scipy.fft.rfft(y_normed)\n",
    "        power = np.abs(yfft)**2\n",
    "        xfreq = scipy.fft.rfftfreq(x_data.size, 1 / fps)\n",
    "\n",
    "        #frequency cutoff at the lower end \n",
    "        #as these frequencies are not relevant  \n",
    "        power[np.abs(xfreq) < lowendcutoffinHz] = 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 = xfreq[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,0]\n",
    "        for p in range(5):\n",
    "            #find the nearest value in our list to\n",
    "            bucket = np.argmax(xfreq>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",
    "\n",
    "            #plot the vertical movement of this hand\n",
    "            plt.subplot(1,2,1)\n",
    "            plt.plot(x_time, y_normed)\n",
    "            plt.xlabel('Time (s)')\n",
    "            plt.ylabel('Vertical movement')\n",
    "            \n",
    "            #now plot the Fourier power \n",
    "            plt.subplot(1,2,2)\n",
    "            plt.plot(xfreq, power)\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.show()\n",
    "\n",
    "\n",
    "        if savedata:\n",
    "            #save all values to dataframe\n",
    "            vals = []\n",
    "            vals.append(frqY)  #max frequency\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(\"LeftHand.FixedBinFreq.xlsx\")"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "for vid in videos:\n",
    "    for cam in videos[vid]:\n",
    "        videos[vid][cam][\"side\"] = \"right\"\n",
    "        videos[vid][cam][\"start\"] = 0\n",
    "        videos[vid][cam][\"end\"] = videos[vid][cam][\"frames\"]"
   ]
  },
  {
   "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"
   ]
  }
 ],
 "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
}