glmnet_cf_ensemble.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glmnet_python\n",
    "from glmnet import glmnet\n",
    "\n",
    "# Import relevant modules and setup for calling glmnet\n",
    "%matplotlib widget\n",
    "# %matplotlib inline\n",
    "\n",
    "import time\n",
    "import sys\n",
    "import os\n",
    "import re\n",
    "from itertools import compress\n",
    "from glob import glob\n",
    "import pickle\n",
    "\n",
    "import tensorflow as tf # data is in TFRecord format\n",
    "\n",
    "import scipy, importlib, pprint, matplotlib.pyplot as plt, warnings\n",
    "from scipy.io import loadmat\n",
    "\n",
    "import numpy as np\n",
    "from glmnet import glmnet; from glmnetPlot import glmnetPlot; from glmnetSet import glmnetSet\n",
    "from glmnetPrint import glmnetPrint; from glmnetCoef import glmnetCoef; from glmnetPredict import glmnetPredict\n",
    "from cvglmnet import cvglmnet; from cvglmnetCoef import cvglmnetCoef\n",
    "from cvglmnetPlot import cvglmnetPlot; from cvglmnetPredict import cvglmnetPredict\n",
    "\n",
    "# glmnet has many deprecation warnings from using scipy.* instad of numpy.*\n",
    "warnings.filterwarnings(\"ignore\", category=DeprecationWarning) \n",
    "\n",
    "def plot_distributions(y, y_hat):\n",
    "\n",
    "    fig, axs = plt.subplots(2, 4, sharey=False, sharex=False, tight_layout=True)\n",
    "    fig.set_size_inches(20,10)\n",
    "\n",
    "    for ic in range(0,8):\n",
    "        idx = np.unravel_index(ic, axs.shape)\n",
    "        if ic < 4:\n",
    "            axs[idx].set_title(\"Real %s\"%names[ic])\n",
    "            x = y[:,ic]        \n",
    "        else:\n",
    "            axs[idx].set_title(\"Predicted %s\"%names[ic-4])\n",
    "            x = y_hat[:,ic-4]\n",
    "        axs[idx].hist(x, bins=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define constants\n",
    "ncores = 28 # 56\n",
    "names = ['DA', '5HT', 'pH', 'NE']\n",
    "results_prefix = '/mnt/nfs/proj/in-vitro/Leonardo/glmnet/fits/'\n",
    "\n",
    "if not (os.path.exists(results_prefix)):\n",
    "    os.makedirs(results_prefix, exist_ok=True)\n",
    "\n",
    "data_prefix = '/mnt/nfs/proj/in-vitro/Mark/four_analyte/slow/allin'\n",
    "# good_probes = ['CF025']\n",
    "# good_probes = ['CF025', 'CF027']\n",
    "good_probes = ['CF025', 'CF027', 'CF057', 'CF064', 'CF066', 'CF078', 'CF081', 'CF082']\n",
    "n_records_per_label_per_probe = 1 # Loads N sweeps per concentration per probe\n",
    "# n_records_per_label_per_probe = -1 # Loads all records\n",
    "val_ratio = .1 # save some random percentage of probes different from holdout for validation \n",
    "# val_ratio = .0 # validation set will be empty\n",
    "hold_probe = 0\n",
    "# hold_probe = -1 # turn the validation set in holdout set (validation set will be empty)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2438, 999) (2438, 4)\n",
      "(387, 999) (387, 4)\n",
      "(271, 999) (271, 4)\n",
      "--- 5.317294597625732 seconds ---\n"
     ]
    }
   ],
   "source": [
    "# Load data\n",
    "\n",
    "start_time = time.time()\n",
    "\n",
    "def natural_keys(text):\n",
    "    '''\n",
    "    alist.sort(key=natural_keys) sorts in human order\n",
    "    http://nedbatchelder.com/blog/200712/human_sorting.html\n",
    "    (See Toothy's implementation in the comments)\n",
    "    '''\n",
    "    def atoi(text):\n",
    "        return int(text) if text.isdigit() else text\n",
    "    return [atoi(c) for c in re.split(r'(\\d+)', text)]\n",
    "\n",
    "def preprocess(serialized_example):\n",
    "    features = tf.io.parse_example(\n",
    "        serialized_example,\n",
    "        features={\n",
    "            'gram': tf.io.FixedLenFeature([], tf.string),\n",
    "            'label': tf.io.FixedLenFeature([], tf.string)\n",
    "        })\n",
    "    data = tf.io.decode_raw(features['gram'], tf.float32)\n",
    "    label = tf.io.decode_raw(features['label'], tf.float32)\n",
    "    data.set_shape((None, 999))\n",
    "    label.set_shape((None, 4))\n",
    "    return data, label\n",
    "\n",
    "def merge_datasets(vfiles, projy=lambda x: x, asnumpy=False):\n",
    "    yv = []\n",
    "    yl = []\n",
    "    for filename in vfiles:\n",
    "        ds = tf.data.TFRecordDataset(filename)\n",
    "        ds = ds.batch(batch_size=2**13)\n",
    "        ds = ds.map(map_func=preprocess)\n",
    "        for v,l in ds:\n",
    "            v = np.array(v).astype(np.float64)\n",
    "            l = np.array(l).astype(np.float64)\n",
    "            l = np.apply_along_axis(projy, axis=1, arr=l) \n",
    "\n",
    "            # get n unique records per label (tuple), or all data (-1)\n",
    "            if n_records_per_label_per_probe > -1:\n",
    "                # just to make sure there is only once concertraiton tuple in this record\n",
    "                _, ulidx = np.unique(l, return_index=True, axis=0)\n",
    "                u_v = []\n",
    "                u_l = []\n",
    "                for idx in ulidx:\n",
    "                    # TODO: I think this can come from unique\n",
    "                    one_label_idxs = np.where((l == l[idx,:]).all(axis=1))[0]\n",
    "                    u_v.append(v[one_label_idxs[:n_records_per_label_per_probe], :])\n",
    "                    u_l.append(l[one_label_idxs[:n_records_per_label_per_probe], :])\n",
    "                v = np.concatenate(u_v)\n",
    "                l = np.concatenate(u_l)\n",
    "\n",
    "            yv.append(v)\n",
    "            yl.append(l)\n",
    "\n",
    "    x = np.vstack(yv)\n",
    "    y = np.vstack(yl)\n",
    "\n",
    "    if asnumpy:\n",
    "        return x,y\n",
    "    else:\n",
    "        d = tf.data.Dataset.from_tensor_slices((x, y))\n",
    "        return d\n",
    "\n",
    "x_val, y_val = np.array([]), np.array([])\n",
    "if hold_probe < 0:\n",
    "    all_files = sum([\n",
    "        sorted(tf.io.gfile.glob(os.path.join(data_prefix, probe, 'total_records', '*')),\n",
    "               key=natural_keys) for probe in good_probes\n",
    "    ], [])\n",
    "\n",
    "    x, y = merge_datasets(all_files, asnumpy=True)\n",
    "    \n",
    "    idxs = np.random.permutation(x.shape[0])\n",
    "    lim = int(x.shape[0]*(1-val_ratio))\n",
    "    d1idx = idxs[idxs[:lim]]\n",
    "    d2idx = idxs[idxs[lim:]]\n",
    "    x_train, y_train, x_test, y_test = x[d1idx,:], y[d1idx,:], x[d2idx,:], y[d2idx,:]\n",
    "else:\n",
    "    hold_probe_str = good_probes.pop(hold_probe)\n",
    "    train_val_files = sum([\n",
    "        sorted(tf.io.gfile.glob(os.path.join(data_prefix, probe, 'total_records', '*')),\n",
    "               key=natural_keys) for probe in good_probes\n",
    "    ], [])\n",
    "    x, y = merge_datasets(train_val_files, asnumpy=True)\n",
    "\n",
    "    if val_ratio > 0.:\n",
    "        idxs = np.random.permutation(x.shape[0])\n",
    "        lim = int(x.shape[0]*(1-val_ratio))\n",
    "        d1idx = idxs[idxs[:lim]]\n",
    "        d2idx = idxs[idxs[lim:]]\n",
    "        x_train, y_train, x_val, y_val = x[d1idx,:], y[d1idx,:], x[d2idx,:], y[d2idx,:]\n",
    "    else:\n",
    "        x_train, y_train = x, y\n",
    "\n",
    "    hold_files = sum([\n",
    "        sorted(tf.io.gfile.glob(os.path.join(data_prefix, probe, 'total_records', '*')),\n",
    "               key=natural_keys) for probe in [hold_probe_str]\n",
    "    ], [])\n",
    "    x_test, y_test = merge_datasets(hold_files, asnumpy=True)\n",
    "\n",
    "print(x_train.shape, y_train.shape)\n",
    "print(x_test.shape, y_test.shape)\n",
    "if x_val.size > 0:\n",
    "    print(x_val.shape, y_val.shape)\n",
    "\n",
    "print(\"--- %s seconds ---\" % (time.time() - start_time))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Spliting training data into cells per concentration...  done.\n"
     ]
    }
   ],
   "source": [
    "print('Spliting training data into cells per concentration... ', end='')\n",
    "\n",
    "# nint = 5\n",
    "nint = 2\n",
    "analytes = [0,1,3]\n",
    "concentrations = [[]]*len(analytes)\n",
    "ranges = [[]]*len(analytes)\n",
    "for ic in range(len(analytes)):\n",
    "    concentrations[ic] = np.linspace(min(y_train[:,analytes[ic]]), max(y_train[:,analytes[ic]]), nint+1)\n",
    "#     print(c[ic])\n",
    "    ranges[ic] = [concentrations[ic][ir] + np.mean(np.diff(concentrations[ic]))/2 for ir in range(len(concentrations[ic])-1)]\n",
    "#     print(r[ic])\n",
    "\n",
    "# use itertools product instead of nested for loops to make it easier to change the number of analytes\n",
    "from itertools import product\n",
    "split_data = {x: np.array([]) for x in product(range(nint), repeat=len(analytes))}\n",
    "error = []\n",
    "# nsamples = []\n",
    "# print(len(split_data))\n",
    "for x in split_data:\n",
    "\n",
    "    # build binary mask for this cell\n",
    "    good_data = np.full((y_train.shape[0],), True)\n",
    "#     print(good_data.shape, good_data[0].__class__)\n",
    "    for (ix,nx) in enumerate(x):\n",
    "#         print(sum(good_data))\n",
    "#         print(a[ix], nx)\n",
    "        zeros = y_train[:,analytes[ix]] == 0.0\n",
    "        above = y_train[:,analytes[ix]] > concentrations[ix][nx]\n",
    "#         print(above.shape, above[0].__class__, c[ix][nx])\n",
    "        below = y_train[:,analytes[ix]] < concentrations[ix][nx+1]\n",
    "#         print(below.shape, below[0].__class__, c[ix][nx+1])\n",
    "        interval = np.bitwise_or(np.bitwise_and(above, below), zeros)\n",
    "#         print(interval.shape, interval[0].__class__, sum(interval))\n",
    "        good_data = np.bitwise_and(good_data, interval)\n",
    "\n",
    "#     print(sum(good_data))\n",
    "    # save data\n",
    "    split_data[x] = (x_train[good_data,:], y_train[good_data,:])\n",
    "    \n",
    "    y = split_data[x][1]\n",
    "    means = np.array([np.mean(y[y[:,x]>0, x]) for x in range(4)])\n",
    "#     print(means)\n",
    "    expected_means = np.array([ranges[i][x[i]] for i in range(len(analytes))])\n",
    "#     print(expected_means)\n",
    "    error.append(np.sqrt( (means[np.array(analytes)] - expected_means)**2 ))\n",
    "#     nsamples.append(y.shape[0])\n",
    "\n",
    "print(' done.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1127.5, 3382.5], [1123.0, 3369.0], [1129.25, 3387.75]]\n",
      "[array([   0., 2255., 4510.]), array([   0., 2246., 4492.]), array([   0. , 2258.5, 4517. ])]\n"
     ]
    }
   ],
   "source": [
    "print(ranges)\n",
    "print(concentrations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 0, 0) 1103\n",
      "(0, 0, 1) 1163\n",
      "(0, 1, 0) 1160\n",
      "(0, 1, 1) 1066\n",
      "(1, 0, 0) 1163\n",
      "(1, 0, 1) 1067\n",
      "(1, 1, 0) 1070\n",
      "(1, 1, 1) 972\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2d0a842d696345f7901834ea6a94944d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# inspect split data\n",
    "nsamples = np.zeros(len(split_data))\n",
    "for (ix, x) in enumerate(split_data):\n",
    "    print(x, split_data[x][1].shape[0])\n",
    "    \n",
    "e = np.array(error)\n",
    "fig, axs = plt.subplots(1, 3, figsize=(20,5))\n",
    "for idx in range(3):\n",
    "    axs[idx].hist(e[:,idx])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.00000000e-04 4.64158883e-04 2.15443469e-03 1.00000000e-02\n",
      " 4.64158883e-02 2.15443469e-01 1.00000000e+00 4.64158883e+00\n",
      " 2.15443469e+01 1.00000000e+02]\n"
     ]
    }
   ],
   "source": [
    "lambdas = np.logspace(-4, 2, 10)\n",
    "# lambdas = []\n",
    "nlambdas = len(lambdas)\n",
    "if nlambdas:\n",
    "    print(lambdas)\n",
    "    min_lambda = min(lambdas)\n",
    "    max_lambda = max(lambdas)\n",
    "else:\n",
    "    min_lambda = 0\n",
    "    max_lambda = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/mnt/nfs/proj/in-vitro/Leonardo/glmnet/fits/cf_alpha_1.0_CF025_1_0.0001_100.0_10_2x2x2.pickle\n",
      "Computing 10 cross-validations for model 0 of 8 ((0, 0, 0)) ... fitting model with 10 pre-defined lambdas\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[status]\tParallel glmnet cv with 28 cores\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " took 390.96596217155457 seconds\n",
      "Computing 10 cross-validations for model 1 of 8 ((0, 0, 1)) ... fitting model with 10 pre-defined lambdas\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[status]\tParallel glmnet cv with 28 cores\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " took 437.1633641719818 seconds\n",
      "Computing 10 cross-validations for model 2 of 8 ((0, 1, 0)) ... fitting model with 10 pre-defined lambdas\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[status]\tParallel glmnet cv with 28 cores\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " took 217.40631675720215 seconds\n",
      "Computing 10 cross-validations for model 3 of 8 ((0, 1, 1)) ... fitting model with 10 pre-defined lambdas\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[status]\tParallel glmnet cv with 28 cores\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " took 233.76714086532593 seconds\n",
      "Computing 10 cross-validations for model 4 of 8 ((1, 0, 0)) ... fitting model with 10 pre-defined lambdas\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[status]\tParallel glmnet cv with 28 cores\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " took 263.94223618507385 seconds\n",
      "Computing 10 cross-validations for model 5 of 8 ((1, 0, 1)) ... fitting model with 10 pre-defined lambdas\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[status]\tParallel glmnet cv with 28 cores\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " took 370.0534858703613 seconds\n",
      "Computing 10 cross-validations for model 6 of 8 ((1, 1, 0)) ... fitting model with 10 pre-defined lambdas\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[status]\tParallel glmnet cv with 28 cores\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " took 175.5290012359619 seconds\n",
      "Computing 10 cross-validations for model 7 of 8 ((1, 1, 1)) ... fitting model with 10 pre-defined lambdas\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[status]\tParallel glmnet cv with 28 cores\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " took 180.25086522102356 seconds\n",
      "All cvs took 2269.0804686546326 seconds\n"
     ]
    }
   ],
   "source": [
    "# fit GLMNET in parallel (ncores) with cross validation to find lambda\n",
    "\n",
    "output_file = os.path.join(results_prefix, f'cf_alpha_1.0_{hold_probe_str}_{n_records_per_label_per_probe}_{min_lambda}_{max_lambda}_{nlambdas}_{nint}x{nint}x{nint}.pickle')\n",
    "print(output_file)\n",
    "\n",
    "# alphas = [.9, 1.]\n",
    "alpha = 1.0\n",
    "models = {x: [] for x in split_data}\n",
    "ncross = 10\n",
    "nx = len(split_data)\n",
    "all_start_time = time.time()\n",
    "for (ix, x) in enumerate(split_data):\n",
    "    this_x_train, this_y_train = split_data[x]\n",
    "    print(f'Computing {ncross} cross-validations for model {ix} of {nx} ({x}) ... ', end='')\n",
    "    start_time = time.time()\n",
    "    if nlambdas:\n",
    "        print(f'fitting model with {nlambdas} pre-defined lambdas')\n",
    "        models[x] = cvglmnet(x = this_x_train.copy(), y = this_y_train.copy(), family='mgaussian', parallel=ncores, ptype = 'mse', nfolds = ncross, alpha=alpha, lambdau=lambdas)\n",
    "    else:\n",
    "        print('fitting model with 100 lambdas decided by the first run with all the data')\n",
    "        models[x] = cvglmnet(x = this_x_train.copy(), y = this_y_train.copy(), family='mgaussian', parallel=ncores, ptype = 'mse', nfolds = ncross, alpha=alpha)        \n",
    "    with open(output_file, 'wb') as f:\n",
    "        pickle.dump((models, ranges, analytes), f)\n",
    "    print(\" took %s seconds\" % (time.time() - start_time))\n",
    "print(\"All cvs took %s seconds\" % (time.time() - all_start_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DEBUG\n",
    "# _, uidx = np.unique(y_train, return_index=True, axis=0)\n",
    "# for (idx, x) in enumerate(y_train[uidx, :]):\n",
    "#     if x[0] < 2000 and x[3] < 2000 and x[1] > 3000:\n",
    "#         print(idx, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating predictions for each one of the 8 models\n",
      "0 of 8 ((0, 0, 0))...  took 0.09232568740844727 seconds\n",
      "1 of 8 ((0, 0, 1))...  took 0.07277297973632812 seconds\n",
      "2 of 8 ((0, 1, 0))...  took 0.07445716857910156 seconds\n",
      "3 of 8 ((0, 1, 1))...  took 0.07344317436218262 seconds\n",
      "4 of 8 ((1, 0, 0))...  took 0.07306981086730957 seconds\n",
      "5 of 8 ((1, 0, 1))...  took 0.08481907844543457 seconds\n",
      "6 of 8 ((1, 1, 0))...  took 0.07306742668151855 seconds\n",
      "7 of 8 ((1, 1, 1))...  took 0.08271408081054688 seconds\n"
     ]
    }
   ],
   "source": [
    "# Compute the predictions of each model in each one of the new/hold out records\n",
    "\n",
    "x_to_pred = x_test\n",
    "y_to_pred = y_test\n",
    "\n",
    "# x_to_pred = x_val\n",
    "# y_to_pred = y_val\n",
    "\n",
    "# # DEBUG\n",
    "# x_to_pred = x_train\n",
    "# y_to_pred = y_train\n",
    "\n",
    "# x_to_pred = x_train[uidx[307], :][None, :]\n",
    "# y_to_pred = y_train[uidx[307], :][None, :]\n",
    "# x_to_pred = x_train[[uidx[307], uidx[210]], :]\n",
    "# y_to_pred = y_train[[uidx[307], uidx[210]], :]\n",
    "\n",
    "# generate the predictions for each model\n",
    "y_hats = np.zeros((len(models), y_to_pred.shape[0], y_to_pred.shape[1]))\n",
    "nx = len(models)\n",
    "print(f'Generating predictions for each one of the {nx} models')\n",
    "for (ix, x) in enumerate(models):\n",
    "    print(f'{ix} of {nx} ({x})... ', end='')\n",
    "    model = models[x]\n",
    "    start_time = time.time()\n",
    "#     y_hat = cvglmnetPredict(fit, newx = x_to_pred, s='lambda_min')\n",
    "    y_hat = cvglmnetPredict(model, newx = x_to_pred, s='lambda_1se')\n",
    "    \n",
    "    # DEBUG\n",
    "#     y_hat = y_test[:,:,None]*(1.+0.1*np.random.randn(y_test.shape[0], y_test.shape[1], 1))\n",
    "\n",
    "    y_hats[ix,:,:] = y_hat[:,:,0]\n",
    "    print(\" took %s seconds\" % (time.time() - start_time))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing differences between predictions of the models and the center of the range in which the models were trained... done\n"
     ]
    }
   ],
   "source": [
    "print('Computing differences between predictions of the models and the center of the range in which the models were trained... ', end='')\n",
    "# compute the differences to original intervals\n",
    "diff_y_hats = np.zeros(y_hats[:,:,analytes].shape)\n",
    "diff_y_true = np.zeros(y_hats[:,:,analytes].shape)\n",
    "for (ix,x) in enumerate(models):\n",
    "#     print(x)\n",
    "    for ia in range(len(analytes)):\n",
    "#         print(ix, ia, analytes[ia], x[ia], ranges[ia][x[ia]])\n",
    "        diff_y_hats[ix,:,ia] = (y_hats[ix,:,analytes[ia]] - ranges[ia][x[ia]])**2\n",
    "        # save that for comparing how often the predicted model matches the \"ideal\" model \n",
    "        diff_y_true[ix,:,ia] = (y_to_pred[:,analytes[ia]] - ranges[ia][x[ia]])**2 \n",
    "print('done')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DA: 893.91162\n",
      "5HT: 518.89891\n",
      "pH: 0.18095\n",
      "NE: 1092.24035\n"
     ]
    }
   ],
   "source": [
    "# Choose N models: For each new/hold out records, chose the model which has the prediction closest to the range in which the model was trained \n",
    "\n",
    "# find the model with sum of predictions closest to original intervals\n",
    "model_e = np.sqrt(np.sum(diff_y_hats,axis=2)) # take sqrt to make it easier to debug, doesn't change the ordering\n",
    "min_idx = model_e.argmin(axis=0)\n",
    "\n",
    "# compute the RMSE for each sample (there is certainly a more pythonic way to do that, but that works)\n",
    "rmse = np.zeros((4,), dtype=np.float64)\n",
    "for (sample_idx, model_idx) in enumerate(min_idx):\n",
    "    rmse += (y_hats[model_idx, sample_idx, :] - y_to_pred[sample_idx, :])**2/y_to_pred.shape[0]\n",
    "\n",
    "rmse = np.sqrt(rmse) \n",
    "\n",
    "for (armse, name) in zip(rmse,names):\n",
    "    print('%s: %4.5f'%(name,armse), flush=True)\n",
    "\n",
    "# DA: 776.18963\n",
    "# 5HT: 342.67418\n",
    "# pH: 0.07271\n",
    "# NE: 920.62202\n",
    "\n",
    "# DA: 791.33012\n",
    "# 5HT: 418.37242\n",
    "# pH: 0.10728\n",
    "# NE: 1066.39777"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "How often the predicted model matches the ideal model? 0.6685615848406546\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fd740ada6e62439daa170e04ee9ff0f1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# find the model with sum of predictions closest to original intervals\n",
    "real_model_e = np.sqrt(np.sum(diff_y_true,axis=2)) # take sqrt to make it easier to debug, doesn't change the ordering\n",
    "real_min_idx = real_model_e.argmin(axis=0)\n",
    "ratio_same = sum(real_min_idx == min_idx)/min_idx.shape[0]\n",
    "print(f'How often the predicted model matches the expected model? {ratio_same}')\n",
    "# fig, axs = plt.subplots(1, 2, sharey=False, sharex=False, tight_layout=True, figsize=(20,10))\n",
    "x = np.array(list(split_data.keys()))\n",
    "d = np.zeros(min_idx.shape)\n",
    "for sample in range(len(d)):\n",
    "    d[sample] = sum(abs(x[min_idx[sample], :] - x[real_min_idx[sample], :]))\n",
    "fig, axs = plt.subplots(figsize=(10,5))\n",
    "axs.hist(d)\n",
    "axs.set_xticks(range(0, x.shape[1]+1))\n",
    "axs.set_title('L1 distance between predicted model and expected model')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DEBUG\n",
    "# print(y_to_pred)\n",
    "# print(y_hats)\n",
    "# x = np.sqrt(diff_y_hats)\n",
    "# print(x)\n",
    "# print(model_e)\n",
    "# print(min_idx)\n",
    "# # compute the RMSE for each sample (there is certainly a more pythonic way to do that, but that works)\n",
    "# rmse = np.zeros((4,), dtype=np.float64)\n",
    "# for (sample_idx, model_idx) in enumerate(min_idx):\n",
    "#     print(sample_idx, model_idx)\n",
    "#     rmse += (y_hats[model_idx, sample_idx, :] - y_to_pred[sample_idx, :])**2/y_to_pred.shape[0]\n",
    "\n",
    "# rmse = np.sqrt(rmse) \n",
    "# print(rmse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DEBUG\n",
    "# x1=(np.array([1969.4478911  , 905.53209513 ,  7.41054667 ,3656.49362621]) - np.array([2232.    ,     729.     ,      7.4000001, 3727.  ]))**2\n",
    "# x2=(np.array([ -55.492296  ,  4099.18583448   ,  7.40349391 ,   80.2708897]) - np.array([    0.    ,    4099.       ,    7.4000001  ,  0.           ]))**2\n",
    "# print(np.sqrt((x1+x2)/y_to_pred.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose only ONE model: The one with the minimum TOTAL error across ALL new/hold out records\n",
    "\n",
    "# find the model with sum of predictions closest to original intervals\n",
    "model_e = np.sqrt(np.sum(np.sum(diff_y_hats,axis=2), axis=1)) # take sqrt to make it easier to debug, doesn't change the ordering\n",
    "print(model_e)\n",
    "min_idx = model_e.argmin(axis=0)\n",
    "print(min_idx)\n",
    "\n",
    "# compute the RMSE for each sample (there is certainly a more pythonic way to do that, but that works)\n",
    "rmse = np.sqrt(np.mean((np.squeeze(y_hats[min_idx,:,:])-y_to_pred)**2,axis=0))\n",
    "for (armse, name) in zip(rmse,names):\n",
    "    print('%s: %4.5f'%(name,armse), flush=True)\n",
    "\n",
    "# DA: 833.85557\n",
    "# 5HT: 499.82708\n",
    "# pH: 0.11135\n",
    "# NE: 1047.89369\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# inspect the weights for a model in given interval of concentrations (cell in the grid) \n",
    "fig, axs = plt.subplots(2, 2, sharey=False, sharex=False, tight_layout=True, figsize=(20,10))\n",
    "for cell in models:\n",
    "    l = np.where([x == models[cell]['lambda_1se'][0] for x in models[cell]['lambdau']])[0][0]\n",
    "    for ia in range(4):\n",
    "        idx = np.unravel_index(ia, axs.shape)\n",
    "        x = models[cell]['glmnet_fit']['beta'][ia][:, l]\n",
    "        axs[idx].plot(x)\n",
    "        axs[idx].set_title(\"Coeficients %s\"%names[ia])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "784e482424ef484aa86b0990c05d325a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# plt.close('all')\n",
    "clrs = sns.color_palette(\"husl\", len(models))\n",
    "fig, axs = plt.subplots(1, 2, tight_layout=True, figsize=(15,5))\n",
    "for (iproj, proj) in enumerate([lambda x: x, np.log]):\n",
    "    for (icell, cell) in enumerate(models):\n",
    "        e = proj(models[cell]['cvm'])\n",
    "        estd = proj(models[cell]['cvm'])\n",
    "        ls = proj(models[cell]['lambdau'])\n",
    "        l1 = proj(models[cell]['lambda_1se'][0])\n",
    "        l1idx = np.where([x == l1 for x in ls])[0][0]\n",
    "        lm = proj(models[cell]['lambda_min'][0])\n",
    "        lmidx = np.where([x == lm for x in ls])[0][0]\n",
    "\n",
    "        with sns.axes_style(\"darkgrid\"):\n",
    "#             epochs = list(range(len(e)))\n",
    "            axs[iproj].plot(ls, e, label=str(cell)+' '+('%1.3f'%l1), c=clrs[icell])\n",
    "#             if icell == 1:\n",
    "#                 axs[iproj].fill_between(epochs, e-estd, e+estd ,alpha=0.3, facecolor=clrs[icell])\n",
    "            axs[iproj].plot(ls[l1idx], e[l1idx], 'k*')\n",
    "            axs[iproj].plot(ls[lmidx], e[lmidx], 'r*')\n",
    "    axs[iproj].legend()\n",
    "    axs[iproj].set_title(proj.__name__)\n",
    "    axs[iproj].set_ylabel('Average error across cross-validation')\n",
    "    axs[iproj].set_xlabel('Different lambdas')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}