20200319-picard.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gcsfs\n",
    "from scipy.io import loadmat\n",
    "from scipy.signal import resample\n",
    "from scipy.signal import decimate\n",
    "from scipy.stats import zscore\n",
    "import matplotlib.pyplot as plt\n",
    "import os.path as op\n",
    "import numpy as np\n",
    "from neurodsp.filt.filter import filter_signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = gcsfs.GCSFileSystem(project='learning-2-learn-221016')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = 'learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-1/sub-SA_sess-20180918_chan-1.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-2/sub-SA_sess-20180918_chan-2.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-3/sub-SA_sess-20180918_chan-3.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-4/sub-SA_sess-20180918_chan-4.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-5/sub-SA_sess-20180918_chan-5.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-6/sub-SA_sess-20180918_chan-6.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-7/sub-SA_sess-20180918_chan-7.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-8/sub-SA_sess-20180918_chan-8.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-9/sub-SA_sess-20180918_chan-9.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-10/sub-SA_sess-20180918_chan-10.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-11/sub-SA_sess-20180918_chan-11.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-12/sub-SA_sess-20180918_chan-12.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-13/sub-SA_sess-20180918_chan-13.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-14/sub-SA_sess-20180918_chan-14.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-15/sub-SA_sess-20180918_chan-15.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-16/sub-SA_sess-20180918_chan-16.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-17/sub-SA_sess-20180918_chan-17.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-18/sub-SA_sess-20180918_chan-18.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-19/sub-SA_sess-20180918_chan-19.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-20/sub-SA_sess-20180918_chan-20.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-21/sub-SA_sess-20180918_chan-21.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-22/sub-SA_sess-20180918_chan-22.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-23/sub-SA_sess-20180918_chan-23.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-24/sub-SA_sess-20180918_chan-24.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-25/sub-SA_sess-20180918_chan-25.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-26/sub-SA_sess-20180918_chan-26.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-27/sub-SA_sess-20180918_chan-27.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-28/sub-SA_sess-20180918_chan-28.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-29/sub-SA_sess-20180918_chan-29.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-30/sub-SA_sess-20180918_chan-30.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-31/sub-SA_sess-20180918_chan-31.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-32/sub-SA_sess-20180918_chan-32.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-33/sub-SA_sess-20180918_chan-33.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-34/sub-SA_sess-20180918_chan-34.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-35/sub-SA_sess-20180918_chan-35.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-36/sub-SA_sess-20180918_chan-36.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-37/sub-SA_sess-20180918_chan-37.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-38/sub-SA_sess-20180918_chan-38.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-39/sub-SA_sess-20180918_chan-39.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-40/sub-SA_sess-20180918_chan-40.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-41/sub-SA_sess-20180918_chan-41.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-42/sub-SA_sess-20180918_chan-42.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-43/sub-SA_sess-20180918_chan-43.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-44/sub-SA_sess-20180918_chan-44.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-45/sub-SA_sess-20180918_chan-45.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-46/sub-SA_sess-20180918_chan-46.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-47/sub-SA_sess-20180918_chan-47.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-48/sub-SA_sess-20180918_chan-48.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-49/sub-SA_sess-20180918_chan-49.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-50/sub-SA_sess-20180918_chan-50.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-51/sub-SA_sess-20180918_chan-51.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-52/sub-SA_sess-20180918_chan-52.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-53/sub-SA_sess-20180918_chan-53.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-54/sub-SA_sess-20180918_chan-54.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-55/sub-SA_sess-20180918_chan-55.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-56/sub-SA_sess-20180918_chan-56.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-57/sub-SA_sess-20180918_chan-57.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-58/sub-SA_sess-20180918_chan-58.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-59/sub-SA_sess-20180918_chan-59.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-60/sub-SA_sess-20180918_chan-60.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-61/sub-SA_sess-20180918_chan-61.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-62/sub-SA_sess-20180918_chan-62.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-63/sub-SA_sess-20180918_chan-63.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-64/sub-SA_sess-20180918_chan-64.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-65/sub-SA_sess-20180918_chan-65.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-66/sub-SA_sess-20180918_chan-66.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-67/sub-SA_sess-20180918_chan-67.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-68/sub-SA_sess-20180918_chan-68.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-69/sub-SA_sess-20180918_chan-69.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-70/sub-SA_sess-20180918_chan-70.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-71/sub-SA_sess-20180918_chan-71.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-72/sub-SA_sess-20180918_chan-72.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-73/sub-SA_sess-20180918_chan-73.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-74/sub-SA_sess-20180918_chan-74.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-75/sub-SA_sess-20180918_chan-75.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-76/sub-SA_sess-20180918_chan-76.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-77/sub-SA_sess-20180918_chan-77.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-78/sub-SA_sess-20180918_chan-78.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-79/sub-SA_sess-20180918_chan-79.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-80/sub-SA_sess-20180918_chan-80.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-81/sub-SA_sess-20180918_chan-81.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-82/sub-SA_sess-20180918_chan-82.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-83/sub-SA_sess-20180918_chan-83.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-84/sub-SA_sess-20180918_chan-84.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-85/sub-SA_sess-20180918_chan-85.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-86/sub-SA_sess-20180918_chan-86.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-87/sub-SA_sess-20180918_chan-87.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-88/sub-SA_sess-20180918_chan-88.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-89/sub-SA_sess-20180918_chan-89.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-90/sub-SA_sess-20180918_chan-90.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-91/sub-SA_sess-20180918_chan-91.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-92/sub-SA_sess-20180918_chan-92.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-93/sub-SA_sess-20180918_chan-93.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-94/sub-SA_sess-20180918_chan-94.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-95/sub-SA_sess-20180918_chan-95.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-96/sub-SA_sess-20180918_chan-96.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-97/sub-SA_sess-20180918_chan-97.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-98/sub-SA_sess-20180918_chan-98.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-99/sub-SA_sess-20180918_chan-99.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-100/sub-SA_sess-20180918_chan-100.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-101/sub-SA_sess-20180918_chan-101.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-102/sub-SA_sess-20180918_chan-102.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-103/sub-SA_sess-20180918_chan-103.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-104/sub-SA_sess-20180918_chan-104.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-105/sub-SA_sess-20180918_chan-105.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-106/sub-SA_sess-20180918_chan-106.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-107/sub-SA_sess-20180918_chan-107.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-108/sub-SA_sess-20180918_chan-108.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-109/sub-SA_sess-20180918_chan-109.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-110/sub-SA_sess-20180918_chan-110.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-111/sub-SA_sess-20180918_chan-111.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-112/sub-SA_sess-20180918_chan-112.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-113/sub-SA_sess-20180918_chan-113.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-114/sub-SA_sess-20180918_chan-114.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-115/sub-SA_sess-20180918_chan-115.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-116/sub-SA_sess-20180918_chan-116.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-117/sub-SA_sess-20180918_chan-117.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-118/sub-SA_sess-20180918_chan-118.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-119/sub-SA_sess-20180918_chan-119.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-120/sub-SA_sess-20180918_chan-120.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-121/sub-SA_sess-20180918_chan-121.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-122/sub-SA_sess-20180918_chan-122.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-123/sub-SA_sess-20180918_chan-123.mat\n",
      "Opening learning2learn/Buffalo/wcst-nhp/SA/sess-20180918/ephys/chan-124/sub-SA_sess-20180918_chan-124.mat\n"
     ]
    }
   ],
   "source": [
    "data = np.empty((124, 11796177))\n",
    "\n",
    "for chan in range(1, 125):\n",
    "    fname = op.join(path, f'chan-{chan}', f'sub-SA_sess-20180918_chan-{chan}.mat')\n",
    "    print(f\"Opening {fname}\")\n",
    "    with fs.open(fname) as f_chan:\n",
    "        f_chan = loadmat(f_chan)\n",
    "        data[chan-1] = f_chan['data'].squeeze()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# debug"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "resampled_data = resample(data, data.shape[-1] // 5, axis=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('resampled', resampled_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "decimated_data = decimate(data, 5, axis=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('decimated', decimated_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from picard import picard"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mixing(W, K):\n",
    "    w = np.dot(W, K)\n",
    "    return np.dot(w.T, np.linalg.inv(np.dot(w, w.T)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_components = 15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-25-4afba9e951d6>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mK\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mW\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mY\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpicard\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresampled_data\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_components\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mn_components\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/opt/tljh/user/lib/python3.7/site-packages/picard/solver.py\u001b[0m in \u001b[0;36mpicard\u001b[0;34m(X, fun, n_components, ortho, extended, whiten, return_X_mean, return_n_iter, centering, max_iter, tol, m, ls_tries, lambda_min, check_fun, w_init, fastica_it, random_state, verbose)\u001b[0m\n\u001b[1;32m    204\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    205\u001b[0m     \u001b[0;32mwith\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0merrstate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'raise'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# So code breaks if overflow/Nans\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 206\u001b[0;31m         \u001b[0mY\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mW\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minfos\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcore_picard\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    207\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    208\u001b[0m     \u001b[0;32mdel\u001b[0m \u001b[0mX1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.7/site-packages/picard/_core_picard.py\u001b[0m in \u001b[0;36mcore_picard\u001b[0;34m(X, density, ortho, extended, m, max_iter, tol, lambda_min, ls_tries, verbose, covariance)\u001b[0m\n\u001b[1;32m     86\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mn\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmax_iter\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     87\u001b[0m         \u001b[0;31m# Compute the score function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 88\u001b[0;31m         \u001b[0mpsiY\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpsidY\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdensity\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscore_and_der\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mY\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     89\u001b[0m         \u001b[0;31m# Compute the relative gradient and the Hessian off-diagonal\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     90\u001b[0m         \u001b[0mG\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minner\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpsiY\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mY\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mT\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.7/site-packages/picard/densities.py\u001b[0m in \u001b[0;36mscore_and_der\u001b[0;34m(self, Y)\u001b[0m\n\u001b[1;32m     46\u001b[0m         \u001b[0malpha\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0malpha\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     47\u001b[0m         \u001b[0mscore\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mne\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mevaluate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'tanh(alpha * Y)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 48\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mscore\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0malpha\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0malpha\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mscore\u001b[0m \u001b[0;34m**\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     49\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     50\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "K, W, Y = picard(resampled_data, n_components=n_components)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = mixing(W, K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.matshow(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(data[:, 200000:210000] + np.arange(0, data.shape[0] * 200, 200)[:, None])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = plt.plot((data[:20, 200000:200500] + np.arange(0, 20 * 5000, 5000)[:, None]).T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(Y[12, 200000:210000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_components = 3\n",
    "mean_corr = [0]\n",
    "break_it = False\n",
    "\n",
    "while not break_it:   \n",
    "    K1, W1, Y1 = picard(data[:, :data.shape[-1]//2], n_components=n_components)\n",
    "    K2, W2, Y2 = picard(data[:, data.shape[-1]//2:], n_components=n_components)\n",
    "    mixing1 = mixing(W1, K1)\n",
    "    mixing2 = mixing(W2, K2)\n",
    "    \n",
    "    corr = np.empty((mixing1.shape[-1], mixing2.shape[-1]))\n",
    "    for ii in range(corr.shape[0]):\n",
    "        for jj in range(corr.shape[1]):\n",
    "            corr[ii, jj] = np.abs(np.corrcoef(mixing1[:, ii], mixing2[:, jj])[0, 1])\n",
    "    \n",
    "    new_mixing2 = np.empty_like(mixing2)\n",
    "    max_corr = np.zeros(mixing1.shape[-1])\n",
    "    for it in range(mixing1.shape[-1]):\n",
    "        max_corr[it] = np.nanmax(corr)\n",
    "        idx = np.where(corr == max_corr[it])\n",
    "        new_mixing2[:, idx[0]] = mixing2[:, idx[1]]\n",
    "        corr[idx[0], :] = np.nan\n",
    "        corr[:, idx[1]] = np.nan\n",
    "   \n",
    "    this_mean_corr = np.mean(max_corr)\n",
    "    if this_mean_corr < mean_corr[-1] and (mean_corr[-1] - this_mean_corr) > 0.05 * mean_corr[-1]:\n",
    "        mean_corr.append(this_mean_corr)\n",
    "        break_it = True\n",
    "    else:\n",
    "        mean_corr.append(this_mean_corr)\n",
    "    n_components = n_components + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(mean_corr[1:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2)\n",
    "ax[0].matshow(new_mixing2[:10])\n",
    "ax[1].matshow(mixing1[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Another approach might be to rely on the fact that:\n",
    "# Y = WKX \n",
    "# To derive a cross-validated prediction of Y\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}