SCL_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This is the notebook we use to generate the Fig.1 in our paper, this notebook is inspired by \n",
    "#https://keras.io/examples/vision/supervised-contrastive-learning/#1-supervised-contrastive-learning-loss-function\n",
    "#to run this notebook, you need install tensorflow, tensorflow_addons as well as sklearn\n",
    "import os \n",
    "from pathlib import Path\n",
    "import random\n",
    "\n",
    "import tensorflow as tf \n",
    "import tensorflow_addons as tfa\n",
    "from tensorflow import keras \n",
    "\n",
    "import numpy as np \n",
    "import matplotlib.pyplot as plt \n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "import os \n",
    "%matplotlib inline\n",
    "\n",
    "#We assume the Feature dict to have this structure\n",
    "#  Feature/\n",
    "#    |- Train\n",
    "#       |- Genuine\n",
    "#           /- LA_T_3100997.npz (This is the LFCC feature of the utterance LA_T_3100997.flac extracted by the baseline program of ASVspoof 2019 LA)\n",
    "#           /- ...\n",
    "#  \t    |- Spoof\n",
    "#           /- LA_T_3756331.npz\n",
    "#           /- ...\n",
    "#    |- Val\n",
    "#       |- Genuine\n",
    "#           /- ...\n",
    "#  \t    |- Spoof\n",
    "#           /- ...\n",
    "#    |- Test\n",
    "#       |- Genuine\n",
    "#           /- ...\n",
    "#  \t    |- Spoof\n",
    "#           /- ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_path='/data/Feature/Train'\n",
    "protocol_file='/data/ASVspoof2019.LA.cm.train.trn.txt'\n",
    "genuine_file_path=os.path.join(base_path, 'Genuine')\n",
    "spoof_file_path=os.path.join(base_path, 'Spoof')\n",
    "with open(protocol_file) as f:\n",
    "    infos=f.readlines()\n",
    "mapping={info.split(' ')[1]:info.replace('\\n', '').split(' ') for info in infos}\n",
    "way_mapping={\n",
    "    'A01':1, \n",
    "    'A01':2, \n",
    "    'A03':3, \n",
    "    'A04':4, \n",
    "    'A05':5, \n",
    "    'A06':6, \n",
    "    'A07':7, \n",
    "    'A08':8, \n",
    "    'A09':9, \n",
    "    'A10':10, \n",
    "    'A11':11, \n",
    "    'A12':12, \n",
    "    'A13':13, \n",
    "    'A14':14, \n",
    "    'A15':15, \n",
    "    'A16':16, \n",
    "    'A17':17, \n",
    "}\n",
    "genuine_file_path=list(Path(genuine_file_path).glob('*.npz'))\n",
    "spoof_file_path=list(Path(spoof_file_path).glob('*.npz'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "counter_spoof_file=random.sample(list(spoof_file_path), k=len(list(genuine_file_path)))\n",
    "#choose the same number of LFCC features of spoofed utterances as that of genuine utterances to make a balanced dataset "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "67"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#to read .npz file into memory and find the shortest length of the LFCC features in our selected features\n",
    "spoof_list=[]\n",
    "for i in counter_spoof_file:\n",
    "    name=i.name.replace('.npz', '')\n",
    "    feature=np.load(i)['arr_0']\n",
    "    bundle=(name, feature)\n",
    "    spoof_list.append(bundle)\n",
    "\n",
    "genuine_list=[]\n",
    "for i in list(genuine_file_path):\n",
    "    name=i.name.replace('.npz', '')\n",
    "    feature=np.load(i)['arr_0']\n",
    "    bundle=(name, feature)\n",
    "    genuine_list.append(bundle)\n",
    "\n",
    "short=512\n",
    "for spoof in spoof_list:\n",
    "    length=spoof[1].shape[0]\n",
    "    if length < short:\n",
    "        short=length\n",
    "for genuine in genuine_list:\n",
    "    length=genuine[1].shape[0]\n",
    "    if length < short:\n",
    "        short=length\n",
    "short"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define the SCL, we directly copy this from https://keras.io/examples/vision/supervised-contrastive-learning/#1-supervised-contrastive-learning-loss-function\n",
    "class SupervisedContrastiveLoss(keras.losses.Loss):\n",
    "    def __init__(self, temperature=1, name=None):\n",
    "        super(SupervisedContrastiveLoss, self).__init__(name=name)\n",
    "        self.temperature = temperature\n",
    "\n",
    "    def __call__(self, labels, feature_vectors, sample_weight=None):\n",
    "        # Normalize feature vectors\n",
    "        feature_vectors_normalized = tf.math.l2_normalize(feature_vectors, axis=1)\n",
    "        # Compute logits\n",
    "        logits = tf.divide(\n",
    "            tf.matmul(\n",
    "                feature_vectors_normalized, tf.transpose(feature_vectors_normalized)\n",
    "            ),\n",
    "            self.temperature,\n",
    "        )\n",
    "        return tfa.losses.npairs_loss(tf.squeeze(labels), logits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss=SupervisedContrastiveLoss(temperature=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute SCL of datasets of different length\n",
    "def compute_loss(pca_dim, shortest):\n",
    "    los=[]\n",
    "    for i in range(int(pca_dim/60), shortest+1):\n",
    "        for_plot=[]\n",
    "\n",
    "        for bundle in spoof_list:\n",
    "            name=bundle[0]\n",
    "            feature=bundle[1]\n",
    "            frames=feature.shape[0]\n",
    "            index_start=0\n",
    "            using_feature=feature[index_start:index_start+i, :]\n",
    "            target=1\n",
    "            bundle=(using_feature, target)\n",
    "            for_plot.append(bundle)\n",
    "\n",
    "\n",
    "        for bundle in genuine_list:\n",
    "            name=bundle[0]\n",
    "            feature=bundle[1]\n",
    "            frames=feature.shape[0]\n",
    "            index_start=0\n",
    "            using_feature=feature[index_start:index_start+i, :]\n",
    "            target=0\n",
    "            bundle=(using_feature, target)\n",
    "            for_plot.append(bundle)\n",
    "\n",
    "        a, b=[], []\n",
    "        for j in range(len(for_plot)):\n",
    "            a.append(for_plot[j][0])\n",
    "            b.append(for_plot[j][1])\n",
    "        arr=np.array(a, dtype=np.float32).reshape(5160, -1)\n",
    "        new_a=PCA(n_components=pca_dim).fit_transform(arr)\n",
    "        tensor_a=tf.constant(new_a, dtype=tf.float32)\n",
    "        tenosr_b=tf.constant(b, dtype=tf.float32)\n",
    "        feature=tf.reshape(tensor_a, shape=(5160, -1))\n",
    "        value=loss(tenosr_b, feature).numpy()\n",
    "        los.append(value)\n",
    "    return np.array(los)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate this loss for many times to avoid possible numerical errors\n",
    "def sampling(pca_dim, shortest, samples=60):\n",
    "    total_loss=np.zeros(shape=(shortest+1-pca_dim//60, ), dtype=np.float32)\n",
    "    for _ in range(samples):\n",
    "        los=compute_loss(pca_dim, shortest)\n",
    "        total_loss+=los\n",
    "    return total_loss/samples\n",
    "\n",
    "# sampling(180, 65, 2).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f864111b1d0>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAnvElEQVR4nO3deXxV1bn/8c+TOSEzGYQwTyIqg0YQsVapWuu1Fay1dlCsVa7V9qfW3ta2t8Pt7WBvW61eW+uAihWnOhSvVq0XccABDfMsYRQISRgzAAlJnt8fZ8ONmJAEQk6S/X2/XueVvdc+a5/n8Dqc5+y11l7L3B0REQmfmGgHICIi0aEEICISUkoAIiIhpQQgIhJSSgAiIiEVF+0A2iInJ8cHDBgQ7TBERLqUefPmbXP33EPLu1QCGDBgAEVFRdEOQ0SkSzGzDU2VqwlIRCSklABEREJKCUBEJKSUAEREQkoJQEQkpJQARERCSglARCSkQpEAXltZyp9fL452GCIinUooEsBbq7fx59lroh2GiEinEooEkJ+eRFVNHVU1ddEORUSk02gxAZhZXzObbWbLzWyZmd0YlH8p2G8ws8LD1F9vZkvMbKGZFTUqzzazV81sdfA3q33e0iflpycCUFax71i9hIhIl9OaK4A64BZ3HwGcDtxgZiOApcAlwJutOMc57j7a3RsniluBWe4+FJgV7B8TeWlJAJRV1hyrlxAR6XJaTADuXuLu84PtSmAFUODuK9x91VG89sXA9GB7OjDpKM51WAeuAEp1BSAiclCb+gDMbAAwBpjbhmoO/NPM5pnZ1Ebl+e5eEmxvBfKbec2pZlZkZkXl5eVtCfegvPTgCqBCVwAiIge0OgGYWSrwDHCTu1e04TXOdPdTgM8RaT4669AnuLsTSRSf4O73uXuhuxfm5n5iOutWSUuMIzk+VlcAIiKNtCoBmFk8kS//Ge7+bFtewN03B3/LgOeAscGhUjPrFZy/F1DWlvO2hZmRl56oPgARkUZaMwrIgGnACne/vS0nN7MeZpZ2YBs4n0jnMcDzwJRgewowsy3nbqv8tCRdAYiINNKaK4AJwBXAxGAo50Izu9DMJpvZJmA88KKZvQJgZr3N7B9B3XxgjpktAt4HXnT3l4NjtwHnmdlq4Nxg/5jRFYCIyMe1uCSku88BrJnDzzXx/C3AhcH2WmBUM+fdDnym1ZEepby0JGZXHLNWJhGRLicUdwJDZChodW297gYWEQmEKAFEhoKqH0BEJCI0CSBPN4OJiHxMeBJAMB1EuTqCRUSAECUATQchIvJxoUkAqYlxpCTEUqrpIEREgBAlADMjL033AoiIHBCaBACRSeHUBCQiEhGqBJCfnqRFYUREAuFKAGmJlFbUEJl8VEQk3EKVAPLSE9m7X3cDi4hAyBLA/90NrI5gEZFQJYCDawOrH0BEJFwJ4MDNYBoKKiISsgSQpwnhREQOClUCSE2Mo4fuBhYRAUKWACDSEVxaqSsAEZHQJYDctETKdQUgIhK+BKArABGRiBYTgJn1NbPZZrbczJaZ2Y1B+ZeC/QYzK2xL3eDYz81sc+OF5tvvbTUvPz2R0op9uhtYREKvxUXhgTrgFnefb2ZpwDwzexVYClwC3NvWuu6+PDh+h7v//mjeQFvlpyexb38DlTV1pCfFd+RLi4h0Ki1eAbh7ibvPD7YrgRVAgbuvcPdVR1L36MM+crlpwb0AGgoqIiHXpj4AMxsAjAHmtvWFmqn7bTNbbGYPmllWM/WmmlmRmRWVl5e39WU/QdNBiIhEtDoBmFkq8Axwk7tXtOVFmql7DzAYGA2UAH9oqq673+fuhe5emJub25aXbdKBBFCmjmARCblWJQAziyfyBT7D3Z9tyws0V9fdS9293t0bgPuBsW0575HKSzuwNrCuAEQk3FozCsiAacAKd7+9LSc/XF0z69VodzKRTuVjrkdiHKmJcZoOQkRCrzVXABOAK4CJjYdsmtlkM9sEjAdeNLNXAMyst5n943B1g2P/ZWZLzGwxcA5wc7u+s8PIS0+kTFcAIhJyLQ4Ddfc5gDVz+Lkmnr8FuLCluu5+RevDbF/5aUnqAxCR0AvdncAQuQJQH4CIhF0oE0B+epLuBhaR0AtlAshLS6SmroGKfVobWETCK5QJ4OC9ABoJJCIhFsoEoHsBRERCmgDytTSkiEg4E0CeFocXEQlnAkhJiCMvLZGVW9s0pZGISLcSygQAMH5wT95Zs11DQUUktEKbAM4Y3JPyyhpWl1VFOxQRkagIcQLIAeCd4m1RjkREJDpCmwD6ZqfQLzuFt9dsj3YoIiJREdoEAJFmoPfWbqeuviHaoYiIdLhwJ4AhOVTuq2PpFo0GEpHwCXcCGNwTgHfWqB9ARMIn1AkgJzWR4cel8U6x+gFEJHxCnQAgcj/AB+t3sG9/fbRDERHpUKFPABMG51BT18D8jTujHYqISIcKfQIYNyib2BjjXQ0HFZGQaTEBmFlfM5ttZsvNbJmZ3RiUfynYbzCzwsPUv8DMVplZsZnd2qh8oJnNDcqfNLOE9nlLbZOWFM/JBRm8rRvCRCRkWnMFUAfc4u4jgNOBG8xsBLAUuAR4s7mKZhYL/An4HDAC+EpQF+C3wB3uPgTYCXzziN/FUZowpCeLNu2mct/+aIUgItLhWkwA7l7i7vOD7UpgBVDg7ivcfVUL1ccCxe6+1t1rgSeAi83MgInA08HzpgOTjvA9HLUJg3Oob3DeX7cjWiGIiHS4NvUBmNkAYAwwt5VVCoCPGu1vCsp6Arvcve6Q8qZec6qZFZlZUXl5eVvCbbVT+meREBfD2xoOKiIh0uoEYGapwDPATe7eYbfOuvt97l7o7oW5ubnH5DWS4mMp7J+lG8JEJFRalQDMLJ7Il/8Md3+2DeffDPRttN8nKNsOZJpZ3CHlUTNhSA4rt1ayflt1NMMQEekwrRkFZMA0YIW7397G838ADA1G/CQAlwPPe2QVltnApcHzpgAz23judjV5TAFpSXF896mFmhxOREKhNVcAE4ArgIlmtjB4XGhmk81sEzAeeNHMXgEws95m9g+AoI3/28ArRDqPn3L3ZcF5fwB818yKifQJTGvXd9ZGvTOT+dXkk5m/cRd3vVYczVBERDqEdaUlEQsLC72oqOiYvsYtTy3iuQWbeGLqeMYOzD6mryUi0hHMbJ67f+J+rdDfCXyo/7j4RPpmp3DzkwvZvVf3BYhI96UEcIjUxDjuvHwMpRX7+PFzS7RovIh0W0oATRjdN5ObzxvGC4tLeOz9jdEOR0TkmFACaMZ1nx7MmUNy+PFzS/n+04s0TYSIdDtKAM2IjTGmXVXI9WcP5ul5m7jgj29pxlAR6VaUAA4jMS6W718wnL9ddwbxscZX7n+PX/zPcjZsr1bfgIh0eRoG2kp7auu47aWVPPLuBiCynGRh/yxO7Z/F2IHZnFSQQWyMRSU2EZHDaW4YqBJAGxWXVfHe2u3M27CTeRt2snHHHgCyUuI5a1gunx6Wy1nDcslJTYxqnCIiBygBHCNlFft4d+123lhVzhsflrO9uhaAs4/P5TsTh3Jq/6woRygiYacE0AEaGpxlWyp4dUUpf313PTv37OdTQ3P4zsShuqtYRKJGCaCDVdfUMWPuBu57cy3bqmoZP6gnf7hsFL0zk6MdmoiEjKaC6GA9EuOYetZg3vr+RH560QiWbt7Npfe8Q3FZZbRDExEBlACOueSEWK4+cyBP/Ovp1NY7l/7lXRZs3BntsERElAA6yom9M3jmW+PJSI7nq/fP5fVVZdEOSURCTgmgA/Xv2YOnrzuDgTk9uGZ6ETMXRnURNBEJOSWADpablsiT/3o6hQOyuPnJhTy/aEu0QxKRkFICiIK0pHgeumoshQOyufnJhby8dGu0QxKREFICiJLkhFgevOo0RvXJ4DuPz+e1laXRDklEQkYJIIpSE+N4+OqxnNArnesenc9bq8ujHZKIhEiLCcDM+prZbDNbbmbLzOzGoDzbzF41s9XB30/MeWBm5zRaSH6hme0zs0nBsYfNbF2jY6Pb+811BelJ8Txy9VgG5fTg2keKmLtWU06LSMdozRVAHXCLu48ATgduMLMRwK3ALHcfCswK9j/G3We7+2h3Hw1MBPYA/2z0lH87cNzdFx7dW+m6MlMSmHHNOPpkpXDN9CKWbt4d7ZBEJARaTADuXuLu84PtSmAFUABcDEwPnjYdmNTCqS4FXnL3PUccbTfWMzWRv35zLOnJ8Vz54PsUl1VFOyQR6eba1AdgZgOAMcBcIN/dS4JDW4H8FqpfDjx+SNmvzGyxmd1hZk3On2xmU82syMyKysu7dxt5r4xkHr1mHDFmXDFtLpt2KleKyLHT6gRgZqnAM8BN7l7R+JhHZpRrdlY5M+sFnAy80qj4h8Bw4DQgG/hBU3Xd/T53L3T3wtzc3NaG22UNzOnBI1ePpaqmjiumvU95ZU20QxKRbqpVCcDM4ol8+c9w92eD4tLgi/3AF/zh5ja4DHjO3Q+urB40Lbm71wAPAWOP5A10RyN6p/PQVadRsnsvV0yby0c7dCUgIu2vNaOADJgGrHD32xsdeh6YEmxPAWYe5jRf4ZDmn0bJw4j0HyxtddQhUDggm/uvLGTzrr1ceOdbumNYRNpda64AJgBXABMbDdm8ELgNOM/MVgPnBvuYWaGZPXCgctBv0Bd445DzzjCzJcASIAf45dG+me7mU0Nz+cf/+xRD81P5f48v4N/+tog9tXXRDktEugktCNMF1NU3cOes1dw9u5iBOT34w5dGMaaflpoUkdbRgjBdWFxsDLecfzwzrhlHdU0dk//8DtfPmMeacg0VFZEjpwTQhZwxOIdZt5zNTecO5Y1V5Zx/x5v88NnFlOzeG+3QRKQLUgLoYlIT47jp3GG88f1zuHJ8f56et4mzf/c6z87fFO3QRKSLUQLoonJSE/nZ50/ktVvO5pR+WXz3qUX8/pVVNDR0nT4dEYkuJYAurm92CtOvHsuXC/ty9+xivvP4AvbW1kc7LBHpAuKiHYAcvYS4GG774skMyUvl1y+tYNPOPdx/ZSF56UnRDk1EOjFdAXQTZsa1Zw3i3q+fyoelVXzuzrf463sb2F/fEO3QRKSTUgLoZs4/8Tieu+EMBuel8pO/L+WCP77Jq8tL6Ur3e4hIx9CNYN2Uu/O/K8r4zUsrWFtezbiB2UweU8CI3ukMy08jKT422iGKSAdp7kYwJYBubn99A0988BF3/u9qtlVFZhaNMRiUm8qIXumM7JPB6L6ZnNg7g+QEJQWR7kgJIOQaGpyNO/awcmsFy0sqWVFSwbLNu9myex8AsTHGsPw0Ti5IZ0SvdEb0zmB4rzTSk+KjHLmIHC0lAGlSWeU+Fn+0m8WbdrHgo10s31LB9urag8cLMpPJS08kMzmerJQEMlMSGJTbg0ljCkhN1CAyka5ACUBaxd0pq6xheUkFy7dU8GFpJTuqa9m5p5ad1fvZtaeW6tp60pLi+Oq4flx1xgB6ZSRHO2wROQwlAGk3Czbu5IE563hpSQkxZlw0shdTzxrMiN7p0Q5NRJqgBCDt7qMde3jo7fU8+cFGqmvrOfv4XL716cGMHZhNZJ0fEekMlADkmNm9dz+PvreBB+esY3t1Laf2z+LaTw1k/OAcMpLViSwSbUoAcszt21/PU0Ufce8ba9m8KzJF9cCcHozsk8GoPpmc2j+LkwoyiI3R1YFIR1ICkA6zv76BuWt3sGjTLhZ9tIvFm3aztSIy3DQ9KY7TB/VkwpAcJgzJYUheapSjFen+mksALY7jM7O+wCNAPuDAfe5+p5llA08CA4D1wGXuvrOJ+vVE1v0F2OjuXwjKBwJPAD2BecAV7l57aH3peuJjYzhzaA5nDs05WLZ19z7eX7+Dd4q3Mad4G/9cXgpAYf8svnX2YCYOz1O/gUgHa/EKwMx6Ab3cfb6ZpRH5sp4EXAXscPfbzOxWIMvdf9BE/Sp3/8TPPDN7CnjW3Z8ws78Ai9z9nsPFoiuA7mPj9j28uqKUB+esY/OuvQw/Lo3rPj2Yi0b2Ii5WU1SJtKd2awIys5nA3cHjbHcvCZLE6+5+fBPP/0QCsMhPvXLgOHevM7PxwM/d/bOHe20lgO5nf30D/7NoC/e8vobVZVX0ykjicyf14oKTjuPU/lnqLxBpB0fcBHTISQYAY4C5QL67lwSHthJpImpKkpkVAXXAbe7+dyLNPrvcvS54ziagoC2xSPcQHxvDJaf0YdLoAmatLOPJDzby6NwNPPj2OnJSE/nsiflMGlNAYf8sNRGJtLNWJwAzSwWeAW5y94rG/xnd3c2suUuJ/u6+2cwGAa+Z2RJgdxtedyowFaBfv36trSZdTEyMcd6IfM4bkU9VTR2zV5bx8tKtPLdgMzPmbmRYfipfG9efyacUaH4ikXbSqiYgM4sHXgBecffbg7JVtKIJ6JDzPByc5xnUBCStsKe2jhcWlTBj7gYWbdpNcnwsXxjVm2vPGqQRRCKt1FwTUIu9bUF7/TRgxYEv/8DzwJRgewows4m6WWaWGGznABOA5R7JOrOBSw9XXyQlIY7LTuvLzG+fyQvfOZNJYwp4ftEWzrvjDW54bD4rt1ZEO0SRLqs1o4DOBN4iMpTzwPqCPyLSD/AU0A/YQGQY6A4zKwSuc/drzOwM4N6gXgzwR3efFpx3EJFhoNnAAuDr7l5zuFh0BSAA26tqmDZnHY+8u4GqmjrOH5HPzecN44RemotIpCm6EUy6nV17annw7fU89PY6auoauP2yUVw0sne0wxLpdI64CUiks8pMSeC75w3j9e+dzciCDL792AL+NLtY6x+LtJISgHR5PVMTefSacXxhVG9+98oqfvDMYmrrGlquKBJyWtJJuoWk+FjuvHw0A3qmcNdrxWzauZd7vn6qZiMVOQxdAUi3YWZ89/zj+cOXRvHB+h3cMGM+dfW6EhBpjhKAdDtfPLUPv5p0MnOKt/HLF1dEOxyRTktNQNItXXZaX1ZureTBt9cx/Lg0Lh+ru8hFDqUrAOm2fnThcD41NIefzFzKB+t3RDsckU5HCUC6rbjYGO7+yikUZCZz3V/nHVylTEQilACkW8tIieeBKYXU1jUw9ZEidQqLNKIEIN3ekLw0bvviSJZtqWDmwi3RDkek01ACkFC48OTjOKFXOn+aXUx9g+4UFgElAAkJM+M7E4ewdls1Ly4pabmCSAgoAUhoXHDicQzJS+Xu11bToKsAESUACY+YGOPb5wzhw9Iq/rl8a7TDEYk6JQAJlYtG9mJAzxT++zXNGiqiBCChEhcbw/XnDGHZlgpmryqLdjgiUaUEIKEzeUwBfbKSuWuWrgIk3JQAJHTiY2P41tmDWfjRLuYUb4t2OCJRowQgoXTpqX3olZHEf76wnH3766MdjkhUtJgAzKyvmc02s+VmtszMbgzKs83sVTNbHfzNaqLuaDN7N6i32My+3OjYw2a2zswWBo/R7frORA4jMS6W2744kg9Lq/jli8ujHY5IVLTmCqAOuMXdRwCnAzeY2QjgVmCWuw8FZgX7h9oDXOnuJwIXAH80s8xGx//N3UcHj4VH8T5E2uzTw3KZetYgHn1vIy8v1bBQCZ8WE4C7l7j7/GC7ElgBFAAXA9ODp00HJjVR90N3Xx1sbwHKgNx2iVykHXzv/OMZ2SeDHzyzWLOFSui0qQ/AzAYAY4C5QL67H7infiuQ30LdsUACsKZR8a+CpqE7zCyxmXpTzazIzIrKy8vbEq5IixLiYrjr8jHU1Tdw8xMLNVuohEqrE4CZpQLPADe5e0XjYx4ZS9fseDoz6wX8FfiGux/4H/ZDYDhwGpAN/KCpuu5+n7sXunthbq4uHqT9DcjpwS8nn8T763dw9+ziaIcj0mFalQDMLJ7Il/8Md382KC4NvtgPfME3eVeNmaUDLwI/dvf3DpQHTUvu7jXAQ8DYI38bIkdn8pg+XDKmgLtmreafy9QfIOHQmlFABkwDVrj77Y0OPQ9MCbanADObqJsAPAc84u5PH3LsQPIwIv0HS48gfpF284tJJ3Fyn0yunzGfl5dqxlDp/lpzBTABuAKY2GjI5oXAbcB5ZrYaODfYx8wKzeyBoO5lwFnAVU0M95xhZkuAJUAO8Mt2e1ciRyA1MY6/fnMsI/tkcMNjC3hhsRaPke7NutKt8IWFhV5UVBTtMKSbq6qp4xsPvc+8DTu548ujuXh0QbRDEjkqZjbP3QsPLdedwCKHSE2M4+FvjGXswGxufnIhz8zbFO2QRI4JJQCRJvRIjOOhq8YyfnBPvvf0Iv6+YHO0QxJpd0oAIs1ITojlgStP4/SBPfnuUwv5n0XqE5DuRQlA5DCSE2KZdlUhhQOyuenJhbyk9YSlG1ECEGlBSkIcD151GqP7ZvKdxxfoPgHpNpQARFoh0jF8GicVZHDDY/N5dXlptEMSOWpKACKtlJYUz/SrxzKidwbXPTqPmQvVMSxdmxKASBtkJMcz45pxFPbP4qYnF/LY3I3RDknkiCkBiLRRamIc068eyznH5/Gj55Zw7xtrWq4k0gkpAYgcgaT4WO694lQuGtmL37y0kt+/skoLzEuXExftAES6qvjYGO68fAypiXHcPbuYvfvr+fd/OYHI/IYinZ8SgMhRiI0xfnPJySTFxzJtzjpq6ur5xRdOIiZGSUA6PyUAkaNkZvzs8yNIjIvh3jfXUlvXwG8uGUmskoB0ckoAIu3AzLj1c8NJjIvhrteKqa1r4PdfGkVcrLrZpPNSAhBpJ2bGd88/nsT4WH73yioaHP745dFqDpJOSwlApJ3dcM4QAH73yiry0hL594tGRDkikaYpAYgcA9efPZjyyhoemLOO4zKSuOZTg6IdksgnKAGIHANmxk8uGkFpxT5++eIK8tOT+Pyo3tEOS+Rj1EMlcozExhh3fHk0Ywdkc8tTi3h3zfZohyTyMS0mADPra2azzWy5mS0zsxuD8mwze9XMVgd/s5qpPyV4zmozm9Ko/FQzW2JmxWZ2l+nuGemGkuJjuf/KQvr3TGHqI0Ws2loZ7ZBEDmrNFUAdcIu7jwBOB24wsxHArcAsdx8KzAr2P8bMsoGfAeOAscDPGiWKe4BrgaHB44KjfC8inVJGSjwPXz2W5IRYrn74A8oq9kU7JBGgFQnA3UvcfX6wXQmsAAqAi4HpwdOmA5OaqP5Z4FV33+HuO4FXgQvMrBeQ7u7veWQClUeaqS/SLRRkJvPgVaexc08t1zxSxJ7aumiHJNK2PgAzGwCMAeYC+e5+YH28rUB+E1UKgI8a7W8KygqC7UPLm3rNqWZWZGZF5eXlbQlXpFM5qSCD//7KGJZu3s2NTyykvkGTx0l0tToBmFkq8Axwk7tXND4W/Io/Jp9md7/P3QvdvTA3N/dYvIRIh/nMCfn89KIRvLq8lF//Y0W0w5GQa1UCMLN4Il/+M9z92aC4NGjKIfhb1kTVzUDfRvt9grLNwfah5SLd3lUTBnLVGQOYNmcd099ZH+1wJMRaMwrIgGnACne/vdGh54EDo3qmADObqP4KcL6ZZQWdv+cDrwRNRxVmdnpw/iubqS/SLf3kohGce0I+P3t+GX9+vVhrCUhUtOYKYAJwBTDRzBYGjwuB24DzzGw1cG6wj5kVmtkDAO6+A/hP4IPg8YugDOB64AGgGFgDvNR+b0ukc4uNMe7+6hg+P6o3//XyKn703FLq6huiHZaEjHWlXx6FhYVeVFQU7TBE2k1Dg/P7f67iz6+v4dPDcvnT104hNVE36Ev7MrN57l54aLnuBBaJopgY4/sXDOfXk09mTvE2LvvLu5TqPgHpIPqpIdIJfHVcP3pnJnHDjPlcdu+7PHbt6RRkJkc7LGlnDQ3+ieGSdQ0NVNfUU7WvjqqayGP33v3srK5lx55adu6pZWd1Ldd+ahBD89PaNR4lAJFO4uzj83j0mnFcOe19vnzvuzx+7en0zU6JdljSSjuqa1lTXsXa8irWllezpryassp9kS/14Mt9T219m8+bEBdDdkoCk0YXtHsCUB+ASCezeNMuvv7AXNKS4nn82tPp11NJoLPZX9/AipIK5m/YyfyNu5i/cSebdu49eDwhNob+PVPonZlMWlIcqYmRR0piHHGHLBAUG2P0SIglNSme1MRYeiTGkZmcQFaPeLJSEkhJiOVop0prrg9ACUCkE1q6eTdfnzaX5PhYHr/2dAbk9Ih2SKFWVrmP+Rt2seCjnSzYsIvFm3exb39k1FZ+eiKn9s9idN9MhuanMSinB32yUjrVmtBKACJdzPItFXztgfdIiIth+tVjGX5cerRD6tYaGpzNu/ZSXFbFmvIq1pRXH2zS2VZVC0B8rHFi7wzG9MvklH5ZnNI/i94ZSUf9C/1YUwIQ6YJWba3kygfnUl1Tz5+/dgpnDdN0KO1hZ3UtS7fsZsnm3azaWnnwS//Ar3qA7B4JDMrpweDcVIbmpzKmXyYn9s4gKT42ipEfGSUAkS6qZPdevvHQB6wuq+JXk07i8rH9oh1Sl9LQ4KzcWsm7a7fzwbodLNm8m827/q+9viAzmSF5qR97DM5NJbtHQhSjbl/NJQCNAhLp5HplJPO368Zzw2MLuPXZJWzcsYfvnX88MZ2ojbmzqa6p44XFW3h9VTnvrd3Ozj37AeibncyYfplcMb4/JxdkcFLvDDJS4qMcbfQoAYh0AWlJ8UybUshPZy7jz6+vYU15Fb+5ZGS3+pXaHtaWV/HX9zbwdNEmKmvqKMhMZuLwfMYP7sn4wT11b8UhlABEuoj42Bh+PfkkBuf24Lcvr+T8O97gV5NP5rMnHhft0KKqtGIfc1ZvY+aiLbz5YTnxscaFJ/fiyvEDOKVfZqfvoI0m9QGIdEErt1Zwy1OLWLalgkmje/PzL5xIZko4rgbKK2tYsnkXc1ZvZ05xOR+WVgFwXHoSXx3Xj8vH9iUvLSnKUXYu6gQW6Wb21zfwp9nF3P1aMdk9EvjtpSM55/i8aIfVrrbu3seCjTtZsnk3y0sqWLalgvLKGgAS42IYOzCbM4fkcObQHE44Ll39Is1QAhDpppZu3s0tTy1iVWklXx3Xjx9feAI9utiMou5OaUUNq8sqWVlSycKPInfXluyOTIwXG2MMzUtlRO90TuydwYm90xndN7NLDsmMBo0CEummTirIYOa3J3D7qx9y/1trebt4G7dfNopT+2dHO7Rm7a2tZ07xNl5fVcbykgqKS6uorKk7eLxPVjKFA7IZ0zeTMf0yOaFXur7sjwFdAYh0I++t3c4tTy2iZPderj1rENefPYSM5M4xzHHLrr3MXlXGrBVlvF28jZq6BlIT4zi5IIMheZGbrYbkpTI0L43ctMRoh9utqAlIJCQq9+3nP19YzlNFm0hLiuMbEwbyzQkDO3y8+7799by3djtvfriNN1eXU1wW6aztm53MZ4bnc+4J+YwdmE1CnJYlOdaUAERCZtmW3fz3rGJeXraVtMQ4vjFhAJ8f1Zu+2SlH3Zzi7uyormX99mrWllezfns1Jbv3sa2qlu1VNWyrqmF7VS11DX6ws/bTw3I5a1guQ/NSNTSzgykBiITUipIK7pq1mpeWbgXADHpnJNO/Zwr9e/ZgSF4qw/IjTS/56Ymf+HJuaHDWbqtm6ebI3DlLNu9mZUkFFfv+r80+NsbIT0skJy2RnNREevZIID89idMGZjNuYLba76PsiBOAmT0IXASUuftJQdko4C9AKrAe+Jq7VxxS73jgyUZFg4CfuvsfzeznwLVAeXDsR+7+j5behBKAyJFbW17F4k27Wb+9mg3b97B+ezXrtlWzK5gmASAtMY6ctERq6xqoqaunZn8De/fXU9cQ+Z5IjIvhhF7pjOidzuDcVAbmpDAwJ5U+WcnEx6opp7M6mlFADwN3A480KnsA+J67v2FmVwP/BvykcSV3XwWMDl48FtgMPNfoKXe4++/b8B5E5CgMyk1lUG7qx8rcne3VtawuraK4rJIPS6vYtXc/iXExwSOWxPgYBub0ONhZqy/67qPFBODub5rZgEOKhwFvBtuvAq9wSAI4xGeANe6+4UiCFJFjw8zISY0024wf3DPa4UgHO9JUvgy4ONj+EtC3hedfDjx+SNm3zWyxmT1oZlnNVTSzqWZWZGZF5eXlzT1NRETa6EgTwNXA9WY2D0gDapt7opklAF8A/tao+B5gMJEmohLgD83Vd/f73L3Q3Qtzc7UYhohIezmiO4HdfSVwPoCZDQP+5TBP/xww391LG9U/uG1m9wMvHEkcIiJy5I7oCsDM8oK/McC/ExkR1JyvcEjzj5n1arQ7GVh6JHGIiMiRazEBmNnjwLvA8Wa2ycy+CXzFzD4EVgJbgIeC5/Y2s380qtsDOA949pDT/peZLTGzxcA5wM3t8m5ERKTVdCOYiEg319x9ABrQKyISUkoAIiIh1aWagMysHDjSm8lygG3tGE5HUdwdr6vGrrg7VleKu7+7f2IcfZdKAEfDzIqaagPr7BR3x+uqsSvujtVV425MTUAiIiGlBCAiElJhSgD3RTuAI6S4O15XjV1xd6yuGvdBoekDEBGRjwvTFYCIiDSiBCAiElKhSABmdoGZrTKzYjO7NdrxNCdYG6HMzJY2Kss2s1fNbHXwt9m1E6LFzPqa2WwzW25my8zsxqC8U8duZklm9r6ZLQri/o+gfKCZzQ0+L08GU5p3OmYWa2YLzOyFYL/Tx21m64N5wBaaWVFQ1qk/JwBmlmlmT5vZSjNbYWbju0LcLen2CSBYjvJPRKalHkFkIrsR0Y2qWQ8DFxxSdiswy92HArOC/c6mDrjF3UcApwM3BP/GnT32GmCiu48isjbFBWZ2OvBbIkuWDgF2At+MXoiHdSOwotF+V4n7HHcf3WgMfWf/nADcCbzs7sOBUUT+3btC3Ifn7t36AYwHXmm0/0Pgh9GO6zDxDgCWNtpfBfQKtnsBq6IdYyvew0wis8B2mdiBFGA+MI7I3Z1xTX1+OssD6EPkS2cikfU0rIvEvR7IOaSsU39OgAxgHcGgma4Sd2se3f4KACgAPmq0vyko6yry3b0k2N4K5EczmJYE60ePAebSBWIPmlEWAmVE1rdeA+xy97rgKZ318/JH4PtAQ7Dfk64RtwP/NLN5ZjY1KOvsn5OBQDnwUNDk9kAw1X1nj7tFYUgA3YZHfmp02nG7ZpYKPAPc5O4VjY911tjdvd7dRxP5RT0WGB7diFpmZhcBZe4+L9qxHIEz3f0UIk2yN5jZWY0PdtLPSRxwCnCPu48BqjmkuaeTxt2iMCSAzXx80fo+QVlXUXpgBbXgb1mU42mSmcUT+fKf4e4HFgDqErEDuPsuYDaRppNMMzuwXGpn/LxMAL5gZuuBJ4g0A91J548bd98c/C0DniOSdDv752QTsMnd5wb7TxNJCJ097haFIQF8AAwNRkgkAJcDz0c5prZ4HpgSbE8h0r7eqZiZAdOAFe5+e6NDnTp2M8s1s8xgO5lIv8UKIong0uBpnS5ud/+hu/dx9wFEPs+vufvX6ORxm1kPM0s7sE1kXfGldPLPibtvBT4ys+ODos8Ay+nkcbdKtDshOuIBXAh8SKR998fRjucwcT4OlAD7ifzq+CaRtt1ZwGrgf4HsaMfZRNxnErn8XQwsDB4XdvbYgZHAgiDupcBPg/JBwPtAMfA3IDHasR7mPZwNvNAV4g7iWxQ8lh34v9jZPydBjKOBouCz8ncgqyvE3dJDU0GIiIRUGJqARESkCUoAIiIhpQQgIhJSSgAiIiGlBCAiElJKACIiIaUEICISUv8fZyGDPXi2OyIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#samples=1 for smoke-test\n",
    "a=sampling(60, short, 1)\n",
    "plt.plot(a)\n",
    "# array([21.522217, 21.123184, 20.978636, 20.910416, 20.87409 , 20.855097,\n",
    "#     20.844511, 20.826216, 20.794949, 20.757719, 20.715427, 20.683243,\n",
    "#     20.654678, 20.609268, 20.577244, 20.56762 , 20.559587, 20.554028,\n",
    "#     20.548992, 20.532965, 20.515858, 20.511631, 20.508186, 20.494783,\n",
    "#     20.481657, 20.469555, 20.449566, 20.415064, 20.366112, 20.29015 ,\n",
    "#     20.241735, 20.206081, 20.164122, 20.121264, 20.07475 , 20.026148,\n",
    "#     19.974396, 19.933092, 19.890413, 19.852352, 19.822243, 19.795166,\n",
    "#     19.78155 , 19.775568, 19.77166 , 19.771927, 19.778786, 19.78453 ,\n",
    "#     19.789854, 19.801006, 19.814644, 19.822136, 19.831566, 19.843311,\n",
    "#     19.858652, 19.87511 , 19.885893, 19.894947, 19.902826, 19.908686,\n",
    "#     19.917107, 19.926258, 19.934507, 19.936892, 19.936365],\n",
    "#     dtype=float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b=sampling(120, short, 100)\n",
    "# array([21.097296, 20.928377, 20.817253, 20.744236, 20.700178, 20.673655,\n",
    "#        20.64334 , 20.605518, 20.56773 , 20.527596, 20.494265, 20.462645,\n",
    "#        20.42074 , 20.386883, 20.373096, 20.360886, 20.3515  , 20.34234 ,\n",
    "#        20.32536 , 20.306707, 20.30042 , 20.292835, 20.278442, 20.265417,\n",
    "#        20.250805, 20.229948, 20.193758, 20.144403, 20.067535, 20.016518,\n",
    "#        19.977465, 19.933535, 19.8881  , 19.837147, 19.782572, 19.72631 ,\n",
    "#        19.680334, 19.63247 , 19.589724, 19.555426, 19.523794, 19.502132,\n",
    "#        19.490362, 19.482433, 19.481913, 19.487122, 19.49326 , 19.50129 ,\n",
    "#        19.512474, 19.52678 , 19.536047, 19.548637, 19.563725, 19.581848,\n",
    "#        19.600828, 19.615416, 19.628656, 19.641958, 19.653297, 19.667889,\n",
    "#        19.683592, 19.694887, 19.703693, 19.708015], dtype=float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#we choose 60 here to reduce the computation time\n",
    "sampling(180, short, 60)\n",
    "# array([20.928394, 20.817251, 20.739382, 20.687216, 20.649586, 20.609169,\n",
    "#        20.56354 , 20.519043, 20.474443, 20.438   , 20.404505, 20.361128,\n",
    "#        20.326286, 20.310675, 20.295412, 20.284178, 20.272837, 20.254536,\n",
    "#        20.235962, 20.228611, 20.219383, 20.203674, 20.18978 , 20.174437,\n",
    "#        20.152134, 20.115072, 20.064825, 19.987852, 19.936275, 19.896858,\n",
    "#        19.852089, 19.806501, 19.75445 , 19.699453, 19.643366, 19.59709 ,\n",
    "#        19.54877 , 19.50431 , 19.468105, 19.43395 , 19.408   , 19.392189,\n",
    "#        19.380388, 19.376055, 19.376692, 19.379065, 19.38375 , 19.392857,\n",
    "#        19.40464 , 19.412865, 19.424658, 19.43964 , 19.457191, 19.475262,\n",
    "#        19.490425, 19.504164, 19.51871 , 19.531963, 19.547993, 19.5651  ,\n",
    "#        19.577978, 19.588915, 19.59575 ], dtype=float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
  },
  "kernelspec": {
   "display_name": "Python 3.6.9 64-bit",
   "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.6.9"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}