sepsis_def.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",
    "\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",
    "# Matlab to Python Reference \n",
    "# Reference: https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html \n",
    "# Reference: http://mathesaurus.sourceforge.net/matlab-numpy.html\n",
    "#\n",
    "# You need to understand how to handle the pandas index \n",
    "# https://towardsdatascience.com/pandas-index-explained-b131beaf6f7b\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Github Reference: 0669f8907e65503641857ca76aa46938641e513f \n",
    "# Reference File Name: AIClinician_sepsis3_def_160219.py \n",
    "# Source: https://github.com/matthieukomorowski/AI_Clinician\n",
    "\n",
    "## AI Clinician Identifiying MIMIC-III sepsis cohort\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",
    "# IDENTIFIES THE COHORT OF PATIENTS WITH SEPSIS in MIMIC-III\n",
    "\n",
    "# PURPOSE:\n",
    "# ------------------------------\n",
    "# This creates a list of icustayIDs of patients who develop sepsis at some point \n",
    "# in the ICU. records charttime for onset of sepsis. Uses sepsis3 criteria\n",
    "\n",
    "# STEPS:\n",
    "# -------------------------------\n",
    "# IMPORT DATA FROM CSV FILES\n",
    "# FLAG PRESUMED INFECTION\n",
    "# PREPROCESSING\n",
    "# REFORMAT in 4h time slots\n",
    "# COMPUTE SOFA at each time step\n",
    "# FLAG SEPSIS\n",
    "\n",
    "# note: the process generates the same features as the final MDP dataset, most of which are not used to compute SOFA\n",
    "# External files required: Reflabs, Refvitals, sample_and_hold (all saved in reference_matrices.mat file)\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",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "import numpy as np \n",
    "import scipy.io as sio \n",
    "import pandas as pd \n",
    "import csv \n",
    "import math\n",
    "import pickle \n",
    "from tqdm import tqdm\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 Reference_Matrices.mat \n",
    "mat_data = sio.loadmat('reference_matrices.mat')\n",
    "Reflabs = mat_data['Reflabs'] # numpy.ndarray \n",
    "Refvitals = mat_data['Refvitals'] # numpy.ndarray \n",
    "sample_and_hold = mat_data['sample_and_hold'] # numpy.ndarray \n",
    "\n",
    "# preprocessing sample_and_hold \n",
    "for index,s in enumerate(sample_and_hold[0,:]):\n",
    "    sample_and_hold[0,index] = s[0].replace('\\'','') \n",
    "\n",
    "for index,s in enumerate(sample_and_hold[1,:]):\n",
    "    sample_and_hold[1,index] = s[0][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ########################################################################\n",
    "# IMPORT ALL DATA\n",
    "\n",
    "# Based on whos command in Matlab, data type was double (8 bytes) except demog (Table)\n",
    "# In Python, it was stored in float64 \n",
    "\n",
    "print('Load abx')\n",
    "abx = pd.read_csv('D:/exportdir/abx.csv',delimiter='|').values\n",
    "print('Load culture')\n",
    "culture=pd.read_csv('D:/exportdir/culture.csv',delimiter='|').values\n",
    "print('Load microbio')\n",
    "microbio=pd.read_csv('D:/exportdir/microbio.csv',delimiter='|').values\n",
    "print('Load demog')\n",
    "demog=pd.read_csv('D:/exportdir/demog.csv',delimiter='|') # read as DataFrame \n",
    "print('Load ce010')\n",
    "ce010=pd.read_csv('D:/exportdir/ce010000.csv',delimiter='|').values\n",
    "print('Load ce1020')\n",
    "ce1020=pd.read_csv('D:/exportdir/ce1000020000.csv',delimiter='|').values\n",
    "print('Load ce2030')\n",
    "ce2030=pd.read_csv('D:/exportdir/ce2000030000.csv',delimiter='|').values\n",
    "print('Load ce3040')\n",
    "ce3040=pd.read_csv('D:/exportdir/ce3000040000.csv',delimiter='|').values\n",
    "print('Load ce4050')\n",
    "ce4050=pd.read_csv('D:/exportdir/ce4000050000.csv',delimiter='|').values\n",
    "print('Load ce5060')\n",
    "ce5060=pd.read_csv('D:/exportdir/ce5000060000.csv',delimiter='|').values\n",
    "print('Load ce6070')\n",
    "ce6070=pd.read_csv('D:/exportdir/ce6000070000.csv',delimiter='|').values\n",
    "print('Load ce7080')\n",
    "ce7080=pd.read_csv('D:/exportdir/ce7000080000.csv',delimiter='|').values\n",
    "print('Load ce8090')\n",
    "ce8090=pd.read_csv('D:/exportdir/ce8000090000.csv',delimiter='|').values\n",
    "print('Load ce90100') \n",
    "ce90100=pd.read_csv('D:/exportdir/ce90000100000.csv',delimiter='|').values\n",
    "print('Load labU') \n",
    "labU = np.vstack([pd.read_csv('D:/exportdir/labs_ce.csv',delimiter='|').values, pd.read_csv('D:/exportdir/labs_le.csv',delimiter='|').values  ])\n",
    "print('Load MV')\n",
    "MV=pd.read_csv('D:/exportdir/mechvent.csv',delimiter='|').values\n",
    "print('Load inputpreadm')\n",
    "inputpreadm=pd.read_csv('D:/exportdir/preadm_fluid.csv',delimiter='|').values\n",
    "print('Load inputMV')\n",
    "inputMV=pd.read_csv('D:/exportdir/fluid_mv.csv',delimiter='|').values\n",
    "print('Load inputCV')\n",
    "inputCV=pd.read_csv('D:/exportdir/fluid_cv.csv',delimiter='|').values\n",
    "print('Load vasoMV')\n",
    "vasoMV=pd.read_csv('D:/exportdir/vaso_mv.csv',delimiter='|').values\n",
    "print('Load vasoCV')\n",
    "vasoCV=pd.read_csv('D:/exportdir/vaso_cv.csv',delimiter='|').values\n",
    "print('Load UOpreadm')\n",
    "UOpreadm=pd.read_csv('D:/exportdir/preadm_uo.csv',delimiter='|').values\n",
    "print('Load UO') \n",
    "UO=pd.read_csv('D:/exportdir/uo.csv',delimiter='|').values\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ########################################################################\n",
    "#                       INITIAL DATA MANIPULATIONS\n",
    "# #########################################################################\n",
    "# Be Careful: Matlab Index starts from 1 but Python Index starts from 0  \n",
    "# Numpy Cheat Sheet: http://taewan.kim/post/numpy_cheat_sheet/ \n",
    "# Numpy *, multiply(), dot() difference => https://www.tutorialexample.com/understand-numpy-np-multiply-np-dot-and-operation-a-beginner-guide-numpy-tutorial/\n",
    "\n",
    "\n",
    "ii = np.isnan(microbio[:,2]) # if charttime is empty but chartdate isn't\n",
    "microbio[ii,2] = microbio[ii,3]  #copy time\n",
    "microbio = np.delete(microbio,3,1) # delete chardate\n",
    "# Add empty col in microbio (# 3 and #5)\n",
    "microbio = np.insert(microbio,2,0,axis=1)\n",
    "microbio = np.insert(microbio,4,0,axis=1)\n",
    "# Combine both tables for micro events\n",
    "bacterio = np.vstack([microbio, culture])\n",
    "# correct NaNs in DEMOG\n",
    "demog.loc[np.isnan(demog.loc[:,'morta_90']),'morta_90']=0 \n",
    "demog.loc[np.isnan(demog.loc[:,'morta_hosp']),'morta_hosp']=0\n",
    "demog.loc[np.isnan(demog.loc[:,'elixhauser']),'elixhauser']=0\n",
    "# compute normalized rate of infusion\n",
    "# if we give 100 ml of hypertonic fluid (600 mosm/l) at 100 ml/h (given in 1h) it is 200 ml of NS equivalent\n",
    "# so the normalized rate of infusion is 200 ml/h (different volume in same duration)\n",
    "inputMV = np.insert(inputMV,7,np.nan,axis=1) # Initialize the new column with nan \n",
    "ii = inputMV[:,4]!=0 # to avoid divide by zero \n",
    "inputMV[ii,7] = inputMV[ii,6]*inputMV[ii,5]/inputMV[ii,4] \n",
    "\n",
    "# fill-in missing ICUSTAY IDs in bacterio\n",
    "for i in range(bacterio.shape[0]):\n",
    "    if(i%100000==0) :\n",
    "        print(i) \n",
    "    if(bacterio[i,2]==0): # if missing icustayid\n",
    "        o = bacterio[i,3] # charttime\n",
    "        subjectid = bacterio[i,0]\n",
    "        hadmid = bacterio[i,1]\n",
    "        ii = np.where(demog.loc[:,'subject_id']==subjectid)[0]\n",
    "        jj = np.where((demog.loc[:,'subject_id']==subjectid) & (demog.loc[:, 'hadm_id']==hadmid))[0]\n",
    "        for j in range(len(ii)): \n",
    "            if (o >= demog.loc[ii[j],'intime']-48*3600) and (o <= demog.loc[ii[j],'outtime']+48*3600):\n",
    "                bacterio[i,2] = demog.loc[ii[j],'icustay_id']\n",
    "            elif len(ii) == 1 : # if we cant confirm from admission and discharge time but there is only 1 admission: it's the one!!\n",
    "                bacterio[i,2] = demog.loc[ii[j],'icustay_id']\n",
    "                \n",
    "for i in range(bacterio.shape[0]):\n",
    "    if(i%100000==0) :\n",
    "        print(i) \n",
    "    if(bacterio[i,2]==0): # if missing icustayid\n",
    "        subjectid = bacterio[i,0]\n",
    "        hadmid = bacterio[i,1]\n",
    "        jj = np.where((demog.loc[:,'subject_id']==subjectid) & (demog.loc[:,'hadm_id']==hadmid))[0]\n",
    "        if(len(jj)==1): \n",
    "            bacterio[i,2] = demog.loc[jj,'icustay_id']\n",
    "            \n",
    "\n",
    "# fill-in missing ICUSTAY IDs in ABx\n",
    "for i in range(abx.shape[0]):\n",
    "    if(np.isnan(abx[i,1])): \n",
    "        o = abx[i,2] # time of event\n",
    "        hadmid = abx[i,0]\n",
    "        ii = np.where(demog.loc[:,'hadm_id']==hadmid)[0] # row in table demographics\n",
    "        for j in range(len(ii)): \n",
    "            if (o >= demog.loc[ii[j],'intime']-48*3600) and (o <= demog.loc[ii[j],'outtime']+48*3600):\n",
    "                abx[i,1] = demog.loc[ii[j],'icustay_id']\n",
    "            elif len(ii) == 1 : # if we cant confirm from admission and discharge time but there is only 1 admission: it's the one!!\n",
    "                abx[i,1] = demog.loc[ii[j],'icustay_id']            \n",
    "            \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "########################################################################\n",
    "#    find presumed onset of infection according to sepsis3 guidelines\n",
    "########################################################################\n",
    "\n",
    "# METHOD:\n",
    "# I loop through all the ABx given, and as soon as there is a sample present\n",
    "# within the required time criteria I pick this flag and break the loop.\n",
    "\n",
    "\n",
    "onset = np.zeros((100000,3))\n",
    "\n",
    "# In Matlab, for icustayid=1:100000 means 1,2,3,...,100000 \n",
    "\n",
    "for icustayid in range(1,100001):\n",
    "    if(icustayid%10000==0):\n",
    "        print(icustayid)\n",
    "        \n",
    "    ab = abx[abx[:,1]==icustayid+200000,2] # start time of abx for this icustayid\n",
    "    bact = bacterio[bacterio[:,2]==icustayid+200000,3] # time of sample \n",
    "    subj_bact = bacterio[bacterio[:,2]==icustayid+200000,0] # subjectid\n",
    "    \n",
    "    if(ab.size!=0 and bact.size!=0):  # if we have data for both: proceed\n",
    "        D = np.zeros((ab.size,bact.size)) # pairwise distances btw antibio and cultures, in hours\n",
    "        for i,a in enumerate(ab):\n",
    "            for j,b in enumerate(bact): \n",
    "                D[i,j]=(math.sqrt((a-b)*(a-b))/3600) \n",
    "                \n",
    "        for i in range(D.shape[0]): # looping through all rows of AB given, from early to late\n",
    "            I = np.argmin(D[i]) # minimum distance in this row\n",
    "            M = D[i,I]\n",
    "           \n",
    "            ab1 = ab[i] # timestamp of this value in list of antibio\n",
    "            bact1 = bact[I] # timestamp in list of cultures\n",
    "                    \n",
    "            if M<=24 and ab1<=bact1 : # if ab was first and delay < 24h\n",
    "                onset[icustayid-1][0] = subj_bact[0] # subject_id\n",
    "                onset[icustayid-1][1]=icustayid # icustay_id\n",
    "                onset[icustayid-1][2]=ab1 # onset of infection = abx time \n",
    "                break \n",
    "            elif M<=72 and ab1>=bact1 : # elseif sample was first and elay < 72h \n",
    "                onset[icustayid-1][0] = subj_bact[0] # subject_id\n",
    "                onset[icustayid-1][1]=icustayid # icustay_id\n",
    "                onset[icustayid-1][2]=bact1 # onset of infection = sample time \n",
    "                break \n",
    "                \n",
    "                \n",
    "                \n",
    "            \n",
    "## Replacing item_ids with column numbers from reference tables\n",
    "\n",
    "# replace itemid in labs with column number\n",
    "# this will accelerate process later\n",
    "\n",
    "def replace_item_ids(reference,data):\n",
    "    temp = {} \n",
    "    for i in range(data.shape[0]): \n",
    "        key = data[i,2]\n",
    "        if(key not in temp):\n",
    "            temp[key] = np.argwhere(reference==key)[0][0]+1 # +1 because matlab index starts from 1 \n",
    "        data[i,2] = temp[key] \n",
    "\n",
    "replace_item_ids(Reflabs,labU)\n",
    "replace_item_ids(Refvitals,ce010)\n",
    "replace_item_ids(Refvitals,ce1020)\n",
    "replace_item_ids(Refvitals,ce2030)\n",
    "replace_item_ids(Refvitals,ce3040)\n",
    "replace_item_ids(Refvitals,ce4050)\n",
    "replace_item_ids(Refvitals,ce5060)\n",
    "replace_item_ids(Refvitals,ce6070)\n",
    "replace_item_ids(Refvitals,ce7080)\n",
    "replace_item_ids(Refvitals,ce8090)\n",
    "replace_item_ids(Refvitals,ce90100)"
   ]
  },
  {
   "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 -48 -> +24 because that's for sepsis3 cohort defintion!!\n",
    "# I need different time period for the MDP (-24 -> +48)\n",
    "\n",
    "reformat = np.full((2000000,68),np.nan)\n",
    "qstime=np.zeros((100000,4)) \n",
    "winb4=49   #lower limit for inclusion of data (48h before time flag)\n",
    "winaft=25 # upper limit (24h after)\n",
    "irow=0;  #recording row for summary table\n",
    "\n",
    "for icustayid in range(1,100001):\n",
    "    \n",
    "    if(icustayid%10000==0):\n",
    "        print(icustayid)\n",
    "    \n",
    "    qst=onset[icustayid-1,2] #flag for presumed infection\n",
    "    if qst>0 : # if we have a flag \n",
    "        ii = demog.loc[:,'icustay_id']==icustayid+200000\n",
    "        d1 = demog.loc[ii,['age','dischtime']].values[0]\n",
    "    \n",
    "        if d1[0]>6574: # if older than 18 years old \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",
    "    \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",
    "            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 #flag for presumed infection / this is time of sepsis if SOFA >=2 for this patient\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]=d1[1] # dischargetime \n",
    "\n",
    "                \n",
    "reformat = reformat[:irow,:]  #delete extra unused rows\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",
    "    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",
    "\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",
    "\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "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",
    "#ESTIMATE FiO2 /// with use of interface / device (cannula, mask, ventilator....)\n",
    "reformatsah=SAH(reformat,sample_and_hold) # do SAH first to handle this task\n",
    "\n",
    "\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"
   ]
  },
  {
   "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)\n"
   ]
  },
  {
   "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): # -52 until +28 = 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",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ########################################################################\n",
    "#    CONVERT TO TABLE AND DELETE VARIABLES WITH EXCESSIVE MISSINGNESS\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",
    "miss=(np.sum(np.isnan(reformat2),axis=0)/reformat2.shape[0])\n",
    "\n",
    "\n",
    "# if values have less than 70% missing values (over 30% of values present): I keep them\n",
    "reformat3t = reformat2t.iloc[:,np.hstack([np.full(11,True),(miss[11:74]<0.70),np.full(11,True)])]\n",
    " \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#########################################################################\n",
    "#             HANDLING OF MISSING VALUES  &  CREATE REFORMAT4T\n",
    "#########################################################################\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",
    "\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",
    "# 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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##########################################################################\n",
    "#        COMPUTE SOME DERIVED VARIABLES: P/F, Shock Index, SOFA, SIRS...\n",
    "##########################################################################\n",
    "\n",
    "\n",
    "# CORRECT GENDER\n",
    "reformat4t.loc[:,'gender']=reformat4t.loc[:,'gender']-1 \n",
    "\n",
    "# CORRECT AGE > 200 yo\n",
    "ii=reformat4t.loc[:,'age']>150*365.25\n",
    "reformat4t.loc[ii,'age']=91.4*365.25\n",
    "\n",
    "# FIX MECHVENT\n",
    "ii = np.isnan(reformat4t.loc[:,'mechvent'])\n",
    "reformat4t.loc[ii,'mechvent']=0\n",
    "ii = reformat4t.loc[:,'mechvent']>0\n",
    "reformat4t.loc[ii,'mechvent']=1\n",
    "\n",
    "# FIX Elixhauser missing values\n",
    "ii = np.isnan(reformat4t.loc[:,'elixhauser'])\n",
    "reformat4t.loc[ii,'elixhauser']=np.nanmedian(reformat4t.loc[:,'elixhauser'])  #use the median value / only a few missing data points \n",
    "\n",
    "# vasopressors / no NAN\n",
    "a=reformat4t.columns.get_loc('median_dose_vaso') \n",
    "ii=np.isnan(reformat4[:,a])\n",
    "reformat4t.loc[ii,'median_dose_vaso']=np.zeros((np.sum(ii)))\n",
    "a=reformat4t.columns.get_loc('max_dose_vaso') \n",
    "ii=np.isnan(reformat4[:,a])\n",
    "reformat4t.loc[ii,'max_dose_vaso']=np.zeros((np.sum(ii)))\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",
    "\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",
    "\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",
    "\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",
    "# 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",
    "#                            EXCLUSION OF SOME PATIENTS \n",
    "#########################################################################\n",
    "\n",
    "print(np.unique(reformat4t.loc[:,'icustayid']).size)  #count before\n",
    "\n",
    "# check for patients with extreme UO = outliers = to be deleted (>40 litres of UO per 4h!!)\n",
    "a = np.where(reformat4t.loc[:,'output_4hourly']>12000)[0]\n",
    "i = np.unique(reformat4t.loc[a,'icustayid'])\n",
    "i = np.where(np.isin(reformat4t.loc[:,'icustayid'],i))[0]\n",
    "reformat4t.drop(i,inplace=True)\n",
    "reformat4t.reset_index(inplace=True,drop=True)\n",
    "\n",
    "\n",
    "# some have bili = 999999\n",
    "a = np.where(reformat4t.loc[:,'Total_bili']>10000)[0]\n",
    "i = np.unique(reformat4t.loc[a,'icustayid'])\n",
    "i = np.where(np.isin(reformat4t.loc[:,'icustayid'],i))[0]\n",
    "reformat4t.drop(i,inplace=True)\n",
    "reformat4t.reset_index(inplace=True,drop=True)\n",
    "\n",
    "\n",
    "# check for patients with extreme INTAKE = outliers = to be deleted (>10 litres of intake per 4h!!)\n",
    "a=np.where(reformat4t.loc[:,'input_4hourly']>10000)[0]\n",
    "i = np.unique(reformat4t.loc[a,'icustayid'])  # 28 ids\n",
    "i = np.where(np.isin(reformat4t.loc[:,'icustayid'],i))[0]\n",
    "reformat4t.drop(i,inplace=True)\n",
    "reformat4t.reset_index(inplace=True,drop=True)\n",
    "\n",
    "# #### exclude early deaths from possible withdrawals ####\n",
    "# stats per patient\n",
    "\n",
    "q=reformat4t.loc[:,'bloc']==1\n",
    "\n",
    "# fence_posts=find(q(:,1)==1);\n",
    "num_of_trials=np.unique(reformat4t.loc[:,'icustayid']).size # size(fence_posts,1)\n",
    "a=reformat4t.loc[:,['icustayid','mortality_90d', 'max_dose_vaso', 'SOFA']]\n",
    "a.columns = ['id','mortality_90d','vaso','sofa']\n",
    "cnt = a.groupby('id').count()\n",
    "d=a.groupby('id').max() \n",
    "d['GroupCount'] = cnt['mortality_90d']\n",
    "d.reset_index(inplace=True)\n",
    "d.columns = ['id','max_mortality_90d','max_vaso','max_sofa','GroupCount']\n",
    "d=d[['id','GroupCount','max_mortality_90d','max_vaso','max_sofa']]\n",
    "\n",
    "# finds patients who match our criteria\n",
    "e=np.zeros([num_of_trials])\n",
    "\n",
    "for i in range(num_of_trials): \n",
    "    if d.loc[i,'max_mortality_90d']==1:\n",
    "        ii=(reformat4t['icustayid']==d.loc[i,'id']) & (reformat4t['bloc']==d.loc[i,'GroupCount'])  #last row for this patient\n",
    "        e[i]=sum((reformat4t.loc[ii,'max_dose_vaso']==0) & (d.loc[i,'max_vaso']>0.3) & (reformat4t.loc[ii,'SOFA']>=d.loc[i,'max_sofa']/2))>0\n",
    "\n",
    "\n",
    "r=d.loc[(e==1) & (d['GroupCount']<20),'id'] # ids to be removed\n",
    "ii=np.where(np.isin(reformat4t['icustayid'],r))[0]\n",
    "reformat4t.drop(ii,inplace=True)\n",
    "reformat4t.reset_index(inplace=True,drop=True)\n",
    "\n",
    "# exclude patients who died in ICU during data collection period\n",
    "ii=(reformat4t['bloc']==1) & (reformat4t['died_within_48h_of_out_time']==1) & (reformat4t['delay_end_of_record_and_discharge_or_death']<24)\n",
    "ii=np.where(np.isin(icustayidlist,reformat4t.loc[ii,'icustayid']))[0]\n",
    "reformat4t.drop(ii,inplace=True)\n",
    "reformat4t.reset_index(inplace=True,drop=True)\n",
    "\n",
    "print(np.unique(reformat4t.loc[:,'icustayid']).size)  #count after\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #######################################################################\n",
    "#                       CREATE SEPSIS COHORT\n",
    "# ########################################################################\n",
    "\n",
    "# create array with 1 row per icu admission\n",
    "# keep only patients with flagged sepsis (max sofa during time period of interest >= 2)\n",
    "# we assume baseline SOFA of zero (like other publications)\n",
    "\n",
    "sepsis=np.zeros([30000,5]) \n",
    "irow=0\n",
    "\n",
    "for icustayid in range(1,100001): \n",
    "    ii=np.where(np.isin(reformat4t['icustayid'],icustayid))[0]\n",
    "    \n",
    "    if(icustayid%10000==0): \n",
    "        print(icustayid)\n",
    "        \n",
    "    if(ii.size!=0):     \n",
    "        sofa=reformat4t.loc[ii,'SOFA'] \n",
    "        sirs=reformat4t.loc[ii,'SIRS'] \n",
    "        sepsis[irow,0]=icustayid+200000 \n",
    "        sepsis[irow,1]=reformat4t.loc[ii[0],'mortality_90d'] # 90-day mortality\n",
    "        sepsis[irow,2]=max(sofa)\n",
    "        sepsis[irow,3]=max(sirs)\n",
    "        sepsis[irow,4]=qstime[icustayid-1,0]   # time of onset of sepsis\n",
    "        irow=irow+1\n",
    "\n",
    "sepsis = sepsis[:irow,:]\n",
    "\n",
    "sepsis=pd.DataFrame(sepsis,columns = ['icustayid','morta_90d','max_sofa','max_sirs','sepsis_time']) \n",
    "\n",
    "# delete all non-sepsis\n",
    "sepsis.drop(np.where(sepsis['max_sofa']<2)[0],inplace=True) \n",
    "sepsis.reset_index(inplace=True,drop=True)\n",
    "\n",
    "# final count\n",
    "print(sepsis.shape[0])   \n",
    "\n",
    "# save cohort\n",
    "sepsis.to_csv('sepsis_mimiciii.csv',index=False,na_rep='NaN')   \n",
    "\n",
    "\n",
    "# Save to pickle for step 3 \n",
    "with open('step_3_start.pkl', 'wb') as file:\n",
    "    pickle.dump(sample_and_hold,file)\n",
    "    pickle.dump(demog,file)\n",
    "    pickle.dump(ce010,file)\n",
    "    pickle.dump(ce1020,file)\n",
    "    pickle.dump(ce2030,file)\n",
    "    pickle.dump(ce3040,file)\n",
    "    pickle.dump(ce4050,file)\n",
    "    pickle.dump(ce5060,file)\n",
    "    pickle.dump(ce6070,file)\n",
    "    pickle.dump(ce7080,file)\n",
    "    pickle.dump(ce8090,file)\n",
    "    pickle.dump(ce90100,file)\n",
    "    pickle.dump(labU,file)\n",
    "    pickle.dump(MV,file)\n",
    "    pickle.dump(inputpreadm,file)\n",
    "    pickle.dump(inputMV,file)\n",
    "    pickle.dump(inputCV,file)\n",
    "    pickle.dump(vasoMV,file)\n",
    "    pickle.dump(vasoCV,file)\n",
    "    pickle.dump(UOpreadm,file)\n",
    "    pickle.dump(UO,file) \n",
    "    pickle.dump(sepsis,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
}