reproduce_tables_2_3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "81ae39b3",
   "metadata": {},
   "source": [
    "# Code to reproduce Tables 2 and 3 from the paper 'Improving cognitive-state analysis from eye gaze with synthetic eye-movement data'\n",
    "* see: https://www.sciencedirect.com/science/article/pii/S0097849324000281"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7546f270",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import random\n",
    "from tqdm.notebook import tqdm\n",
    "import sys\n",
    "import joblib\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "from typing import Dict\n",
    "from typing import List\n",
    "from typing import Tuple\n",
    "\n",
    "\n",
    "# tensorflow and SP-EyeGAN\n",
    "import tensorflow\n",
    "import tensorflow as tf\n",
    "from sp_eyegan.model import eventGAN as eventGAN\n",
    "from sp_eyegan.model import vae_baseline as vae\n",
    "from scipy.stats import ttest_ind, ttest_1samp\n",
    "\n",
    "%pylab inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d24113b",
   "metadata": {},
   "source": [
    "### Helper Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a17825c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def vel_to_dva(vel_data, x_start = 0,\n",
    "             y_start = 0):\n",
    "    x_vel = vel_data[:,0]\n",
    "    y_vel = vel_data[:,1]\n",
    "    x_px  = []\n",
    "    y_px  = []\n",
    "    cur_x_pos = x_start\n",
    "    cur_y_pos = y_start\n",
    "    for i in range(len(x_vel)):\n",
    "        x_px.append(cur_x_pos + x_vel[i])\n",
    "        y_px.append(cur_y_pos + y_vel[i])\n",
    "        cur_x_pos = x_px[-1]\n",
    "        cur_y_pos = y_px[-1]\n",
    "    return np.concatenate([np.expand_dims(np.array(x_px),axis=1),\n",
    "                           np.expand_dims(np.array(y_px),axis=1)],axis=1)\n",
    "\n",
    "def get_fixation_stats(input_data):\n",
    "    real_x_vels      = []\n",
    "    real_y_vels      = []\n",
    "    real_vels        = []\n",
    "    real_dispersions = []\n",
    "    cur_vel_data     = input_data\n",
    "    for i in tqdm(np.arange(cur_vel_data.shape[0])):\n",
    "        cur_vels = cur_vel_data[i]\n",
    "        cur_dva = vel_to_dva(cur_vels)\n",
    "        try:\n",
    "            end_id = np.where(np.logical_and(cur_vels[:,0] == 0,\n",
    "                                             cur_vels[:,1] ==0 ))[0][0]\n",
    "        except:\n",
    "            end_id = len(cur_vels)\n",
    "        if end_id == 0:\n",
    "            continue\n",
    "        cur_vels = cur_vels[0:end_id]\n",
    "        cur_dva = cur_dva[0:end_id]\n",
    "        x_dva = cur_dva[:,0]\n",
    "        y_dva = cur_dva[:,1]\n",
    "        real_x_vels.append(cur_vels[:,0])\n",
    "        real_y_vels.append(cur_vels[:,1])\n",
    "        \n",
    "        vels =  np.power(np.array(real_x_vels[-1]),2) +\\\n",
    "                np.power(np.array(real_y_vels[-1]),2)\n",
    "        vels = np.sqrt(vels)\n",
    "        real_vels.append(vels)\n",
    "        x_amp = np.abs(np.max(x_dva) - np.min(x_dva))\n",
    "        y_amp = np.abs(np.max(y_dva) - np.min(y_dva))\n",
    "        cur_dispersion = x_amp + y_amp\n",
    "        real_dispersions.append(cur_dispersion)\n",
    "    \n",
    "    return real_vels, real_x_vels,real_y_vels, real_dispersions\n",
    "\n",
    "\n",
    "def get_saccade_stats(input_data, max_velocity = 0.5):\n",
    "    real_x_vels      = []\n",
    "    real_y_vels      = []\n",
    "    real_vels        = []\n",
    "    real_amplitudes  = []\n",
    "    real_x_accs      = []\n",
    "    real_y_accs      = []\n",
    "    real_accs        = []\n",
    "    cur_vel_data     = input_data\n",
    "    for i in tqdm(np.arange(cur_vel_data.shape[0])):\n",
    "        cur_vels = cur_vel_data[i]\n",
    "        cur_dva = vel_to_dva(cur_vels)\n",
    "        try:\n",
    "            end_id = np.where(np.logical_and(cur_vels[:,0] == 0,\n",
    "                                             cur_vels[:,1] ==0 ))[0][0]\n",
    "        except:\n",
    "            end_id = len(cur_vels)\n",
    "        if end_id == 0:\n",
    "            continue\n",
    "        cur_vels = cur_vels[0:end_id]\n",
    "        cur_dva = cur_dva[0:end_id]\n",
    "        x_dva = cur_dva[:,0]\n",
    "        y_dva = cur_dva[:,1]\n",
    "        real_x_vels.append(cur_vels[:,0])\n",
    "        real_y_vels.append(cur_vels[:,1])\n",
    "        \n",
    "        vels =  np.power(np.array(real_x_vels[-1]),2) +\\\n",
    "                np.power(np.array(real_y_vels[-1]),2)\n",
    "        vels = np.sqrt(vels)\n",
    "        real_vels.append(vels)        \n",
    "        \n",
    "        x_accs = cur_vels[1:,0] - cur_vels[:-1,0]\n",
    "        y_accs = cur_vels[1:,1] - cur_vels[:-1,0]\n",
    "        \n",
    "        accs = np.power(np.array(x_accs),2) +\\\n",
    "                np.power(np.array(y_accs),2)\n",
    "        accs = np.sqrt(accs)\n",
    "        \n",
    "        real_x_accs.append(x_accs)\n",
    "        real_y_accs.append(y_accs)\n",
    "        real_accs.append(accs)\n",
    "        cur_complete_vels = np.sqrt(np.power(cur_vels[:,0],2) + np.power(cur_vels[:,1],2))\n",
    "        cur_complete_vels[cur_complete_vels > max_velocity] = max_velocity\n",
    "        cur_amplitude = np.sqrt(np.power(x_dva[0] - x_dva[-1],2) + np.power(y_dva[0] - y_dva[-1],2))\n",
    "        real_amplitudes.append(cur_amplitude)\n",
    "    \n",
    "    return real_vels, real_x_vels, real_y_vels, real_accs, real_x_accs, real_y_accs, np.array(real_amplitudes)\n",
    "\n",
    "\n",
    "# calculate the kl divergence\n",
    "def kl_divergence(p, q):\n",
    "    from math import log2\n",
    "    from math import sqrt\n",
    "    rel_etropies = [p[i] * log2(p[i]/q[i]) for i in range(len(p))]\n",
    "    return np.sum(np.array(rel_etropies,dtype=np.float32))\n",
    " \n",
    "# calculate the js divergence\n",
    "def js_divergence(p, q):\n",
    "    m = 0.5 * (p + q)\n",
    "    return 0.5 * kl_divergence(p, m) + 0.5 * kl_divergence(q, m)\n",
    "\n",
    "\n",
    "\n",
    "def js_divergence_sampling( values_1, values_2,\n",
    "                            bins,\n",
    "                            epsilon =  0.00001,\n",
    "                            iterations = 10,\n",
    "                            number_per_iter = 10000,\n",
    "                            random_state = 42):\n",
    "    np.random.seed(random_state)\n",
    "    values = []\n",
    "    for iter in range(iterations):\n",
    "        idx_1 = np.random.choice(len(values_1), number_per_iter)\n",
    "        idx_2 = np.random.choice(len(values_2), number_per_iter)\n",
    "        counts_1,_ = np.histogram(values_1[idx_1], bins = bins, density=True)\n",
    "        counts_1 /= np.sum(counts_1)\n",
    "        counts_2,_ = np.histogram(values_2[idx_2], bins = bins, density=True)\n",
    "        counts_2 /= np.sum(counts_2)\n",
    "        values.append(js_divergence(np.array(counts_1)  + epsilon,\n",
    "                                      np.array(counts_2) + epsilon))\n",
    "    return values\n",
    "\n",
    "def draw_display(dispsize: Tuple[int, int], imagefile=None):\n",
    "    # construct screen (black background)\n",
    "    # dots per inch\n",
    "    img = image.imread(imagefile)\n",
    "    dpi = 100.0\n",
    "    # determine the figure size in inches\n",
    "    figsize = (dispsize[0]/dpi, dispsize[1]/dpi)\n",
    "    # create a figure\n",
    "    fig = pyplot.figure(figsize=figsize, dpi=dpi, frameon=False)\n",
    "    ax = pyplot.Axes(fig, [0, 0, 1, 1])\n",
    "    ax.set_axis_off()\n",
    "    fig.add_axes(ax)\n",
    "    # plot display\n",
    "    ax.axis([0, dispsize[0], 0, dispsize[1]])\n",
    "    ax.imshow(img)\n",
    "\n",
    "    return fig, ax\n",
    "\n",
    "\n",
    "\n",
    "def pix2deg(pix, screenPX,screenCM,distanceCM, adjust_origin=True):\n",
    "    # Converts pixel screen coordinate to degrees of visual angle\n",
    "    # screenPX is the number of pixels that the monitor has in the horizontal\n",
    "    # axis (for x coord) or vertical axis (for y coord)\n",
    "    # screenCM is the width of the monitor in centimeters\n",
    "    # distanceCM is the distance of the monitor to the retina \n",
    "    # pix: screen coordinate in pixels\n",
    "    # adjust origin: if origin (0,0) of screen coordinates is in the corner of the screen rather than in the center, set to True to center coordinates\n",
    "    pix=np.array(pix)\n",
    "    # center screen coordinates such that 0 is center of the screen:\n",
    "    if adjust_origin: \n",
    "        pix = pix-(screenPX)/2 # pixel coordinates start with (0,0) \n",
    "    # eye-to-screen-distance in pixels of screen\n",
    "    distancePX = distanceCM*(screenPX/screenCM)\n",
    "    return np.arctan2(pix,distancePX) * 180/np.pi #  *180/pi wandelt bogenmass in grad\n",
    "\n",
    "\n",
    "def deg2pix(deg, screenPX, screenCM, distanceCM, adjust_origin = True, offsetCM = 0):\n",
    "    # Converts degrees of visual angle to pixel screen coordinates\n",
    "    # screenPX is the number of pixels that the monitor has in the horizontal\n",
    "    # screenCM is the width of the monitor in centimeters\n",
    "    # distanceCM is the distance of the monitor to the retina \n",
    "    phi = np.arctan2(1,distanceCM)*180/np.pi\n",
    "    pix = deg/(phi/(screenPX/(screenCM)))\n",
    "    if adjust_origin:\n",
    "        pix += (screenPX/2)\n",
    "    if offsetCM != 0:\n",
    "        offsetPX = offsetCM*(screenPX/screenCM)\n",
    "        pix += offsetPX\n",
    "    return pix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ce23a7a",
   "metadata": {},
   "source": [
    "### Params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "df9323b8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Device mapping:\n",
      "/job:localhost/replica:0/task:0/device:GPU:0 -> device: 0, name: A100-SXM4-40GB, pci bus id: 0000:47:00.0, compute capability: 8.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# select gpu to use for inference\n",
    "GPU = 2\n",
    "\n",
    "# paths for data and plots\n",
    "data_dir = 'data/'\n",
    "plot_dir = 'plots/'\n",
    "result_dir = 'results/'\n",
    "\n",
    "# modelnames\n",
    "model_name = 'GAN'\n",
    "real_name = 'real'\n",
    "gauss_name = 'Gauss'\n",
    "\n",
    "# specify if you want to use a GPU\n",
    "flag_train_on_gpu = True\n",
    "# specify if you want to recompute the fixations\n",
    "flag_recompute = False\n",
    "\n",
    "if flag_train_on_gpu:\n",
    "    import tensorflow as tf\n",
    "    # select graphic card\n",
    "    os.environ[\"CUDA_VISIBLE_DEVICES\"] = str(GPU)\n",
    "    os.environ[\"CUDA_DEVICE_ORDER\"] = \"PCI_BUS_ID\"\n",
    "    config = tf.compat.v1.ConfigProto(log_device_placement=True)\n",
    "    config.gpu_options.per_process_gpu_memory_fraction = 1.\n",
    "    config.gpu_options.allow_growth = True\n",
    "    tf_session = tf.compat.v1.Session(config=config)\n",
    "else:\n",
    "    import tensorflow as tf\n",
    "    # select graphic card\n",
    "    os.environ['CUDA_VISIBLE_DEVICES'] = '-1'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0019cb41",
   "metadata": {},
   "source": [
    "### Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "57975aa0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check if data was downloaded and processed\n",
    "if not os.path.exists(data_dir + '/GazeBase/raw/'):\n",
    "    os.system('python -m sp_eyegan.get_data --dataset-name GazeBase')\n",
    "    os.system('python -m sp_eyegan.create_event_data_from_gazebase --stimulus text')\n",
    "    \n",
    "# load data\n",
    "column_dict = joblib.load(data_dir + 'column_dict.joblib')\n",
    "fixation_matrix = np.load(data_dir + 'fixation_matrix_gazebase_vd_text.npy')\n",
    "saccade_matrix = np.load(data_dir  + 'saccade_matrix_gazebase_vd_text.npy')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cb9a5cf",
   "metadata": {},
   "source": [
    "### Set params for SP-EyeGAN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "31b17fbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "fix_window_size = fixation_matrix.shape[1]\n",
    "sac_window_size = saccade_matrix.shape[1]\n",
    "gen_kernel_sizes_fixation = [fix_window_size,8,4,2]\n",
    "gen_kernel_sizes_saccade = [sac_window_size,8,4,2]\n",
    "    \n",
    "# params for NN\n",
    "random_size = 32\n",
    "gen_filter_sizes = [16,8,4,2]\n",
    "channels = 2\n",
    "relu_in_last = False\n",
    "batch_size = 256\n",
    "\n",
    "dis_kernel_sizes = [8,16,32]\n",
    "dis_fiter_sizes = [32,64,128]\n",
    "dis_dropout = 0.3\n",
    "\n",
    "sample_size = 1000\n",
    "max_velocity = 0.5\n",
    "seed = 42\n",
    "np.random.seed(seed)\n",
    "\n",
    "model_config_fixation = {'gen_kernel_sizes':gen_kernel_sizes_fixation,\n",
    "                'gen_filter_sizes':gen_filter_sizes,\n",
    "                'dis_kernel_sizes':dis_kernel_sizes,\n",
    "                'dis_fiter_sizes':dis_fiter_sizes,\n",
    "                'dis_dropout':dis_dropout,\n",
    "                'window_size':fix_window_size,\n",
    "                'channels':channels,\n",
    "                'batch_size':batch_size,\n",
    "                'random_size':random_size,\n",
    "                'relu_in_last':relu_in_last,\n",
    "               }\n",
    "\n",
    "model_config_saccade = {'gen_kernel_sizes':gen_kernel_sizes_saccade,\n",
    "                'gen_filter_sizes':gen_filter_sizes,\n",
    "                'dis_kernel_sizes':dis_kernel_sizes,\n",
    "                'dis_fiter_sizes':dis_fiter_sizes,\n",
    "                'dis_dropout':dis_dropout,\n",
    "                'window_size':sac_window_size,\n",
    "                'channels':channels,\n",
    "                'batch_size':batch_size,\n",
    "                'random_size':random_size,\n",
    "                'relu_in_last':relu_in_last,\n",
    "               }\n",
    "\n",
    "\n",
    "fixation_path  = 'event_model/fixation_model_text'\n",
    "saccade_path   = 'event_model/saccade_model_text'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e3e17d93",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape of fixation_matrix: (393711, 100, 6)\n",
      "shape of saccade_matrix: (103883, 30, 6)\n"
     ]
    }
   ],
   "source": [
    "print('shape of fixation_matrix: ' + str(fixation_matrix.shape))\n",
    "print('shape of saccade_matrix: ' + str(saccade_matrix.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e5daa96",
   "metadata": {},
   "source": [
    "## Results for statistical baseline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8af70cf",
   "metadata": {},
   "source": [
    "#### Load data for statitical baseline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9aec8ac5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number fix points: 17889\n",
      "number sac points: 5294\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "357f98d2b5754a7c8318cc411f03f86a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/56504 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gen_data_statitical = pd.read_csv('data/statistical_baseline_fix_data.txt',sep=';',header = None)\n",
    "gen_data_statitical.columns = ['vel','type','frame']\n",
    "gen_data_statitical.head()\n",
    "\n",
    "\n",
    "fixation_index = 0\n",
    "saccade_index = 1\n",
    "velocities = np.array(gen_data_statitical['vel'])\n",
    "event_type = np.array(gen_data_statitical['type'])\n",
    "print('number fix points: ' + str(np.sum(event_type == fixation_index)))\n",
    "print('number sac points: ' + str(np.sum(event_type == saccade_index)))\n",
    "\n",
    "\n",
    "fix_vels_statistical = []\n",
    "sac_vels_statistical = []\n",
    "sac_acc_statistical  = []\n",
    "\n",
    "# transfrom °/s to °/ms\n",
    "scaling_factor = 1000.\n",
    "prev_event = -1\n",
    "cur_fixation = []\n",
    "cur_saccade  = []\n",
    "for i in tqdm(np.arange(len(velocities))):\n",
    "    cur_event = event_type[i]\n",
    "    cur_vel   = velocities[i]\n",
    "    if cur_event == prev_event:\n",
    "        if cur_event == fixation_index:\n",
    "            cur_fixation.append(cur_vel/scaling_factor)\n",
    "        elif cur_event == saccade_index:\n",
    "            cur_saccade.append(cur_vel/scaling_factor)\n",
    "    else:\n",
    "        if prev_event == fixation_index and len(cur_fixation) > 0:\n",
    "            fix_vels_statistical.append(np.array(cur_fixation))\n",
    "        elif prev_event == saccade_index and len(cur_saccade) > 0:\n",
    "            sac_vels_statistical.append(np.array(cur_saccade))\n",
    "            if len(cur_saccade) > 2:\n",
    "                sac_acc_statistical.append(np.array(cur_saccade)[1:] - np.array(cur_saccade)[0:-1])\n",
    "        \n",
    "        if cur_event == fixation_index:\n",
    "            cur_fixation = [cur_vel/scaling_factor]\n",
    "            cur_saccade = []\n",
    "        elif cur_event == saccade_index:\n",
    "            cur_saccade = [cur_vel/scaling_factor]\n",
    "            cur_fixation = []\n",
    "    prev_event = cur_event"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cca3b52",
   "metadata": {},
   "source": [
    "### Table 2 (Results for fixations)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4371e57",
   "metadata": {},
   "source": [
    "#### Fixations for EyeSyn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "cfc2bd8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.io import loadmat\n",
    "eye_syn_fixation_data = loadmat('data/eyeSyn_velocities.mat')\n",
    "eye_syn_fixation_data = eye_syn_fixation_data['out_data']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59d6b7b3",
   "metadata": {},
   "source": [
    "#### Fixations for SP-EyeGAN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2cbfd6ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "sp_fix_result_path = 'results/speyegan_fixations.joblib'\n",
    "if flag_recompute or not os.path.exists(sp_fix_result_path):\n",
    "    tf.keras.backend.clear_session()\n",
    "    fix_model = eventGAN.eventGAN(model_config_fixation)\n",
    "    fix_model.load_model(fixation_path)\n",
    "\n",
    "    noise = tf.random.normal([sample_size, random_size], seed = seed)\n",
    "    gen_fixations = np.array(fix_model.generator(noise, training=False),dtype=np.float32)\n",
    "    joblib.dump(gen_fixations, sp_fix_result_path, compress=3, protocol=2)\n",
    "else:\n",
    "    gen_fixations = joblib.load(sp_fix_result_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccf4a0f9",
   "metadata": {},
   "source": [
    "#### Fixations for VAE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "16d47a40",
   "metadata": {},
   "outputs": [],
   "source": [
    "vae_fix_result_path = 'results/vae_fixations.joblib'\n",
    "if flag_recompute or not os.path.exists(vae_fix_result_path):\n",
    "    tf.keras.backend.clear_session()\n",
    "    vae_model = vae.VAE(vae.get_vae_encoder(64, 2, 2), vae.get_vae_decoder(2, 2))\n",
    "    vae_model.load_model('event_model/vae_fixation_10')\n",
    "    noise = tf.random.normal([sample_size * 2, 2], seed = seed)\n",
    "    vae_fixations = np.array(vae_model.decoder(noise, training=False),dtype=np.float32)\n",
    "    vae_fixations = np.concatenate([vae_fixations[0:sample_size],\n",
    "                                    vae_fixations[sample_size:]],axis=1)\n",
    "    vae_fixations = vae_fixations[:,0:100,:] -0.5\n",
    "    joblib.dump(vae_fixations, vae_fix_result_path, compress=3, protocol=2)\n",
    "else:\n",
    "    vae_fixations = joblib.load(vae_fix_result_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38111d58",
   "metadata": {},
   "source": [
    "#### Original fixations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0ec81dbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "original_fix_result_path = 'results/original_fixations.joblib'\n",
    "if flag_recompute or not os.path.exists(original_fix_result_path):\n",
    "    rand_ids        = np.random.permutation(np.arange(fixation_matrix.shape[0]))\n",
    "    rand_ids_1      = rand_ids[0:sample_size]\n",
    "    rand_ids_2      = rand_ids[sample_size:sample_size+sample_size]\n",
    "    orig_fixations  = fixation_matrix[rand_ids_1,:,4:6]\n",
    "    orig_fixations[orig_fixations > max_velocity] = max_velocity\n",
    "    orig_fixations[orig_fixations < -max_velocity] = -max_velocity\n",
    "\n",
    "    orig_fixations_2 = fixation_matrix[rand_ids_2,:,4:6]\n",
    "    orig_fixations_2[orig_fixations_2 > max_velocity] = max_velocity\n",
    "    orig_fixations_2[orig_fixations_2 < -max_velocity] = -max_velocity\n",
    "    joblib.dump({'orig_fixations':orig_fixations,\n",
    "                 'orig_fixations_2':orig_fixations_2,\n",
    "                }, original_fix_result_path, compress=3, protocol=2)\n",
    "else:\n",
    "    orig_fix_data = joblib.load(original_fix_result_path)\n",
    "    orig_fixations_2 = orig_fix_data['orig_fixations_2']\n",
    "    orig_fixations = orig_fix_data['orig_fixations']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bc1db84",
   "metadata": {},
   "source": [
    "#### Create Table 2 with results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "df1cfa23",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9cea5716d6e74959856501c534ca2e74",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "30fdee301a984f19a41ab8a8a8124a33",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9c26a74c6fa84c5bb6c0d58b6bdb82a2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1477b2a3dd164c21993b94af23093b42",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "28f8af73921f4450973b7a39f045dc3e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "22151fcb2b374a5cabf9dbd2677ea8e4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fixation velocity\n",
      "JS(real || stat): 0.286 $\\pm$ 0.001\n",
      "JS(real || VAE): 0.204 $\\pm$ 0.001\n",
      "JS(real || EyeSyn): 0.065 $\\pm$ 0.001\n",
      "JS(real || GAN): \\textbf{0.03 $\\pm$ 0.001}$^*$\n",
      "JS(real || Real): 0.001 $\\pm$ 0.0\n",
      "Fixation mean velocity\n",
      "JS(real || stat): 0.692 $\\pm$ 0.004\n",
      "JS(real || VAE): 0.958 $\\pm$ 0.003\n",
      "JS(real || EyeSyn): 0.793 $\\pm$ 0.006\n",
      "JS(real || GAN): \\textbf{0.321 $\\pm$ 0.007}$^*$\n",
      "JS(real || Real): 0.073 $\\pm$ 0.003\n",
      "Fixation dispersion\n",
      "JS(real || stat): -\n",
      "JS(real || VAE): 0.739 $\\pm$ 0.006\n",
      "JS(real || EyeSyn): 0.991 $\\pm$ 0.002\n",
      "JS(real || GAN): \\textbf{0.315 $\\pm$ 0.007}$^*$\n",
      "JS(real || Real): 0.109 $\\pm$ 0.006\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Method</th>\n",
       "      <th>Velocity</th>\n",
       "      <th>Mean velocity</th>\n",
       "      <th>Disperson</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Statistical model</td>\n",
       "      <td>0.286 $\\pm$ 0.001</td>\n",
       "      <td>0.692 $\\pm$ 0.004</td>\n",
       "      <td>-</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>VAE</td>\n",
       "      <td>0.204 $\\pm$ 0.001</td>\n",
       "      <td>0.958 $\\pm$ 0.003</td>\n",
       "      <td>0.739 $\\pm$ 0.006</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>EyeSyn</td>\n",
       "      <td>0.065 $\\pm$ 0.001</td>\n",
       "      <td>0.793 $\\pm$ 0.006</td>\n",
       "      <td>0.991 $\\pm$ 0.002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>SpEyeGAN</td>\n",
       "      <td>0.03 $\\pm$ 0.001</td>\n",
       "      <td>0.321 $\\pm$ 0.007</td>\n",
       "      <td>0.315 $\\pm$ 0.007</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Real</td>\n",
       "      <td>0.001 $\\pm$ 0.0</td>\n",
       "      <td>0.073 $\\pm$ 0.003</td>\n",
       "      <td>0.109 $\\pm$ 0.006</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              Method           Velocity      Mean velocity          Disperson\n",
       "0  Statistical model  0.286 $\\pm$ 0.001  0.692 $\\pm$ 0.004                  -\n",
       "1                VAE  0.204 $\\pm$ 0.001  0.958 $\\pm$ 0.003  0.739 $\\pm$ 0.006\n",
       "2             EyeSyn  0.065 $\\pm$ 0.001  0.793 $\\pm$ 0.006  0.991 $\\pm$ 0.002\n",
       "3           SpEyeGAN   0.03 $\\pm$ 0.001  0.321 $\\pm$ 0.007  0.315 $\\pm$ 0.007\n",
       "4               Real    0.001 $\\pm$ 0.0  0.073 $\\pm$ 0.003  0.109 $\\pm$ 0.006"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fix_result_dict = dict()\n",
    "fix_result_dict['Method'] = ['Statistical model',\n",
    "                            'VAE',\n",
    "                            'EyeSyn',\n",
    "                            'SpEyeGAN',\n",
    "                            'Real']\n",
    "fix_result_dict['Velocity'] = []\n",
    "fix_result_dict['Mean velocity'] = []\n",
    "fix_result_dict['Disperson'] = []\n",
    "\n",
    "real_vels, real_x_vels, real_y_vels, real_dispersions = get_fixation_stats(orig_fixations)\n",
    "real_vels_2, real_x_vels_2, real_y_vels_2, real_dispersions_2 = get_fixation_stats(orig_fixations_2)\n",
    "fake_vels, fake_x_vels, fake_y_vels, fake_dispersions = get_fixation_stats(gen_fixations)\n",
    "eyesyn_vels, eyesyn_x_vels, eyesyn_y_vels, eyesyn_dispersions =  get_fixation_stats(eye_syn_fixation_data)\n",
    "vae_vels, vae_x_vels, vae_y_vels, vae_dispersions = get_fixation_stats(vae_fixations)\n",
    "\n",
    "\n",
    "x_vel_mean = np.nanmean([item for sublist in real_x_vels for item in sublist])\n",
    "x_vel_std = np.nanstd([item for sublist in real_x_vels for item in sublist])\n",
    "\n",
    "y_vel_mean = np.nanmean([item for sublist in real_y_vels for item in sublist])\n",
    "y_vel_std = np.nanstd([item for sublist in real_y_vels for item in sublist])\n",
    "gauss_fixations = np.concatenate([np.random.normal(x_vel_mean,x_vel_std,(sample_size,fix_window_size,1)),\n",
    "                                np.random.normal(y_vel_mean,y_vel_std,(sample_size,fix_window_size,1))],axis=2)\n",
    "gauss_vels, gauss_x_vels, gauss_y_vels, gauss_dispersions = get_fixation_stats(gauss_fixations)\n",
    "\n",
    "x_vel_mean = np.nanmean([item for sublist in real_x_vels for item in sublist])\n",
    "x_vel_std = np.nanstd([item for sublist in real_x_vels for item in sublist])\n",
    "#print(x_vel_mean, x_vel_std)\n",
    "\n",
    "y_vel_mean = np.nanmean([item for sublist in real_y_vels for item in sublist])\n",
    "y_vel_std = np.nanstd([item for sublist in real_y_vels for item in sublist])\n",
    "#print(y_vel_mean, y_vel_std)\n",
    "\n",
    "## velocity\n",
    "real   = np.array([item for sublist in real_vels for item in sublist])\n",
    "real_2 = np.array([item for sublist in real_vels_2 for item in sublist])\n",
    "fake   = np.array([item for sublist in fake_vels for item in sublist])\n",
    "gauss  = np.array([item for sublist in gauss_vels for item in sublist])\n",
    "eyesyn = np.array([item for sublist in eyesyn_vels for item in sublist])\n",
    "stat   = np.array([item for sublist in fix_vels_statistical for item in sublist])\n",
    "f_vae  = np.array([item for sublist in vae_vels for item in sublist])\n",
    "\n",
    "min_val = np.min(list(real) + list(fake) + list(gauss) + list(real_2) + list(eyesyn) + list(stat) + list(f_vae))\n",
    "max_val = np.max(list(real) + list(fake) + list(gauss) + list(real_2) + list(eyesyn) + list(stat) + list(f_vae))\n",
    "num_bins = 100\n",
    "epsilon = 0.00001\n",
    "number_per_iter = 10000\n",
    "iterations = 10\n",
    "bins = np.linspace(min_val,max_val,num_bins)\n",
    "\n",
    "stat_vals = js_divergence_sampling( real, stat, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "vae_vals = js_divergence_sampling( real, f_vae, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "eyesyn_vals = js_divergence_sampling( real, eyesyn, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "gan_vals = js_divergence_sampling( real, fake, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "real_vals = js_divergence_sampling( real, real_2, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "\n",
    "model_vals = [stat_vals, vae_vals, eyesyn_vals, gan_vals]\n",
    "model_means = [np.mean(a) for a in model_vals]\n",
    "prefixes = ['' for a in range(len(model_means))]\n",
    "suffixes = ['' for a in range(len(model_means))]\n",
    "arg_sort = np.argsort(model_means)\n",
    "prefixes[arg_sort[0]] = '\\\\textbf{'\n",
    "suffixes[arg_sort[0]] = '}'\n",
    "\n",
    "tt_test_pvalue = ttest_ind(model_vals[arg_sort[0]],model_vals[arg_sort[1]],alternative='two-sided').pvalue\n",
    "if tt_test_pvalue <= 0.05:\n",
    "    suffixes[arg_sort[0]] += '$^*$'\n",
    "\n",
    "                                        \n",
    "print('Fixation velocity')\n",
    "print('JS(real || stat): ' + prefixes[0] + str(np.round(np.mean(stat_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(stat_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[0])\n",
    "fix_result_dict['Velocity'].append(str(np.round(np.mean(stat_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(stat_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || VAE): ' + prefixes[1] +  str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[1])\n",
    "fix_result_dict['Velocity'].append(str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || EyeSyn): ' + prefixes[2] +  str(np.round(np.mean(eyesyn_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(eyesyn_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[2])\n",
    "fix_result_dict['Velocity'].append(str(np.round(np.mean(eyesyn_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(eyesyn_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || GAN): ' + prefixes[3] +  str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[3])\n",
    "fix_result_dict['Velocity'].append(str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || Real): ' + str(np.round(np.mean(real_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "fix_result_dict['Velocity'].append(str(np.round(np.mean(real_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "\n",
    "    \n",
    "# mean velocity\n",
    "real  = np.array([np.mean(a) for a in real_vels])\n",
    "real_2 = np.array([np.mean(a) for a in real_vels_2])\n",
    "fake  = np.array([np.mean(a) for a in fake_vels])\n",
    "gauss = np.array([np.mean(a) for a in gauss_vels])\n",
    "eyesyn = np.array([np.mean(a) for a in eyesyn_vels])\n",
    "f_vae  = np.array([np.mean(a) for a in vae_vels])\n",
    "\n",
    "min_val = np.min(list(real) + list(fake) + list(gauss) + list(real_2) + list(eyesyn) + list(stat) + list(f_vae))\n",
    "max_val = np.max(list(real) + list(fake) + list(gauss) + list(real_2) + list(eyesyn) + list(stat) + list(f_vae))\n",
    "num_bins = 100\n",
    "epsilon = 0.00001\n",
    "number_per_iter = 500\n",
    "iterations = 10\n",
    "bins = np.linspace(min_val,max_val,num_bins)\n",
    "\n",
    "\n",
    "stat_vals = js_divergence_sampling( real, stat, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "vae_vals = js_divergence_sampling( real, f_vae, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "eyesyn_vals = js_divergence_sampling( real, eyesyn, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "gan_vals = js_divergence_sampling( real, fake, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "real_vals = js_divergence_sampling( real, real_2, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "                                        \n",
    "model_vals = [stat_vals, vae_vals, eyesyn_vals, gan_vals]\n",
    "model_means = [np.mean(a) for a in model_vals]\n",
    "prefixes = ['' for a in range(len(model_means))]\n",
    "suffixes = ['' for a in range(len(model_means))]\n",
    "arg_sort = np.argsort(model_means)\n",
    "prefixes[arg_sort[0]] = '\\\\textbf{'\n",
    "suffixes[arg_sort[0]] = '}'\n",
    "\n",
    "tt_test_pvalue = ttest_ind(model_vals[arg_sort[0]],model_vals[arg_sort[1]],alternative='two-sided').pvalue\n",
    "if tt_test_pvalue <= 0.05:\n",
    "    suffixes[arg_sort[0]] += '$^*$'\n",
    "\n",
    "                                        \n",
    "print('Fixation mean velocity')\n",
    "print('JS(real || stat): ' + prefixes[0] + str(np.round(np.mean(stat_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(stat_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[0])\n",
    "fix_result_dict['Mean velocity'].append(str(np.round(np.mean(stat_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(stat_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || VAE): ' + prefixes[1] +  str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[1])\n",
    "fix_result_dict['Mean velocity'].append(str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || EyeSyn): ' + prefixes[2] +  str(np.round(np.mean(eyesyn_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(eyesyn_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[2])\n",
    "fix_result_dict['Mean velocity'].append(str(np.round(np.mean(eyesyn_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(eyesyn_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || GAN): ' + prefixes[3] +  str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[3])\n",
    "fix_result_dict['Mean velocity'].append(str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || Real): ' + str(np.round(np.mean(real_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "fix_result_dict['Mean velocity'].append(str(np.round(np.mean(real_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "\n",
    "\n",
    "\n",
    "# dispersion\n",
    "real   = np.array(real_dispersions)\n",
    "real_2 = np.array(real_dispersions_2)\n",
    "fake   = np.array(fake_dispersions)\n",
    "gauss  = np.array(gauss_dispersions)\n",
    "eyesyn = np.array(eyesyn_dispersions)\n",
    "f_vae  = np.array(vae_dispersions)\n",
    "\n",
    "min_val = np.min(list(real) + list(fake) + list(gauss) + list(real_2) + list(eyesyn) + list(f_vae))\n",
    "max_val = np.max(list(real) + list(fake) + list(gauss) + list(real_2) + list(eyesyn) + list(f_vae))\n",
    "num_bins = 100\n",
    "epsilon = 0.00001\n",
    "number_per_iter = 500\n",
    "iterations = 10\n",
    "bins = np.linspace(min_val,max_val,num_bins)\n",
    "\n",
    "\n",
    "vae_vals = js_divergence_sampling( real, f_vae, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "eyesyn_vals = js_divergence_sampling( real, eyesyn, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "gan_vals = js_divergence_sampling( real, fake, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "real_vals = js_divergence_sampling( real, real_2, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "\n",
    "model_vals = [vae_vals, eyesyn_vals, gan_vals]\n",
    "model_means = [np.mean(a) for a in model_vals]\n",
    "prefixes = ['' for a in range(len(model_means))]\n",
    "suffixes = ['' for a in range(len(model_means))]\n",
    "arg_sort = np.argsort(model_means)\n",
    "prefixes[arg_sort[0]] = '\\\\textbf{'\n",
    "suffixes[arg_sort[0]] = '}'\n",
    "\n",
    "tt_test_pvalue = ttest_ind(model_vals[arg_sort[0]],model_vals[arg_sort[1]],alternative='two-sided').pvalue\n",
    "if tt_test_pvalue <= 0.05:\n",
    "    suffixes[arg_sort[0]] += '$^*$'\n",
    "\n",
    "                                        \n",
    "print('Fixation dispersion')\n",
    "print('JS(real || stat): -')\n",
    "fix_result_dict['Disperson'].append('-')\n",
    "print('JS(real || VAE): ' + prefixes[0] +  str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[0])\n",
    "fix_result_dict['Disperson'].append(str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || EyeSyn): ' + prefixes[1] +  str(np.round(np.mean(eyesyn_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(eyesyn_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[1])\n",
    "fix_result_dict['Disperson'].append(str(np.round(np.mean(eyesyn_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(eyesyn_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || GAN): ' + prefixes[2] +  str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[2])\n",
    "fix_result_dict['Disperson'].append(str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || Real): ' + str(np.round(np.mean(real_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "fix_result_dict['Disperson'].append(str(np.round(np.mean(real_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "\n",
    "pd.DataFrame(fix_result_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b0bb9e2",
   "metadata": {},
   "source": [
    "### Table 3 (Results for saccades)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa26650b",
   "metadata": {},
   "source": [
    "#### Saccades for SP-EyeGAN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f3195cb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "sp_sac_result_path = 'results/speyegan_saccades.joblib'\n",
    "if flag_recompute or not os.path.exists(sp_sac_result_path):\n",
    "    tf.keras.backend.clear_session()\n",
    "    sac_model = eventGAN.eventGAN(model_config_saccade)\n",
    "    sac_model.load_model(saccade_path)\n",
    "\n",
    "    noise = tf.random.normal([sample_size, random_size], seed = seed)\n",
    "    gen_saccades = np.array(sac_model.generator(noise, training=False),dtype=np.float32)\n",
    "    joblib.dump(gen_saccades, sp_sac_result_path, compress=3, protocol=2)\n",
    "else:\n",
    "    gen_saccades = joblib.load(sp_sac_result_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2c09983",
   "metadata": {},
   "source": [
    "#### Saccades for VAE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "889ad606",
   "metadata": {},
   "outputs": [],
   "source": [
    "vae_sac_result_path = 'results/vae_saccades.joblib'\n",
    "if flag_recompute or not os.path.exists(vae_sac_result_path):\n",
    "    tf.keras.backend.clear_session()\n",
    "    vae_model = vae.VAE(vae.get_vae_encoder(64, 2, 2), vae.get_vae_decoder(2, 2))\n",
    "    vae_model.load_model('event_model/vae_saccade_10')\n",
    "    noise = tf.random.normal([sample_size, 2], seed = seed)\n",
    "    vae_saccades = np.array(vae_model.decoder(noise, training=False),dtype=np.float32)\n",
    "    vae_saccades = vae_saccades[:,0:30,:] -0.5\n",
    "    joblib.dump(vae_saccades, vae_sac_result_path, compress=3, protocol=2)\n",
    "else:\n",
    "    vae_saccades = joblib.load(vae_sac_result_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fa8ec39",
   "metadata": {},
   "source": [
    "#### Original saccades"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "f60fadb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "original_sac_result_path = 'results/original_saccades.joblib'\n",
    "if flag_recompute or not os.path.exists(original_sac_result_path):\n",
    "    rand_ids        = np.random.permutation(np.arange(saccade_matrix.shape[0]))\n",
    "    rand_ids_1      = rand_ids[0:sample_size]\n",
    "    rand_ids_2      = rand_ids[sample_size:sample_size+sample_size]\n",
    "\n",
    "    rand_ids        = np.arange(saccade_matrix.shape[0])\n",
    "    rand_ids        = np.random.permutation(rand_ids)[0:sample_size]\n",
    "    orig_saccades   = saccade_matrix[rand_ids_1,:,4:6]\n",
    "    orig_saccades[orig_saccades > max_velocity] = max_velocity\n",
    "    orig_saccades[orig_saccades < -max_velocity] = -max_velocity\n",
    "\n",
    "    orig_saccades_2 = saccade_matrix[rand_ids_2,:,4:6]\n",
    "    orig_saccades_2[orig_saccades_2 > max_velocity] = max_velocity\n",
    "    orig_saccades_2[orig_saccades_2 < -max_velocity] = -max_velocity\n",
    "    joblib.dump({'orig_saccades':orig_saccades,\n",
    "                 'orig_saccades_2':orig_saccades_2,\n",
    "                }, original_sac_result_path, compress=3, protocol=2)\n",
    "else:\n",
    "    orig_sac_data = joblib.load(original_sac_result_path)\n",
    "    orig_saccades = orig_sac_data['orig_saccades']\n",
    "    orig_saccades_2 = orig_sac_data['orig_saccades_2']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df1391d3",
   "metadata": {},
   "source": [
    "#### Create resulst for Table 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "6f0c3a02",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8e785ab6710d477cbf7ca41943343150",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fbc2e018205b40af946f9efd76f4b465",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e8b65a769d9e4a24a9bde7808ae2fda0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e50940a893274e4b8a4016a442c4eee0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cbc41416f2f24a7bad71a51a270361c7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saccade peak velocity\n",
      "JS(real || stat): 0.426 $\\pm$ 0.007\n",
      "JS(real || VAE): 0.951 $\\pm$ 0.004\n",
      "JS(real || GAN): \\textbf{0.399 $\\pm$ 0.006}$^*$\n",
      "JS(real || Real): 0.079 $\\pm$ 0.005\n",
      "Saccade mean velocity\n",
      "JS(real || stat): 0.253 $\\pm$ 0.007\n",
      "JS(real || VAE): 0.906 $\\pm$ 0.007\n",
      "JS(real || GAN): \\textbf{0.251 $\\pm$ 0.004}\n",
      "JS(real || Real): 0.052 $\\pm$ 0.003\n",
      "Saccade peak acceleration\n",
      "JS(real || stat): 0.916 $\\pm$ 0.004\n",
      "JS(real || VAE): 0.882 $\\pm$ 0.004\n",
      "JS(real || GAN): \\textbf{0.294 $\\pm$ 0.004}$^*$\n",
      "JS(real || Real): 0.249 $\\pm$ 0.006\n",
      "Saccade mean acceleration\n",
      "JS(real || stat): 0.859 $\\pm$ 0.004\n",
      "JS(real || VAE): 0.871 $\\pm$ 0.006\n",
      "JS(real || GAN): \\textbf{0.24 $\\pm$ 0.004}$^*$\n",
      "JS(real || Real): 0.046 $\\pm$ 0.002\n",
      "Saccade amplitude\n",
      "JS(real || stat): -\n",
      "JS(real || VAE): 0.92 $\\pm$ 0.008\n",
      "JS(real || GAN): \\textbf{0.267 $\\pm$ 0.006}$^*$\n",
      "JS(real || Real): 0.075 $\\pm$ 0.003\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Method</th>\n",
       "      <th>Peak velocity</th>\n",
       "      <th>Mean velocity</th>\n",
       "      <th>Peak acceleration</th>\n",
       "      <th>Mean acceleration</th>\n",
       "      <th>Amplitude</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Statistical model</td>\n",
       "      <td>0.426 $\\pm$ 0.007</td>\n",
       "      <td>0.253 $\\pm$ 0.007</td>\n",
       "      <td>0.916 $\\pm$ 0.004</td>\n",
       "      <td>0.859 $\\pm$ 0.004</td>\n",
       "      <td>-</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>VAE</td>\n",
       "      <td>0.951 $\\pm$ 0.004</td>\n",
       "      <td>0.906 $\\pm$ 0.007</td>\n",
       "      <td>0.882 $\\pm$ 0.004</td>\n",
       "      <td>0.871 $\\pm$ 0.006</td>\n",
       "      <td>0.92 $\\pm$ 0.008</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>SpEyeGAN</td>\n",
       "      <td>0.399 $\\pm$ 0.006</td>\n",
       "      <td>0.251 $\\pm$ 0.004</td>\n",
       "      <td>0.294 $\\pm$ 0.004</td>\n",
       "      <td>0.24 $\\pm$ 0.004</td>\n",
       "      <td>0.267 $\\pm$ 0.006</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Real</td>\n",
       "      <td>0.079 $\\pm$ 0.005</td>\n",
       "      <td>0.052 $\\pm$ 0.003</td>\n",
       "      <td>0.249 $\\pm$ 0.006</td>\n",
       "      <td>0.046 $\\pm$ 0.002</td>\n",
       "      <td>0.075 $\\pm$ 0.003</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              Method      Peak velocity      Mean velocity  Peak acceleration  \\\n",
       "0  Statistical model  0.426 $\\pm$ 0.007  0.253 $\\pm$ 0.007  0.916 $\\pm$ 0.004   \n",
       "1                VAE  0.951 $\\pm$ 0.004  0.906 $\\pm$ 0.007  0.882 $\\pm$ 0.004   \n",
       "2           SpEyeGAN  0.399 $\\pm$ 0.006  0.251 $\\pm$ 0.004  0.294 $\\pm$ 0.004   \n",
       "3               Real  0.079 $\\pm$ 0.005  0.052 $\\pm$ 0.003  0.249 $\\pm$ 0.006   \n",
       "\n",
       "   Mean acceleration          Amplitude  \n",
       "0  0.859 $\\pm$ 0.004                  -  \n",
       "1  0.871 $\\pm$ 0.006   0.92 $\\pm$ 0.008  \n",
       "2   0.24 $\\pm$ 0.004  0.267 $\\pm$ 0.006  \n",
       "3  0.046 $\\pm$ 0.002  0.075 $\\pm$ 0.003  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sac_result_dict = dict()\n",
    "sac_result_dict['Method'] = ['Statistical model',\n",
    "                            'VAE',\n",
    "                            'SpEyeGAN',\n",
    "                            'Real']\n",
    "sac_result_dict['Peak velocity'] = []\n",
    "sac_result_dict['Mean velocity'] = []\n",
    "sac_result_dict['Peak acceleration'] = []\n",
    "sac_result_dict['Mean acceleration'] = []\n",
    "sac_result_dict['Amplitude'] = []\n",
    "\n",
    "real_vels_sac, real_x_vels_sac, real_y_vels_sac, real_accs_sac, real_x_accs_sac, real_y_accs_sac, real_amplitudes_sac      = get_saccade_stats(orig_saccades)\n",
    "real_vels_sac_2, real_x_vels_sac_2, real_y_vels_sac_2, real_accs_sac_2, real_x_accs_sac_2, real_y_accs_sac_2, real_amplitudes_sac_2      = get_saccade_stats(orig_saccades_2)\n",
    "fake_vels_sac, fake_x_vels_sac, fake_y_vels_sac, fake_accs_sac, fake_real_x_accs_sac, fake_y_accs_sac, fake_amplitudes_sac = get_saccade_stats(gen_saccades)\n",
    "vae_vels_sac, vae_x_vels_sac, vae_y_vels_sac, vae_accs_sac, vae_x_accs_sac, vae_y_accs_sac, vae_amplitudes_sac = get_saccade_stats(vae_saccades)\n",
    "\n",
    "\n",
    "x_vel_mean_sac = np.nanmean([item for sublist in real_x_vels_sac for item in sublist])\n",
    "x_vel_std_sac = np.nanstd([item for sublist in real_x_vels_sac for item in sublist])\n",
    "\n",
    "y_vel_mean_sac = np.nanmean([item for sublist in real_y_vels_sac for item in sublist])\n",
    "y_vel_std_sac = np.nanstd([item for sublist in real_y_vels_sac for item in sublist])\n",
    "gauss_saccades = np.concatenate([np.random.normal(x_vel_mean_sac,x_vel_std_sac,(sample_size,sac_window_size,1)),\n",
    "                                np.random.normal(y_vel_mean_sac,y_vel_std_sac,(sample_size,sac_window_size,1))],axis=2)\n",
    "gauss_vels_sac, gauss_x_vels_sac, gauss_y_vels_sac, gauss_accs_sac, gauss_real_x_accs_sac, gauss_y_accs_sac, gauss_amplitudes_sac = get_saccade_stats(gauss_saccades)\n",
    "\n",
    "\n",
    "# peak velocity\n",
    "real = np.array([np.max(a) for a in real_vels_sac])\n",
    "real_2 = np.array([np.max(a) for a in real_vels_sac_2])\n",
    "fake = np.array([np.max(a) for a in fake_vels_sac])\n",
    "gauss = np.array([np.max(a) for a in gauss_vels_sac])\n",
    "stat   = np.array([np.max(a) for a in sac_vels_statistical])\n",
    "f_vae  = np.array([np.max(a) for a in vae_vels_sac])\n",
    "\n",
    "min_val = np.min(list(real) + list(fake)  + list(real_2)  + list(stat) + list(f_vae))\n",
    "max_val = np.max(list(real) + list(fake)  + list(real_2)  + list(stat) + list(f_vae))\n",
    "num_bins = 100\n",
    "epsilon = 0.00001\n",
    "number_per_iter = 500\n",
    "iterations = 10\n",
    "bins = np.linspace(min_val,max_val,num_bins)\n",
    "\n",
    "print('Saccade peak velocity')\n",
    "stat_vals = js_divergence_sampling( real, stat, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "vae_vals = js_divergence_sampling( real, f_vae, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "gan_vals = js_divergence_sampling( real, fake, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "real_vals = js_divergence_sampling( real, real_2, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "\n",
    "model_vals = [stat_vals, vae_vals, gan_vals]\n",
    "model_means = [np.mean(a) for a in model_vals]\n",
    "prefixes = ['' for a in range(len(model_means))]\n",
    "suffixes = ['' for a in range(len(model_means))]\n",
    "arg_sort = np.argsort(model_means)\n",
    "prefixes[arg_sort[0]] = '\\\\textbf{'\n",
    "suffixes[arg_sort[0]] = '}'\n",
    "\n",
    "tt_test_pvalue = ttest_ind(model_vals[arg_sort[0]],model_vals[arg_sort[1]],alternative='two-sided').pvalue\n",
    "if tt_test_pvalue <= 0.05:\n",
    "    suffixes[arg_sort[0]] += '$^*$'\n",
    "                                       \n",
    "print('JS(real || stat): ' + prefixes[0] + str(np.round(np.mean(stat_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(stat_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[0])\n",
    "sac_result_dict['Peak velocity'].append(str(np.round(np.mean(stat_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(stat_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || VAE): ' + prefixes[1] +  str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[1])\n",
    "sac_result_dict['Peak velocity'].append(str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || GAN): ' + prefixes[2] +  str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[2])\n",
    "sac_result_dict['Peak velocity'].append(str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || Real): ' + str(np.round(np.mean(real_vals),decimals = 3)) + ' $\\\\pm$ ' +\\\n",
    "      str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "sac_result_dict['Peak velocity'].append(str(np.round(np.mean(real_vals),decimals = 3)) + ' $\\\\pm$ ' +\\\n",
    "      str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "\n",
    "# mean velocity\n",
    "real = np.array([np.mean(a) for a in real_vels_sac])\n",
    "real_2 = np.array([np.mean(a) for a in real_vels_sac_2])\n",
    "fake = np.array([np.mean(a) for a in fake_vels_sac])\n",
    "gauss = np.array([np.mean(a) for a in gauss_vels_sac])\n",
    "stat   = np.array([np.mean(a) for a in sac_vels_statistical])\n",
    "f_vae  = np.array([np.mean(a) for a in vae_vels_sac])\n",
    "\n",
    "min_val = np.min(list(real) + list(fake)  + list(real_2)  + list(stat) + list(f_vae))\n",
    "max_val = np.max(list(real) + list(fake)  + list(real_2)  + list(stat) + list(f_vae))\n",
    "num_bins = 100\n",
    "epsilon = 0.00001\n",
    "number_per_iter = 500\n",
    "iterations = 10\n",
    "bins = np.linspace(min_val,max_val,num_bins)\n",
    "\n",
    "print('Saccade mean velocity')\n",
    "stat_vals = js_divergence_sampling( real, stat, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "vae_vals = js_divergence_sampling( real, f_vae, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "gan_vals = js_divergence_sampling( real, fake, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "real_vals = js_divergence_sampling( real, real_2, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "                                        \n",
    "model_vals = [stat_vals, vae_vals, gan_vals]\n",
    "model_means = [np.mean(a) for a in model_vals]\n",
    "prefixes = ['' for a in range(len(model_means))]\n",
    "suffixes = ['' for a in range(len(model_means))]\n",
    "arg_sort = np.argsort(model_means)\n",
    "prefixes[arg_sort[0]] = '\\\\textbf{'\n",
    "suffixes[arg_sort[0]] = '}'\n",
    "\n",
    "tt_test_pvalue = ttest_ind(model_vals[arg_sort[0]],model_vals[arg_sort[1]],alternative='two-sided').pvalue\n",
    "if tt_test_pvalue <= 0.05:\n",
    "    suffixes[arg_sort[0]] += '$^*$'\n",
    "                                       \n",
    "print('JS(real || stat): ' + prefixes[0] + str(np.round(np.mean(stat_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(stat_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[0])\n",
    "sac_result_dict['Mean velocity'].append(str(np.round(np.mean(stat_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(stat_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || VAE): ' + prefixes[1] +  str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[1])\n",
    "sac_result_dict['Mean velocity'].append(str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || GAN): ' + prefixes[2] +  str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[2])\n",
    "sac_result_dict['Mean velocity'].append(str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || Real): ' + str(np.round(np.mean(real_vals),decimals = 3)) + ' $\\\\pm$ ' +\\\n",
    "      str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "sac_result_dict['Mean velocity'].append(str(np.round(np.mean(real_vals),decimals = 3)) + ' $\\\\pm$ ' +\\\n",
    "      str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "\n",
    "\n",
    "# peak acceleration\n",
    "real = np.array([np.max(a) for a in real_accs_sac])\n",
    "real_2 = np.array([np.mean(a) for a in real_accs_sac])\n",
    "fake = np.array([np.max(a) for a in fake_accs_sac])\n",
    "gauss = np.array([np.max(a) for a in gauss_accs_sac])\n",
    "stat   = np.array([np.max(a) for a in sac_acc_statistical])\n",
    "f_vae = np.array([np.max(a) for a in vae_accs_sac])\n",
    "\n",
    "min_val = np.min(list(real) + list(fake) + list(gauss) + list(real_2) + list(eyesyn) + list(stat) + list(f_vae))\n",
    "max_val = np.max(list(real) + list(fake) + list(gauss) + list(real_2) + list(eyesyn) + list(stat) + list(f_vae))\n",
    "num_bins = 100\n",
    "epsilon = 0.00001\n",
    "bins = np.linspace(min_val,max_val,num_bins)\n",
    "\n",
    "\n",
    "print('Saccade peak acceleration')\n",
    "stat_vals = js_divergence_sampling( real, stat, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "vae_vals = js_divergence_sampling( real, f_vae, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "gan_vals = js_divergence_sampling( real, fake, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "real_vals = js_divergence_sampling( real, real_2, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "                                        \n",
    "model_vals = [stat_vals, vae_vals, gan_vals]\n",
    "model_means = [np.mean(a) for a in model_vals]\n",
    "prefixes = ['' for a in range(len(model_means))]\n",
    "suffixes = ['' for a in range(len(model_means))]\n",
    "arg_sort = np.argsort(model_means)\n",
    "prefixes[arg_sort[0]] = '\\\\textbf{'\n",
    "suffixes[arg_sort[0]] = '}'\n",
    "\n",
    "tt_test_pvalue = ttest_ind(model_vals[arg_sort[0]],model_vals[arg_sort[1]],alternative='two-sided').pvalue\n",
    "if tt_test_pvalue <= 0.05:\n",
    "    suffixes[arg_sort[0]] += '$^*$'\n",
    "                                       \n",
    "print('JS(real || stat): ' + prefixes[0] + str(np.round(np.mean(stat_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(stat_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[0])\n",
    "sac_result_dict['Peak acceleration'].append( str(np.round(np.mean(stat_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(stat_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || VAE): ' + prefixes[1] +  str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[1])\n",
    "sac_result_dict['Peak acceleration'].append(str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || GAN): ' + prefixes[2] +  str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[2])\n",
    "sac_result_dict['Peak acceleration'].append(str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || Real): ' + str(np.round(np.mean(real_vals),decimals = 3)) + ' $\\\\pm$ ' +\\\n",
    "      str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "sac_result_dict['Peak acceleration'].append(str(np.round(np.mean(real_vals),decimals = 3)) + ' $\\\\pm$ ' +\\\n",
    "      str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "\n",
    "# mean acceleration\n",
    "real = np.array([np.mean(a) for a in real_accs_sac])\n",
    "real_2 = np.array([np.mean(a) for a in real_accs_sac])\n",
    "fake = np.array([np.mean(a) for a in fake_accs_sac])\n",
    "gauss = np.array([np.mean(a) for a in gauss_accs_sac])\n",
    "stat   = np.array([np.mean(a) for a in sac_acc_statistical])\n",
    "f_vae = np.array([np.mean(a) for a in vae_accs_sac])\n",
    "\n",
    "min_val = np.min(list(real) + list(fake) + list(gauss) + list(real_2) + list(eyesyn) + list(stat) + list(f_vae))\n",
    "max_val = np.max(list(real) + list(fake) + list(gauss) + list(real_2) + list(eyesyn) + list(stat) + list(f_vae))\n",
    "num_bins = 100\n",
    "epsilon = 0.00001\n",
    "bins = np.linspace(min_val,max_val,num_bins)\n",
    "\n",
    "\n",
    "\n",
    "print('Saccade mean acceleration')\n",
    "stat_vals = js_divergence_sampling( real, stat, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "vae_vals = js_divergence_sampling( real, f_vae, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "gan_vals = js_divergence_sampling( real, fake, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "real_vals = js_divergence_sampling( real, real_2, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "\n",
    "model_vals = [stat_vals, vae_vals, gan_vals]\n",
    "model_means = [np.mean(a) for a in model_vals]\n",
    "prefixes = ['' for a in range(len(model_means))]\n",
    "suffixes = ['' for a in range(len(model_means))]\n",
    "arg_sort = np.argsort(model_means)\n",
    "prefixes[arg_sort[0]] = '\\\\textbf{'\n",
    "suffixes[arg_sort[0]] = '}'\n",
    "\n",
    "tt_test_pvalue = ttest_ind(model_vals[arg_sort[0]],model_vals[arg_sort[1]],alternative='two-sided').pvalue\n",
    "if tt_test_pvalue <= 0.05:\n",
    "    suffixes[arg_sort[0]] += '$^*$'\n",
    "                                       \n",
    "print('JS(real || stat): ' + prefixes[0] + str(np.round(np.mean(stat_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(stat_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[0])\n",
    "sac_result_dict['Mean acceleration'].append(str(np.round(np.mean(stat_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(stat_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || VAE): ' + prefixes[1] +  str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[1])\n",
    "sac_result_dict['Mean acceleration'].append(str(np.round(np.mean(vae_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || GAN): ' + prefixes[2] +  str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[2])\n",
    "sac_result_dict['Mean acceleration'].append(str(np.round(np.mean(gan_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || Real): ' + str(np.round(np.mean(real_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "sac_result_dict['Mean acceleration'].append(str(np.round(np.mean(real_vals),decimals = 3)) +\\\n",
    "      ' $\\\\pm$ ' + str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "\n",
    "\n",
    "# amplitude\n",
    "real = real_amplitudes_sac\n",
    "real_2 = np.array([np.mean(a) for a in real_amplitudes_sac_2])\n",
    "fake = fake_amplitudes_sac\n",
    "gauss = gauss_amplitudes_sac\n",
    "f_vae = vae_amplitudes_sac\n",
    "\n",
    "min_val = np.min(list(real) + list(fake) + list(gauss) + list(real_2) + list(f_vae))\n",
    "max_val = np.max(list(real) + list(fake) + list(gauss) + list(real_2) + list(f_vae))\n",
    "num_bins = 100\n",
    "epsilon = 0.00001\n",
    "bins = np.linspace(min_val,max_val,num_bins)\n",
    "\n",
    "\n",
    "\n",
    "print('Saccade amplitude')\n",
    "vae_vals = js_divergence_sampling( real, f_vae, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "gan_vals = js_divergence_sampling( real, fake, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "real_vals = js_divergence_sampling( real, real_2, bins, epsilon =  epsilon,\n",
    "                                        iterations = iterations, number_per_iter = number_per_iter)\n",
    "\n",
    "model_vals = [vae_vals, gan_vals]\n",
    "model_means = [np.mean(a) for a in model_vals]\n",
    "prefixes = ['' for a in range(len(model_means))]\n",
    "suffixes = ['' for a in range(len(model_means))]\n",
    "arg_sort = np.argsort(model_means)\n",
    "prefixes[arg_sort[0]] = '\\\\textbf{'\n",
    "suffixes[arg_sort[0]] = '}'\n",
    "\n",
    "tt_test_pvalue = ttest_ind(model_vals[arg_sort[0]],model_vals[arg_sort[1]],alternative='two-sided').pvalue\n",
    "if tt_test_pvalue <= 0.05:\n",
    "    suffixes[arg_sort[0]] += '$^*$'\n",
    "                                       \n",
    "print('JS(real || stat): -')\n",
    "sac_result_dict['Amplitude'].append('-')\n",
    "print('JS(real || VAE): ' + prefixes[0] +  str(np.round(np.mean(vae_vals),decimals = 3)) + ' $\\\\pm$ ' +\\\n",
    "      str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[0])\n",
    "sac_result_dict['Amplitude'].append(str(np.round(np.mean(vae_vals),decimals = 3)) + ' $\\\\pm$ ' +\\\n",
    "      str(np.round(np.std(vae_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || GAN): ' + prefixes[1] +  str(np.round(np.mean(gan_vals),decimals = 3)) + ' $\\\\pm$ ' +\\\n",
    "      str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)) + suffixes[1])\n",
    "sac_result_dict['Amplitude'].append(str(np.round(np.mean(gan_vals),decimals = 3)) + ' $\\\\pm$ ' +\\\n",
    "      str(np.round(np.std(gan_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "print('JS(real || Real): ' + str(np.round(np.mean(real_vals),decimals = 3)) + ' $\\\\pm$ ' +\\\n",
    "      str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "sac_result_dict['Amplitude'].append(str(np.round(np.mean(real_vals),decimals = 3)) + ' $\\\\pm$ ' +\\\n",
    "      str(np.round(np.std(real_vals)/np.sqrt(iterations),decimals = 3)))\n",
    "\n",
    "pd.DataFrame(sac_result_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e635d42",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}