mimic3_dataset.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Re-implmentation of AI Clinician Matlab Code in Python \n",
    "# Author: KyungJoong Kim (GIST, South Korea)\n",
    "# Date: 2020 June 2 \n",
    "# \n",
    "# This code is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;\n",
    "# without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE\n",
    "\n",
    "# Note \n",
    "# \n",
    "# python interp1d() and matlab interp1() are very close each other but not the same \n",
    "# => possible difference at 15 digits or more after the decimal point \n",
    "# => usually, it doesn't make any change on the outcomes but has small impact on sorting function and equal operator (=) \n",
    "#\n",
    "\n",
    "\n",
    "# AI Clinician Building MIMIC-III dataset \n",
    "\n",
    "# (c) Matthieu Komorowski, Imperial College London 2015-2019\n",
    "# as seen in publication: https://www.nature.com/articles/s41591-018-0213-5\n",
    "\n",
    "# version 16 Feb 19\n",
    "# uses the sepsis-3 cohort previously defined\n",
    "# builds the MIMIC-III dataset\n",
    "\n",
    "# This code is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;\n",
    "# without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE\n",
    "\n",
    "import pickle \n",
    "import numpy as np \n",
    "import pandas as pd \n",
    "import scipy.io as sio \n",
    "import matplotlib \n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats.mstats import mquantiles\n",
    "from scipy.interpolate import interp1d \n",
    "from scipy.spatial.distance import pdist, squareform\n",
    "\n",
    "# To ignore 'Runtime Warning: Invalid value encountered in greater' caused by NaN \n",
    "np.warnings.filterwarnings('ignore')\n",
    "\n",
    "# Load pickle for step 3 \n",
    "with open('step_3_start.pkl', 'rb') as file:\n",
    "    sample_and_hold = pickle.load(file)\n",
    "    demog= pickle.load(file)\n",
    "    ce010= pickle.load(file)\n",
    "    ce1020= pickle.load(file)\n",
    "    ce2030= pickle.load(file)\n",
    "    ce3040= pickle.load(file)\n",
    "    ce4050= pickle.load(file)\n",
    "    ce5060= pickle.load(file)\n",
    "    ce6070= pickle.load(file)\n",
    "    ce7080= pickle.load(file)\n",
    "    ce8090= pickle.load(file)\n",
    "    ce90100= pickle.load(file)\n",
    "    labU= pickle.load(file)\n",
    "    MV= pickle.load(file)\n",
    "    inputpreadm= pickle.load(file)\n",
    "    inputMV= pickle.load(file)\n",
    "    inputCV= pickle.load(file)\n",
    "    vasoMV= pickle.load(file)\n",
    "    vasoCV= pickle.load(file)\n",
    "    UOpreadm= pickle.load(file)\n",
    "    UO= pickle.load(file)\n",
    "    sepsis= pickle.load(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ########################################################################\n",
    "#           INITIAL REFORMAT WITH CHARTEVENTS, LABS AND MECHVENT\n",
    "# ########################################################################\n",
    "\n",
    "# gives an array with all unique charttime (1 per row) and all items in columns.\n",
    "# ################## IMPORTANT !!!!!!!!!!!!!!!!!!\n",
    "# Here i use -24 -> +48 because that's for the MDP\n",
    "\n",
    "\n",
    "reformat = np.full((2000000,68),np.nan)\n",
    "qstime=np.zeros((100000,4)) \n",
    "winb4=25;   #lower limit for inclusion of data (48h before time flag)\n",
    "winaft=49;  # upper limit (24h after)\n",
    "irow=0;  #recording row for summary table\n",
    "\n",
    "for icustayidrow in range(1,sepsis.shape[0]+1):\n",
    "    \n",
    "    if(icustayidrow%1000==0):\n",
    "        print(icustayidrow)\n",
    "    \n",
    "    qst=sepsis.loc[icustayidrow-1,'sepsis_time'] #flag for presumed infection\n",
    "    icustayid = int(sepsis.loc[icustayidrow-1,'icustayid'])-200000 \n",
    "    \n",
    "\n",
    "    if icustayid<10000:\n",
    "        temp = ce010[ce010[:,0]==icustayid+200000,:] \n",
    "    elif icustayid>=10000 and icustayid<20000: \n",
    "        temp = ce1020[ce1020[:,0]==icustayid+200000,:]\n",
    "    elif icustayid>=20000 and icustayid<30000: \n",
    "        temp = ce2030[ce2030[:,0]==icustayid+200000,:]\n",
    "    elif icustayid>=30000 and icustayid<40000: \n",
    "        temp = ce3040[ce3040[:,0]==icustayid+200000,:] \n",
    "    elif icustayid>=40000 and icustayid<50000: \n",
    "        temp = ce4050[ce4050[:,0]==icustayid+200000,:]             \n",
    "    elif icustayid>=50000 and icustayid<60000: \n",
    "        temp = ce5060[ce5060[:,0]==icustayid+200000,:] \n",
    "    elif icustayid>=60000 and icustayid<70000: \n",
    "        temp = ce6070[ce6070[:,0]==icustayid+200000,:]           \n",
    "    elif icustayid>=70000 and icustayid<80000: \n",
    "        temp = ce7080[ce7080[:,0]==icustayid+200000,:]             \n",
    "    elif icustayid>=80000 and icustayid<90000: \n",
    "        temp = ce8090[ce8090[:,0]==icustayid+200000,:]             \n",
    "    elif icustayid>=90000: \n",
    "        temp = ce90100[ce90100[:,0]==icustayid+200000,:] \n",
    "                \n",
    "    ii = (temp[:,1]>=qst-(winb4+4)*3600) & (temp[:,1]<=qst+(winaft+4)*3600) # time period of interest -4h and +4h\n",
    "            \n",
    "    temp = temp[ii,:] # only time period of interest \n",
    "    \n",
    "    # LABEVENTS\n",
    "    ii=labU[:,0]==icustayid+200000\n",
    "    temp2=labU[ii,:]\n",
    "    ii=(temp2[:,1]>=qst-(winb4+4)*3600) & (temp2[:,1]<=qst+(winaft+4)*3600) # time period of interest -4h and +4h\n",
    "    temp2=temp2[ii,:]   # only time period of interest\n",
    "    \n",
    "    #Mech Vent + ?extubated\n",
    "    ii=MV[:,0]==icustayid+200000\n",
    "    temp3=MV[ii,:]\n",
    "    ii=(temp3[:,1]>=qst-(winb4+4)*3600) & (temp3[:,1]<=qst+(winaft+4)*3600) # time period of interest -4h and +4h\n",
    "            \n",
    "    temp3=temp3[ii,:]   # only time period of interest \n",
    "                \n",
    "    temp_list = [] \n",
    "    if(temp.size!=0) :\n",
    "        temp_list.append(temp[:,1].reshape(-1,1)) \n",
    "    if(temp2.size!=0) : \n",
    "        temp_list.append(temp2[:,1].reshape(-1,1)) \n",
    "    if(temp3.size!=0) :\n",
    "        temp_list.append(temp3[:,1].reshape(-1,1)) \n",
    "\n",
    "    if(len(temp_list)==0):\n",
    "        t= np.array([]) \n",
    "    else : \n",
    "        t = np.unique(np.vstack(temp_list))  # list of unique timestamps from all 3 sources / sorted in ascending order\n",
    "                \n",
    "\n",
    "    if(t.size!=0):\n",
    "        for i in range(t.size):\n",
    "            # CHARTEVENTS\n",
    "            ii = temp[:,1]==t[i] \n",
    "            col=temp[ii,2].astype('int64')\n",
    "            value=temp[ii,3] \n",
    "            reformat[irow,0]=i+1; # timestep\n",
    "            reformat[irow,1]=icustayid\n",
    "            reformat[irow,2]=t[i]  #charttime\n",
    "            for index,c in enumerate(col): \n",
    "                reformat[irow,2+c]=value[index] # (locb(:,1)); %store available values\n",
    "                            \n",
    "            # LAB VALUES\n",
    "            ii = temp2[:,1]==t[i] \n",
    "            col=temp2[ii,2].astype('int64') \n",
    "            value=temp2[ii,3]\n",
    "            for index,c in enumerate(col): \n",
    "                reformat[irow,30+c]=value[index] # store available values\n",
    "                                \n",
    "            # MV  \n",
    "            ii = temp3[:,1]==t[i] \n",
    "            if np.nansum(ii)>0:\n",
    "                value=temp3[ii,2:4]\n",
    "                reformat[irow,66:68]=value # store available values\n",
    "            else:\n",
    "                reformat[irow,66:68]=np.nan \n",
    "                                \n",
    "            irow=irow+1\n",
    "                \n",
    "                     \n",
    "        qstime[icustayid-1,0]=qst #time of sepsis\n",
    "        # HERE I SAVE FIRST and LAST TIMESTAMPS, in QSTIME, for each ICUSTAYID\n",
    "        qstime[icustayid-1,1]=t[0]  #first timestamp\n",
    "        qstime[icustayid-1,2]=t[-1] # last timestamp\n",
    "        qstime[icustayid-1,3]=demog.iloc[np.where(demog['icustay_id']==icustayid+200000)[0],4] # dischargetime \n",
    "\n",
    "                \n",
    "reformat = reformat[:irow,:]  #delete extra unused rows\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ########################################################################\n",
    "#                                   OUTLIERS \n",
    "# ########################################################################\n",
    "# \n",
    "\n",
    "def deloutabove(reformat, var, thres):\n",
    "    # DELOUTABOVE delete values above the given threshold, for column 'var'\n",
    "    ii = reformat[:,var]>thres\n",
    "    reformat[ii, var] = np.nan \n",
    "    return reformat\n",
    "\n",
    "def deloutbelow(reformat, var, thres):\n",
    "    # DELOUTBELOW delete values below the given threshold, for column 'var'\n",
    "\n",
    "    ii = reformat[:,var]<thres\n",
    "    reformat[ii, var] = np.nan \n",
    "    return reformat\n",
    "\n",
    "\n",
    "# weight\n",
    "reformat = deloutabove(reformat,4,300) # delete outlier above a threshold (300 kg), for variable # 5 \n",
    "# HR\n",
    "reformat=deloutabove(reformat,7,250)\n",
    "\n",
    "#BP\n",
    "reformat=deloutabove(reformat,8,300)\n",
    "reformat=deloutbelow(reformat,9,0)\n",
    "reformat=deloutabove(reformat,9,200)\n",
    "reformat=deloutbelow(reformat,10,0)\n",
    "reformat=deloutabove(reformat,10,200)\n",
    "\n",
    "#RR\n",
    "reformat=deloutabove(reformat,11,80)\n",
    "\n",
    "#SpO2\n",
    "reformat=deloutabove(reformat,12,150)\n",
    "ii=reformat[:,12]>100\n",
    "reformat[ii,12]=100\n",
    "reformat=deloutbelow(reformat,12,50) # added in step 3 \n",
    "\n",
    "#temp\n",
    "ii=(reformat[:,13]>90) & (np.isnan(reformat[:,14]))\n",
    "reformat[ii,14]=reformat[ii,13]\n",
    "reformat=deloutabove(reformat,13,90)\n",
    "reformat=deloutbelow(reformat,13,25) # added in step 3 \n",
    "\n",
    "#interface / is in col 22\n",
    "\n",
    "# FiO2\n",
    "reformat=deloutabove(reformat,22,100)\n",
    "ii=reformat[:,22]<1\n",
    "reformat[ii,22]=reformat[ii,22]*100\n",
    "reformat=deloutbelow(reformat,22,20)\n",
    "reformat=deloutabove(reformat,23,1.5)\n",
    "\n",
    "# O2 FLOW\n",
    "reformat=deloutabove(reformat,24,70)\n",
    "\n",
    "# PEEP\n",
    "reformat=deloutbelow(reformat,25,0)\n",
    "reformat=deloutabove(reformat,25,40)\n",
    "\n",
    "#TV\n",
    "reformat=deloutabove(reformat,26,1800)\n",
    "\n",
    "#MV\n",
    "reformat=deloutabove(reformat,27,50)\n",
    "\n",
    "#K+\n",
    "reformat=deloutbelow(reformat,31,1)\n",
    "reformat=deloutabove(reformat,31,15)\n",
    "\n",
    "#Na\n",
    "reformat=deloutbelow(reformat,32,95)\n",
    "reformat=deloutabove(reformat,32,178)\n",
    "\n",
    "#Cl\n",
    "reformat=deloutbelow(reformat,33,70)\n",
    "reformat=deloutabove(reformat,33,150)\n",
    "\n",
    "#Glc\n",
    "reformat=deloutbelow(reformat,34,1)\n",
    "reformat=deloutabove(reformat,34,1000)\n",
    "\n",
    "#Creat\n",
    "reformat=deloutabove(reformat,36,150)\n",
    "\n",
    "#Mg\n",
    "reformat=deloutabove(reformat,37,10)\n",
    "\n",
    "#Ca\n",
    "reformat=deloutabove(reformat,38,20)\n",
    "\n",
    "#ionized Ca\n",
    "reformat=deloutabove(reformat,39,5)\n",
    "\n",
    "#CO2\n",
    "reformat=deloutabove(reformat,40,120)\n",
    "\n",
    "#SGPT/SGOT\n",
    "reformat=deloutabove(reformat,41,10000)\n",
    "reformat=deloutabove(reformat,42,10000)\n",
    "\n",
    "#Hb/Ht\n",
    "reformat=deloutabove(reformat,49,20)\n",
    "reformat=deloutabove(reformat,50,65)\n",
    "\n",
    "#WBC\n",
    "reformat=deloutabove(reformat,52,500)\n",
    "\n",
    "#plt\n",
    "reformat=deloutabove(reformat,53,2000)\n",
    "\n",
    "#INR\n",
    "reformat=deloutabove(reformat,57,20)\n",
    "\n",
    "#pH\n",
    "reformat=deloutbelow(reformat,58,6.7)\n",
    "reformat=deloutabove(reformat,58,8)\n",
    "\n",
    "#po2\n",
    "reformat=deloutabove(reformat,59,700)\n",
    "\n",
    "#pco2\n",
    "reformat=deloutabove(reformat,60,200)\n",
    "\n",
    "#BE\n",
    "reformat=deloutbelow(reformat,61,-50)\n",
    "\n",
    "#lactate\n",
    "reformat=deloutabove(reformat,62,30)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#####################################################################\n",
    "# some more data manip / imputation from existing values\n",
    "\n",
    "# estimate GCS from RASS - data from Wesley JAMA 2003\n",
    "ii=(np.isnan(reformat[:,5])) & (reformat[:,6]>=0)\n",
    "reformat[ii,5]=15\n",
    "ii=(np.isnan(reformat[:,5])) & (reformat[:,6]==-1)\n",
    "reformat[ii,5]=14\n",
    "ii=(np.isnan(reformat[:,5])) & (reformat[:,6]==-2)\n",
    "reformat[ii,5]=12\n",
    "ii=(np.isnan(reformat[:,5])) & (reformat[:,6]==-3)\n",
    "reformat[ii,5]=11\n",
    "ii=(np.isnan(reformat[:,5])) & (reformat[:,6]==-4)\n",
    "reformat[ii,5]=6\n",
    "ii=(np.isnan(reformat[:,5])) & (reformat[:,6]==-5)\n",
    "reformat[ii,5]=3\n",
    "\n",
    "\n",
    "# FiO2\n",
    "ii=(~np.isnan(reformat[:,22])) & (np.isnan(reformat[:,23]))\n",
    "reformat[ii,23]=reformat[ii,22]/100\n",
    "ii=(~np.isnan(reformat[:,23])) & (np.isnan(reformat[:,22]))\n",
    "reformat[ii,22]=reformat[ii,23]*100\n",
    "\n",
    "\n",
    "# Matthieu Komorowski - Imperial College London 2017 \n",
    "# will copy a value in the rows below if the missing values are within the\n",
    "# hold period for this variable (e.g. 48h for weight, 2h for HR...)\n",
    "# vitalslab_hold = 2x55 cell (with row1 = strings of names ; row 2 = hold time)\n",
    "\n",
    "def SAH(reformat, vitalslab_hold): \n",
    "    temp = reformat.copy() \n",
    "\n",
    "    hold=vitalslab_hold[1,:]\n",
    "    nrow=temp.shape[0]\n",
    "    ncol=temp.shape[1]\n",
    "\n",
    "    lastcharttime=np.zeros(ncol) \n",
    "    lastvalue=np.zeros(ncol)\n",
    "    oldstayid=temp[0,1]\n",
    "\n",
    "    for i in range(3,ncol): \n",
    "        if(i%10==0):\n",
    "            print(i)\n",
    "        for j in range(0,nrow): \n",
    "            \n",
    "            if oldstayid!=temp[j,1]: \n",
    "                lastcharttime=np.zeros(ncol) \n",
    "                lastvalue=np.zeros(ncol) \n",
    "                oldstayid=temp[j,1] \n",
    "                              \n",
    "            if np.isnan(temp[j,i])==0: \n",
    "                lastcharttime[i]=temp[j,2] \n",
    "                lastvalue[i]=temp[j,i]\n",
    "        \n",
    "            if j>0:\n",
    "                if (np.isnan(temp[j,i])) and  (temp[j,1]==oldstayid) and ((temp[j,2]-lastcharttime[i])<=hold[i-3]*3600): #note : hold has 53 cols, temp has 55\n",
    "                    temp[j,i]=lastvalue[i]   \n",
    "            \n",
    "    return temp \n",
    "\n",
    "\n",
    "# reformatsah=SAH(reformat,sample_and_hold)  # do SAH first to handle this task\n",
    "reformatsah=SAH(reformat,sample_and_hold)\n",
    "\n",
    "# NO FiO2, YES O2 flow, no interface OR cannula\n",
    "ii=np.where(np.isnan(reformatsah[:,22]) & (~np.isnan(reformatsah[:,24])) & ((reformatsah[:,21]==0) | (reformatsah[:,21]==2)))[0] \n",
    "reformat[ii[reformatsah[ii,24]<=15],22]=70\n",
    "reformat[ii[reformatsah[ii,24]<=12],22]=62\n",
    "reformat[ii[reformatsah[ii,24]<=10],22]=55\n",
    "reformat[ii[reformatsah[ii,24]<=8],22]=50\n",
    "reformat[ii[reformatsah[ii,24]<=6],22]=44\n",
    "reformat[ii[reformatsah[ii,24]<=5],22]=40\n",
    "reformat[ii[reformatsah[ii,24]<=4],22]=36\n",
    "reformat[ii[reformatsah[ii,24]<=3],22]=32\n",
    "reformat[ii[reformatsah[ii,24]<=2],22]=28\n",
    "reformat[ii[reformatsah[ii,24]<=1],22]=24\n",
    "\n",
    "# NO FiO2, NO O2 flow, no interface OR cannula\n",
    "ii=np.where((np.isnan(reformatsah[:,22])) & (np.isnan(reformatsah[:,24])) & ((reformatsah[:,21]==0) | (reformatsah[:,21]==2)))[0] #no fio2 given and o2flow given, no interface OR cannula\n",
    "reformat[ii,22]=21\n",
    "\n",
    "\n",
    "#NO FiO2, YES O2 flow, face mask OR.... OR ventilator (assume it's face mask)\n",
    "ii=np.where((np.isnan(reformatsah[:,22])) & (~np.isnan(reformatsah[:,24])) & ((reformatsah[:,21]==np.nan) | (reformatsah[:,21]==1) | (reformatsah[:,21]==3) | (reformatsah[:,21]==4) | (reformatsah[:,21]==5) | (reformatsah[:,21]==6) | (reformatsah[:,21]==9) | (reformatsah[:,21]==10)))[0] \n",
    "reformat[ii[reformatsah[ii,24]<=15],22]=75\n",
    "reformat[ii[reformatsah[ii,24]<=12],22]=69\n",
    "reformat[ii[reformatsah[ii,24]<=10],22]=66\n",
    "reformat[ii[reformatsah[ii,24]<=8],22]=58\n",
    "reformat[ii[reformatsah[ii,24]<=6],22]=40\n",
    "reformat[ii[reformatsah[ii,24]<=4],22]=36\n",
    "\n",
    "# NO FiO2, NO O2 flow, face mask OR ....OR ventilator\n",
    "ii=np.where((np.isnan(reformatsah[:,22])) & (np.isnan(reformatsah[:,24])) & ((reformatsah[:,21]==np.nan) | (reformatsah[:,21]==1) | (reformatsah[:,21]==3) | (reformatsah[:,21]==4) | (reformatsah[:,21]==5) | (reformatsah[:,21]==6) | (reformatsah[:,21]==9) | (reformatsah[:,21]==10)))[0]  #no fio2 given and o2flow given, no interface OR cannula\n",
    "reformat[ii,22]=np.nan\n",
    "\n",
    "\n",
    "# NO FiO2, YES O2 flow, Non rebreather mask\n",
    "ii=np.where((np.isnan(reformatsah[:,22])) & (~np.isnan(reformatsah[:,24]) & (reformatsah[:,21]==7)))[0]\n",
    "reformat[ii[reformatsah[ii,24]>=10],22]=90\n",
    "reformat[ii[reformatsah[ii,24]>=15],22]=100\n",
    "reformat[ii[reformatsah[ii,24]<10],22]=80\n",
    "reformat[ii[reformatsah[ii,24]<=8],22]=70\n",
    "reformat[ii[reformatsah[ii,24]<=6],22]=60\n",
    "\n",
    "\n",
    "# NO FiO2, NO O2 flow, NRM\n",
    "ii=np.where((np.isnan(reformatsah[:,22])) & (np.isnan(reformatsah[:,24]) & (reformatsah[:,21]==7)))[0]  #no fio2 given and o2flow given, no interface OR cannula\n",
    "reformat[ii,22]=np.nan \n",
    "\n",
    "# update again FiO2 columns\n",
    "ii=(~np.isnan(reformat[:,22])) & (np.isnan(reformat[:,23]))\n",
    "reformat[ii,23]=reformat[ii,22]/100\n",
    "ii=(~np.isnan(reformat[:,23])) & (np.isnan(reformat[:,22]))\n",
    "reformat[ii,22]=reformat[ii,23]*100\n",
    "\n",
    "\n",
    "# BP\n",
    "ii=(~np.isnan(reformat[:,8])) & (~np.isnan(reformat[:,9])) & (np.isnan(reformat[:,10]))\n",
    "reformat[ii,10]=(3*reformat[ii,9]-reformat[ii,8])/2\n",
    "ii=(~np.isnan(reformat[:,8])) & (~np.isnan(reformat[:,10])) & (np.isnan(reformat[:,9]))\n",
    "reformat[ii,9]=(reformat[ii,8]+2*reformat[ii,10])/3\n",
    "ii=(~np.isnan(reformat[:,9])) & (~np.isnan(reformat[:,10])) & (np.isnan(reformat[:,8]))\n",
    "reformat[ii,8]=3*reformat[ii,9]-2*reformat[ii,10]\n",
    "\n",
    "#TEMP\n",
    "#some values recorded in the wrong column\n",
    "ii=(reformat[:,14]>25) & (reformat[:,14]<45) #tempF close to 37deg??!\n",
    "reformat[ii,13]=reformat[ii,14]\n",
    "reformat[ii,14]=np.nan\n",
    "ii=reformat[:,13]>70  #tempC > 70?!!! probably degF\n",
    "reformat[ii,14]=reformat[ii,13]\n",
    "reformat[ii,13]=np.nan\n",
    "ii=(~np.isnan(reformat[:,13])) & (np.isnan(reformat[:,14]))\n",
    "reformat[ii,14]=reformat[ii,13]*1.8+32\n",
    "ii=(~np.isnan(reformat[:,14])) & (np.isnan(reformat[:,13]))\n",
    "reformat[ii,13]=(reformat[ii,14]-32)/1.8\n",
    "\n",
    "\n",
    "# Hb/Ht\n",
    "ii=(~np.isnan(reformat[:,49])) & (np.isnan(reformat[:,50]))\n",
    "reformat[ii,50]=(reformat[ii,49]*2.862)+1.216\n",
    "ii=(~np.isnan(reformat[:,50])) & (np.isnan(reformat[:,49]))\n",
    "reformat[ii,49]=(reformat[ii,50]-1.216)/2.862\n",
    "\n",
    "#BILI\n",
    "ii=(~np.isnan(reformat[:,43])) & (np.isnan(reformat[:,44]))\n",
    "reformat[ii,44]=(reformat[ii,43]*0.6934)-0.1752\n",
    "ii=(~np.isnan(reformat[:,44])) & (np.isnan(reformat[:,43]))\n",
    "reformat[ii,43]=(reformat[ii,44]+0.1752)/0.6934\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "########################################################################\n",
    "#                      SAMPLE AND HOLD on RAW DATA\n",
    "#########################################################################\n",
    "\n",
    "reformat=SAH(reformat[:,0:68],sample_and_hold)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#########################################################################\n",
    "#                             DATA COMBINATION\n",
    "#########################################################################\n",
    "\n",
    "#WARNING: the time window of interest has been defined above (here -48 -> +24)! \n",
    "\n",
    "timestep=4  # resolution of timesteps, in hours\n",
    "irow=0\n",
    "icustayidlist=np.unique(reformat[:,1].astype('int64'))\n",
    "reformat2=np.full((reformat.shape[0],85),np.nan)  # output array\n",
    "npt=icustayidlist.size  # number of patients\n",
    "# Adding 2 empty cols for future shock index=HR/SBP and P/F\n",
    "reformat= np.insert(reformat,68,np.nan,axis=1)\n",
    "reformat= np.insert(reformat,69,np.nan,axis=1)\n",
    "\n",
    "\n",
    "for i in range(npt): \n",
    "    \n",
    "    if(i%10000==0):\n",
    "        print(i)\n",
    "        \n",
    "    icustayid=icustayidlist[i]  #1 to 100000, NOT 200 to 300K!\n",
    "     \n",
    "    # CHARTEVENTS AND LAB VALUES\n",
    "    temp=reformat[reformat[:,1]==icustayid,:]   #subtable of interest\n",
    "    beg=temp[0,2]   #timestamp of first record\n",
    "    \n",
    "    # IV FLUID STUFF\n",
    "    iv=np.where((inputMV[:,0]==icustayid+200000))[0]   #rows of interest in inputMV\n",
    "    input=inputMV[iv,:]    # subset of interest\n",
    "    iv=np.where((inputCV[:,0]==icustayid+200000))[0]   #rows of interest in inputCV\n",
    "    input2=inputCV[iv,:]    #subset of interest\n",
    "    startt=input[:,1] # start of all infusions and boluses\n",
    "    endt=input[:,2] # end of all infusions and boluses\n",
    "    rate=input[:,7]  #rate of infusion (is NaN for boluses) || corrected for tonicity\n",
    "        \n",
    "    pread=inputpreadm[inputpreadm[:,0]==icustayid+200000,1] #preadmission volume\n",
    "    if(pread.size!=0):              # store the value, if available\n",
    "        totvol=np.nansum(pread)\n",
    "    else: \n",
    "        totvol=0   # if not documented: it's zero\n",
    "       \n",
    "    # compute volume of fluid given before start of record!!!\n",
    "    t0=0\n",
    "    t1=beg\n",
    "    #input from MV (4 ways to compute)\n",
    "    \n",
    "    infu = np.nansum(rate*(endt-startt)*((endt<=t1)&(startt>=t0))/3600   +    rate*(endt-t0)*((startt<=t0)&(endt<=t1)&(endt>=t0))/3600 +     rate*(t1-startt)*((startt>=t0)&(endt>=t1)&(startt<=t1))/3600 +      rate*(t1-t0)*((endt>=t1)&(startt<=t0))   /3600) \n",
    "    # all boluses received during this timestep, from inputMV (need to check rate is NaN) and inputCV (simpler):\n",
    "    bolus=np.nansum(input[(np.isnan(input[:,5])) & (input[:,1]>=t0) & (input[:,1]<=t1),6]) + np.nansum(input2[(input2[:,1]>=t0) & (input2[:,1]<=t1),4])  \n",
    "    totvol=np.nansum(np.array([totvol,infu,bolus])) \n",
    "    \n",
    "\n",
    "    # VASOPRESSORS    \n",
    "    iv=np.where(vasoMV[:,0]==icustayid+200000)[0]   #rows of interest in vasoMV\n",
    "    vaso1=vasoMV[iv,:]    # subset of interest\n",
    "    iv=np.where(vasoCV[:,0]==icustayid+200000)[0]   #rows of interest in vasoCV\n",
    "    vaso2=vasoCV[iv,:]    # subset of interest\n",
    "    startv=vaso1[:,2]     # start of VP infusion\n",
    "    endv=vaso1[:,3]       # end of VP infusions\n",
    "    ratev=vaso1[:,4]      #rate of VP infusion\n",
    "            \n",
    "   \n",
    "    # DEMOGRAPHICS / gender, age, elixhauser, re-admit, died in hosp?, died within\n",
    "    # 48h of out_time (likely in ICU or soon after), died within 90d after admission?        \n",
    "    demogi=np.where(demog.loc[:,'icustay_id']==icustayid+200000)[0]  \n",
    "    dem=np.array([ demog.loc[demogi,'gender'], demog.loc[demogi,'age'],demog.loc[demogi,'elixhauser'], demog.loc[demogi,'adm_order']>1,demog.loc[demogi,'morta_hosp'], np.abs(demog.loc[demogi,'dod']-demog.loc[demogi,'outtime'])<(24*3600*2), demog.loc[demogi,'morta_90'], (qstime[icustayid-1,3]-qstime[icustayid-1,2])/3600 ])     \n",
    "    \n",
    "\n",
    "    # URINE OUTPUT\n",
    "    iu=np.where(UO[:,0]==icustayid+200000)[0] #rows of interest in inputMV\n",
    "    output=UO[iu,:]    # subset of interest\n",
    "    pread=UOpreadm[UOpreadm[:,0]==icustayid,3] # preadmission UO\n",
    "    if pread.size!=0:     #store the value, if available\n",
    "        UOtot=np.nansum(pread)\n",
    "    else :\n",
    "        UOtot=0;\n",
    "    \n",
    "    # adding the volume of urine produced before start of recording!    \n",
    "    UOnow=np.nansum(output[(output[:,1]>=t0) & (output[:,1]<=t1),3]) #t0 and t1 defined above\n",
    "    UOtot=np.nansum(np.array([UOtot,UOnow])) \n",
    "    \n",
    "    \n",
    "    for j in range(0,80,timestep): # -28 until +52 = 80 hours in total\n",
    "        t0=3600*j+ beg # left limit of time window\n",
    "        t1=3600*(j+timestep)+beg # right limit of time window\n",
    "        ii=(temp[:,2]>=t0) & (temp[:,2] <=t1)  #index of items in this time period\n",
    "        \n",
    "        if np.sum(ii)>0 : \n",
    "            \n",
    "            \n",
    "            #ICUSTAY_ID, OUTCOMES, DEMOGRAPHICS\n",
    "            reformat2[irow,0]=(j/timestep)+1   #'bloc' = timestep (1,2,3...)\n",
    "            reformat2[irow,1]=icustayid        #icustay_ID\n",
    "            reformat2[irow,2]=3600*j+ beg      #t0 = lower limit of time window\n",
    "            reformat2[irow,3:11]=dem           #demographics and outcomes\n",
    "            \n",
    "        \n",
    "            #CHARTEVENTS and LAB VALUES (+ includes empty cols for shock index and P/F)\n",
    "            value=temp[ii,:]#records all values in this timestep\n",
    " \n",
    "        \n",
    "              # #####################   DISCUSS ADDING STUFF HERE / RANGE, MIN, MAX ETC   ################\n",
    "        \n",
    "            if np.sum(ii)==1:   #if only 1 row of values at this timestep\n",
    "                reformat2[irow,11:78]=value[:,3:]\n",
    "            else : \n",
    "                reformat2[irow,11:78]=np.nanmean(value[:,3:],axis=0) #mean of all available values\n",
    "        \n",
    "            #VASOPRESSORS\n",
    "                # for CV: dose at timestamps.\n",
    "                # for MV: 4 possibles cases, each one needing a different way to compute the dose of VP actually administered:\n",
    "                #----t0---start----end-----t1----\n",
    "                #----start---t0----end----t1----\n",
    "                #-----t0---start---t1---end\n",
    "                #----start---t0----t1---end----\n",
    "\n",
    "        \n",
    "            #MV\n",
    "            v =((endv>=t0) & (endv<=t1)) | ((startv>=t0) & (endv<=t1)) | ((startv>=t0) & (startv<=t1)) | ((startv<=t0) & (endv>=t1))\n",
    "\n",
    "            #CV\n",
    "            v2=vaso2[(vaso2[:,2]>=t0) & (vaso2[:,2]<=t1),3] \n",
    "            \n",
    "            temp_list = [] \n",
    "            if(ratev[v].size!=0):\n",
    "                temp_list.append(ratev[v].reshape(-1,1))\n",
    "            if(v2.size!=0):\n",
    "                temp_list.append(v2.reshape(-1,1))\n",
    "                \n",
    "            if(len(temp_list)!=0):\n",
    "                rv = np.vstack(temp_list)\n",
    "            else : \n",
    "                rv = np.array([])\n",
    "            \n",
    "            v1=np.nanmedian(rv)\n",
    "\n",
    "            if(rv.size!=0):\n",
    "                v2=np.nanmax(rv)\n",
    "            else:\n",
    "                v2 =np.array([]) \n",
    "              \n",
    "            if v1.size!=0 and ~np.isnan(v1) and v2.size!=0 and ~np.isnan(v2): \n",
    "                reformat2[irow,78]=v1    #median of dose of VP\n",
    "                reformat2[irow,79]=v2    #max dose of VP\n",
    "        \n",
    "            #INPUT FLUID\n",
    "            #input from MV (4 ways to compute)\n",
    "            infu = np.nansum(rate*(endt-startt)*((endt<=t1)&(startt>=t0))/3600   +    rate*(endt-t0)*((startt<=t0)&(endt<=t1)&(endt>=t0))/3600 +     rate*(t1-startt)*((startt>=t0)&(endt>=t1)&(startt<=t1))/3600 +      rate*(t1-t0)*((endt>=t1)&(startt<=t0))   /3600) \n",
    "            #all boluses received during this timestep, from inputMV (need to check rate is NaN) and inputCV (simpler):\n",
    "            bolus=np.nansum(input[(np.isnan(input[:,5])) & (input[:,1]>=t0) & (input[:,1]<=t1),6]) + np.nansum(input2[(input2[:,1]>=t0) & (input2[:,1]<=t1),4])  \n",
    "                \n",
    "            #sum fluid given\n",
    "            totvol=np.nansum(np.array([totvol,infu,bolus]))\n",
    "            \n",
    "            reformat2[irow,80]=totvol    #total fluid given\n",
    "            reformat2[irow,81]=np.nansum(np.array([infu,bolus]))   #fluid given at this step\n",
    "        \n",
    "            #UO\n",
    "            UOnow=np.nansum(output[(output[:,1]>=t0) & (output[:,1]<=t1),3])  \n",
    "            UOtot=np.nansum(np.array([UOtot, UOnow]))\n",
    "            reformat2[irow,82]=UOtot    #total UO\n",
    "            reformat2[irow,83]=np.nansum(UOnow)   #UO at this step\n",
    "\n",
    "            #CUMULATED BALANCE\n",
    "            reformat2[irow,84]=totvol-UOtot    #cumulated balance\n",
    "\n",
    "            irow=irow+1 \n",
    "\n",
    "\n",
    "reformat2 = reformat2[:irow,:]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#########################################################################\n",
    "#             CONVERT TO TABLE AND KEEP ONLY WANTED VARIABLE\n",
    "#########################################################################\n",
    "\n",
    "# dataheaders \n",
    "dataheaders=sample_and_hold[0,:].tolist()+['Shock_Index', 'PaO2_FiO2']\n",
    "dataheaders = ['bloc','icustayid','charttime','gender','age','elixhauser','re_admission', 'died_in_hosp', 'died_within_48h_of_out_time','mortality_90d','delay_end_of_record_and_discharge_or_death']+dataheaders\n",
    "dataheaders = dataheaders+  [ 'median_dose_vaso','max_dose_vaso','input_total','input_4hourly','output_total','output_4hourly','cumulated_balance']\n",
    "\n",
    "reformat2t=pd.DataFrame(reformat2.copy(),columns = dataheaders) \n",
    "\n",
    "# headers I want to keep\n",
    "dataheaders5 = ['bloc','icustayid','charttime','gender','age','elixhauser','re_admission', 'died_in_hosp', 'died_within_48h_of_out_time','mortality_90d','delay_end_of_record_and_discharge_or_death','SOFA','SIRS',\\\n",
    "    'Weight_kg','GCS','HR','SysBP','MeanBP','DiaBP','RR','SpO2','Temp_C','FiO2_1','Potassium','Sodium','Chloride','Glucose',\\\n",
    "    'BUN','Creatinine','Magnesium','Calcium','Ionised_Ca','CO2_mEqL','SGOT','SGPT','Total_bili','Albumin','Hb','WBC_count','Platelets_count','PTT','PT','INR',\\\n",
    "    'Arterial_pH','paO2','paCO2','Arterial_BE','HCO3','Arterial_lactate','mechvent','Shock_Index','PaO2_FiO2',\\\n",
    "    'median_dose_vaso','max_dose_vaso','input_total','input_4hourly','output_total','output_4hourly','cumulated_balance']\n",
    "\n",
    "ii=np.where(np.isin(reformat2t.columns,dataheaders5))[0] \n",
    "reformat3t = reformat2t.iloc[:,ii] \n",
    "\n",
    "## SOME DATA MANIP BEFORE IMPUTATION\n",
    "\n",
    "# CORRECT GENDER\n",
    "reformat3t.loc[:,'gender']=reformat3t.loc[:,'gender']-1 \n",
    "\n",
    "# CORRECT AGE > 200 yo\n",
    "ii=reformat3t.loc[:,'age']>150*365.25\n",
    "reformat3t.loc[ii,'age']=91.4*365.25\n",
    "\n",
    "# FIX MECHVENT\n",
    "ii = np.isnan(reformat3t.loc[:,'mechvent'])\n",
    "reformat3t.loc[ii,'mechvent']=0\n",
    "ii = reformat3t.loc[:,'mechvent']>0\n",
    "reformat3t.loc[ii,'mechvent']=1\n",
    "\n",
    "# FIX Elixhauser missing values\n",
    "ii = np.isnan(reformat3t.loc[:,'elixhauser'])\n",
    "reformat3t.loc[ii,'elixhauser']=np.nanmedian(reformat3t.loc[:,'elixhauser'])  #use the median value / only a few missing data points \n",
    "\n",
    "\n",
    "# vasopressors / no NAN\n",
    "ii=np.isnan(reformat3t.loc[:,'median_dose_vaso'])\n",
    "reformat3t.loc[ii,'median_dose_vaso']=np.zeros((np.sum(ii)))\n",
    "ii=np.isnan(reformat3t.loc[:,'max_dose_vaso'])\n",
    "reformat3t.loc[ii,'max_dose_vaso']=np.zeros((np.sum(ii)))\n",
    "\n",
    "# check prop of missingness here\n",
    "miss=np.sum(np.isnan(reformat3t)/reformat3t.shape[0]) \n",
    "plt.bar(np.arange(miss.shape[0]), miss)\n",
    "# plt.show() \n",
    "\n",
    "# I fill the values temporarily with zeros, otherwise kNN imp doesnt work\n",
    "reformat3t['Shock_Index']=np.zeros((reformat3t.shape[0],1))\n",
    "reformat3t['PaO2_FiO2']=np.zeros((reformat3t.shape[0],1))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#########################################################################\n",
    "#             HANDLING OF MISSING VALUES  &  CREATE REFORMAT4T\n",
    "#########################################################################\n",
    "\n",
    "def fixgaps(x): \n",
    "    # FIXGAPS Linearly interpolates gaps in a time series\n",
    "    # YOUT=FIXGAPS(YIN) linearly interpolates over NaN\n",
    "    # in the input time series (may be complex), but ignores\n",
    "    # trailing and leading NaN.\n",
    "    # R. Pawlowicz 6/Nov/99\n",
    "    \n",
    "    y=x\n",
    "    bd = np.isnan(x)\n",
    "    gd = np.where(~bd)[0]\n",
    "    \n",
    "    bd[0:min(gd)]=0 \n",
    "    bd[max(gd)+1:]=0 \n",
    "    \n",
    "    y[bd] = interp1d(gd,x[gd])(np.where(bd)[0])\n",
    "                          \n",
    "    return y \n",
    "\n",
    "\n",
    "\n",
    "# K=1\n",
    "# distance = seuclidean \n",
    "# Reference: matlab's knnimpute.m code \n",
    "# Reference: https://github.com/ogeidix/kddcup09/blob/master/utilities/knnimpute.m\n",
    "\n",
    "def knnimpute(data):\n",
    "    \n",
    "    K=1 \n",
    "    userWeights = False \n",
    "    useWMean = True  \n",
    "    # create a copy of data for output\n",
    "    imputed = data.copy() \n",
    "    \n",
    "    # identify missing vals\n",
    "    nanVals = np.isnan(data)\n",
    "    \n",
    "    # use rows without nans for calculation of nearest neighbors\n",
    "    noNans = (np.sum(nanVals,axis=1) == 0)    \n",
    "    dataNoNans = data[noNans,:] \n",
    "    \n",
    "    distances = pdist(np.transpose(dataNoNans),'seuclidean')    \n",
    "    \n",
    "    SqF = squareform(distances)\n",
    "    \n",
    "    temp = SqF - np.identity(SqF.shape[0])\n",
    "    \n",
    "    dists= np.transpose(np.sort(temp))\n",
    "    \n",
    "    ndx = np.transpose(np.argsort(temp,kind='stable'))\n",
    "    \n",
    "    equalDists = np.vstack([np.diff(dists[1:,:],axis=0)==0.0,np.full(dists.shape[1],False)])\n",
    "    \n",
    "    rows= np.where(np.transpose(nanVals))[1]         \n",
    "    cols = np.where(np.transpose(nanVals))[0] \n",
    "\n",
    "    for count in range(rows.size):\n",
    "        for nearest in range(1,ndx.shape[0]-K+1):\n",
    "            L = np.where(equalDists[nearest+K-2:,cols[count]]==0)[0][0] \n",
    "            dataVals = data[rows[count],ndx[nearest:nearest+K+L,cols[count]]]\n",
    "            if(useWMean):\n",
    "                if(~userWeights):\n",
    "                    weights = 1/dists[1:K+L+1,cols[count]]\n",
    "                val = wnanmean(dataVals, weights)\n",
    "            if(~np.isnan(val)):\n",
    "                imputed[rows[count],cols[count]] = val \n",
    "                break \n",
    "    return imputed \n",
    "                \n",
    "def wnanmean(x,weights):\n",
    "    x = x.copy()\n",
    "    weights = weights.copy() \n",
    "    nans = np.isnan(x)\n",
    "    infs = np.isinf(weights)\n",
    "\n",
    "    if all(nans): \n",
    "        return np.nan\n",
    "    if any(infs):\n",
    "        return np.nanmean(x[infs])\n",
    "    \n",
    "    x[nans] = 0 \n",
    "    weights[nans] = 0 \n",
    "    weights = weights/np.sum(weights)\n",
    "    return np.dot(np.transpose(weights),x)\n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "# Do linear interpol where missingness is low (kNN imputation doesnt work if all rows have missing values)\n",
    "\n",
    "reformat3=reformat3t.values.copy() \n",
    "miss=(np.sum(np.isnan(reformat3),axis=0)/reformat3.shape[0]) \n",
    "\n",
    "ii=(miss>0) & (miss<0.05) # less than 5% missingness\n",
    "\n",
    "mechventcol = reformat3t.columns.get_loc('mechvent')\n",
    "\n",
    "for i in range(10,mechventcol): #correct col by col, otherwise it does it wrongly\n",
    "    if ii[i]==1: \n",
    "        reformat3[:,i]=fixgaps(reformat3[:,i])\n",
    "\n",
    "        \n",
    "reformat3t.iloc[:,10:mechventcol]=reformat3[:,10:mechventcol]\n",
    "\n",
    "\n",
    "# KNN IMPUTATION -  Done on chunks of 10K records.\n",
    "\n",
    "ref=reformat3[:,10:mechventcol].copy()  #columns of interest\n",
    "\n",
    "\n",
    "for i in range(0,(reformat3.shape[0]-9999),10000):   #dataset divided in 5K rows chunks (otherwise too large)\n",
    "    print(i)\n",
    "    ref[i:i+10000,:]=np.transpose(knnimpute(np.transpose(ref[i:i+10000,:]))) \n",
    "\n",
    "ref[-10000:,:]=np.transpose(knnimpute(np.transpose(ref[-10000:,:]))) #the last bit is imputed from the last 10K rows\n",
    "\n",
    "\n",
    "# I paste the data interpolated, but not the demographics and the treatments\n",
    "reformat3t.iloc[:,10:mechventcol]=ref   \n",
    "\n",
    "reformat4t=reformat3t.copy() \n",
    "reformat4=reformat4t.values.copy() \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##########################################################################\n",
    "#        COMPUTE SOME DERIVED VARIABLES: P/F, Shock Index, SOFA, SIRS...\n",
    "##########################################################################\n",
    "\n",
    "\n",
    "# re-compute P/F with no missing values...\n",
    "p=reformat4t.columns.get_loc('paO2') \n",
    "f=reformat4t.columns.get_loc('FiO2_1') \n",
    "reformat4t.loc[:,'PaO2_FiO2']=reformat4[:,p]/reformat4[:,f]   \n",
    "\n",
    "# recompute SHOCK INDEX without NAN and INF\n",
    "p=reformat4t.columns.get_loc('HR') \n",
    "f=reformat4t.columns.get_loc('SysBP') \n",
    "a=reformat4t.columns.get_loc('Shock_Index') \n",
    "reformat4[:,a]=reformat4[:,p]/reformat4[:,f]  \n",
    "ii = np.isinf(reformat4[:,a])\n",
    "reformat4[ii,a]=np.nan \n",
    "d=np.nanmean(reformat4[:,a])\n",
    "ii = np.isnan(reformat4[:,a])\n",
    "reformat4[ii,a]=d  #replace NaN with average value ~ 0.8\n",
    "reformat4t.loc[:,'Shock_Index']=reformat4[:,a]\n",
    "# Reference: https://stackoverflow.com/questions/13733034/equavelent-python-command-for-quantile-in-matlab \n",
    "ii=reformat4t['Shock_Index']>=mquantiles(reformat4t['Shock_Index'], 0.999, alphap=0.5, betap=0.5)[0] # replace outliers with 99.9th percentile\n",
    "reformat4t.loc[ii,'Shock_Index']=mquantiles(reformat4t['Shock_Index'],0.999,alphap=0.5, betap=0.5)[0]\n",
    "\n",
    "\n",
    "# SOFA - at each timepoint\n",
    "# need (in this order):  P/F  MV  PLT  TOT_BILI  MAP  NORAD(max)  GCS  CR  UO\n",
    "\n",
    "s=reformat4t.loc[:,['PaO2_FiO2','Platelets_count','Total_bili','MeanBP','max_dose_vaso','GCS','Creatinine','output_4hourly']].values  \n",
    "p=np.array([0, 1, 2, 3, 4]) \n",
    "\n",
    "s1=np.transpose(np.array([s[:,0]>400, (s[:,0]>=300) & (s[:,0]<400), (s[:,0]>=200) & (s[:,0]<300), (s[:,0]>=100) & (s[:,0]<200), s[:,0]<100 ]))  #count of points for all 6 criteria of sofa\n",
    "s2=np.transpose(np.array([s[:,1]>150, (s[:,1]>=100) & (s[:,1]<150), (s[:,1]>=50) & (s[:,1]<100), (s[:,1]>=20) & (s[:,1]<50), s[:,1]<20 ]))\n",
    "s3=np.transpose(np.array([s[:,2]<1.2, (s[:,2]>=1.2) & (s[:,2]<2), (s[:,2]>=2) &(s[:,2]<6), (s[:,2]>=6) & (s[:,2]<12), s[:,2]>12 ]))\n",
    "s4=np.transpose(np.array([s[:,3]>=70, (s[:,3]<70) & (s[:,3]>=65), s[:,3]<65, (s[:,4]>0) & (s[:,4]<=0.1), s[:,4]>0.1 ]))\n",
    "s5=np.transpose(np.array([s[:,5]>14, (s[:,5]>12) & (s[:,5]<=14), (s[:,5]>9) & (s[:,5]<=12), (s[:,5]>5) & (s[:,5]<=9), s[:,5]<=5 ]))\n",
    "s6=np.transpose(np.array([s[:,6]<1.2, (s[:,6]>=1.2) & (s[:,6]<2), (s[:,6]>=2) & (s[:,6]<3.5), ((s[:,6]>=3.5) & (s[:,6]<5)) | (s[:,7]<84), (s[:,6]>5) | (s[:,7]<34) ]))\n",
    "\n",
    "nrcol=reformat4.shape[1]   # nr of variables in data\n",
    "reformat4 = np.hstack([reformat4,np.zeros((reformat4.shape[0],7))])\n",
    "\n",
    "for i in range(0,reformat4.shape[0]):   \n",
    "    p_s1 = p[s1[i,:]]\n",
    "    p_s2 = p[s2[i,:]]\n",
    "    p_s3 = p[s3[i,:]]\n",
    "    p_s4 = p[s4[i,:]]\n",
    "    p_s5 = p[s5[i,:]]\n",
    "    p_s6 = p[s6[i,:]] \n",
    "    \n",
    "    if(p_s1.size==0 or p_s2.size==0 or p_s3.size==0 or p_s4.size==0 or p_s5.size==0 or p_s6.size==0):\n",
    "        t = 0 \n",
    "    else : \n",
    "        t=max(p_s1)+max(p_s2)+max(p_s3)+max(p_s4)+max(p_s5)+max(p_s6)  #SUM OF ALL 6 CRITERIA        \n",
    "    \n",
    "    if(t):\n",
    "        reformat4[i,nrcol:nrcol+7]=np.array([max(p_s1), max(p_s2), max(p_s3), max(p_s4), max(p_s5), max(p_s6), t])\n",
    "        \n",
    "        \n",
    "# SIRS - at each timepoint |  need: temp HR RR PaCO2 WBC \n",
    "s=reformat4t.loc[:,['Temp_C','HR','RR','paCO2','WBC_count']].values  #indices of vars used in SOFA\n",
    "\n",
    "\n",
    "s1=np.transpose(np.array([(s[:,0]>=38) | (s[:,0]<=36) ],dtype=np.int64))   #count of points for all criteria of SIRS\n",
    "s2=np.transpose(np.array([s[:,1]>90 ],dtype=np.int64)) \n",
    "s3=np.transpose(np.array([(s[:,2]>=20)|(s[:,3]<=32)],dtype=np.int64))\n",
    "s4=np.transpose(np.array([(s[:,4]>=12)|(s[:,4]<4)],dtype=np.int64))\n",
    "\n",
    "reformat4 = np.insert(reformat4,nrcol+7,0,axis=1)\n",
    "reformat4[:,nrcol+7]=np.transpose(s1+s2+s3+s4)\n",
    "\n",
    "# adds 2 cols for SOFA and SIRS, if necessary\n",
    "if 'SIRS' not in reformat4t.columns : \n",
    "    reformat4t['SOFA'] = 0 \n",
    "    reformat4t['SIRS'] = 0 \n",
    "    \n",
    "# more IO corrections\n",
    "ii=reformat4t['input_total']<0\n",
    "reformat4t.loc[ii,'input_total']=0\n",
    "ii=reformat4t['input_4hourly']<0\n",
    "reformat4t.loc[ii,'input_4hourly']=0\n",
    "\n",
    "# records values\n",
    "reformat4t.loc[:,'SOFA']=reformat4[:,-2]\n",
    "reformat4t.loc[:,'SIRS']=reformat4[:,-1]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#########################################################################\n",
    "#                     CREATE FINAL MIMIC_TABLE\n",
    "#########################################################################\n",
    "\n",
    "\n",
    "MIMICtable = reformat4t\n",
    "\n",
    "MIMICtable.to_csv('mimictable.csv',index=False,na_rep='NaN')  \n",
    "\n",
    "# Save pickle \n",
    "with open('step_4_start.pkl', 'wb') as file:\n",
    "     pickle.dump(MIMICtable,file)"
   ]
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}