<jupyter_start><jupyter_text>Function Definitions<jupyter_code>def run(featurize, featurizeThese, columns):
    
    featurizedData = pd.DataFrame(columns=columns)

    i = 0
    print("Featurizing (sync)...\n")
    for row in featurizeThese.iterrows():
        i += 1
        print(f"{i} of {featurizeThese.shape[0]}")
        newRow = featurize(row)
        if newRow is not None:
            featurizedData = featurizedData.append(pd.DataFrame([newRow], columns=featurizedData.columns), ignore_index=True)
        else:
            print("None returned. Skipping!")
                
    return finalizeDataset(featurizedData)

def run_async(featurize, featurizeThese, columns, numProcs):
    
    featurizedData = pd.DataFrame(columns=columns)

    with Pool(numProcs) as p:

        i = 0
        num = len(featurizeThese)

        print("Featurizing (async)...\n")
        it = p.imap_unordered(featurize, featurizeThese.iterrows())
        while i < num:

            i += 1

            try:
                #newRow = it.next(timeout=600)
                newRow = it.next()
                print(f"{i} of {num} returned.")
                print("New row:", newRow, sep="\n")
                if newRow is not None:
                    featurizedData = featurizedData.append(pd.DataFrame([newRow], columns=featurizedData.columns), ignore_index=True)
                else:
                    print("None returned. Skipping!")
            except TimeoutError as te:
                print(f"Timed out. {te}")
                traceback.print_exc()
            except StopIteration:
                print("Pool iteration ended.")
                break
            except Exception as e:
                print(f"Other exception... {e}")
                traceback.print_exc()
            
    return finalizeDataset(featurizedData)
            
def finalizeDataset(ds):

    # Replace any inf values with nan
    return ds.replace([np.inf, -np.inf], np.nan)

# Open the file in both audata and h5py
def openFiles(fn):
    return aud.File.open(files_dir+fn, readonly=True), h5.File(files_dir+fn, 'r')
    
# Validate that data is at least minSamplesToBeValid, not a 0 timespan
def validateDataset(data):
    return False if data.shape[0] < minSamplesToBeValid or data.iloc[-1].time - data.iloc[0].time == 0 else True
   
# Convert a list/vector(?) response from R into a dict
def getDictFromR(o):
    return dict(zip(o.names, map(list,list(o))))

def isWaveform(s):
    return s.startswith('/data/waveforms')

# Return basic statistics for a pandas series
def getBasicStats(ds):
    
    features = []
    
    _mean = ds.mean()
    _median = ds.median()
    _std = ds.std()
    _min = ds.min()
    _max = ds.max()
    
    features.append(_mean) # mean
    features.append(_std) # sd
    features.append(_std/_mean if _mean != 0 else np.nan)# cv
    features.append((ds-_mean).median()) # mad
    features.append(hrrr.shape[0]) # n
    features.append(_min) # min
    features.append(_max) # max
    features.append(_median) # median
    features.append(_max-_min) # range
    features.append((_max-_min)/_median if _median != 0 else np.nan)# range_ratio
    
    return features

# Get best numeric series from available options
# f=file, so=series options, left=left timestamp
def getBestSeries(f, h5f, so, left):
    
    printout and print(f"getBestSeries(f, h5f, {so}, {left}) started.")
    
    finalSeries = None
    finalSeriesSampleCount = -1
    finalSeriesTREF = None
    finalSeriesDF = None
    finalSeriesData = None

    # Iterate through the series options
    for s in so:

        # If the series is not available in the file, skip
        if s not in f:
            continue
            
        tref = f[s].time_reference.timestamp()

        if isWaveform(s):
            
            df = h5f[s]
            data = pd.DataFrame(df[(df['time'] > left-tref) & (df['time'] < left+60*3-tref)])
        
        else:
            
            df = pd.DataFrame(f[s].get(raw=True)).drop_duplicates(subset=['time'], ignore_index=True)

            # Downsample if requested. Note we do not downsample NBP because it is already a low-rate measurement.
            if downsample_to_ppinnc1 and sk in ['hr', 'rr', 'spo2']:
                df = df.iloc[::20].reset_index(drop=True)

            data = df[ (df.time > left-tref) & (df.time < left+60*3-tref) ]

        # If data validates and has more samples than current final series, it's our new final series
        if validateDataset(data) and data.shape[0] > finalSeriesSampleCount:

            finalSeries = s
            finalSeriesSampleCount = data.shape[0]
            finalSeriesTREF = tref
            finalSeriesDF = df
            finalSeriesData = data
            
    printout and print(f"getBestSeries(f, h5f, {so}, {left}) returning.")

    return finalSeries, finalSeriesTREF, finalSeriesDF, finalSeriesData

def featurize(row):

    # Grab the row content
    row = row[1]
    
    printout and print(f"File: {row['filename']}\nAlert Series: {row['series']}\nLabel: {'REAL' if row['real_vs_art']==1 else 'ARTIFACT'}")

    # Open the file in both audata and h5py
    (f, h5f) = openFiles(row['filename'])
    
    # This will hold this alert's computed features
    features = [row['filename'], row['series'], row['left'], row['right'], row['real_vs_art']]
    
    # Iterate through the signal types (each signal type e.g. hr, spo2 will get one set of features)
    for sk, so in seriesOptions.items():
        
        # TODO(gus): If this sk is for the label target series category, use the label target series.
        
        (series, tref, df, data) = getBestSeries(f, h5f, so, row['left'])
        
        # If we have no available series for this signal type, add empty values for the features and move on.
        if series is None:
            print(f"No data available for {sk}!")
            
            # TODO(gus): temp
            if sk == 'pleth' or sk == 'ecgii':
                # this is not temp, just pull it out of the if
                features = features + ( [np.nan]*(len(featuresList) + (3 if sk == 'pleth' else 0)) )
            
            # If we have no data for our CRI target series (e.g. if we've invalidated it), skip this CRI alert entirely
            if row['series'].split(':')[0] in seriesOptions[sk]:
                print(f"Not enough / invalid data for for CRI target series. Discluding this CRI alert!")
                return None
            
            continue
        
        # TODO(gus): temp
        if sk != 'pleth' and sk != 'ecgii':
            print(f"Temporarily skipping {sk}")
            continue
        
        printout and print(f"Featurizing series {series}, timespan: {round((data.iloc[-1].time - data.iloc[0].time) / 60, 2)} minutes, sample count: {data.shape[0]}")
        
        _mean = data.value.mean()
        _median = data.value.median()
        _std = data.value.std()
        _min = data.value.min()
        _max = data.value.max()
        features.append(_mean) # mean
        features.append(_std) # sd
        features.append(_std/_mean if _mean != 0 else np.nan) # cv
        features.append((data.value - _mean).median()) # mad
        features.append(data.shape[0]) # n
        features.append(_min) # min
        features.append(_max) # max
        features.append(_median) # median
        features.append(_max-_min) # range
        features.append((_max-_min)/_median if _median != 0 else np.nan) # range_ratio
        
        printout and print("Retrieving trail.")
        trail = df[ (df['time'] < row['left']-tref) & (df['time'] > row['left']-tref-15*60) ]
        if isWaveform(series):
            trail = pd.DataFrame(trail)
        printout and print("Retrieving trail complete.")

        features.append( (data.shape[0])/(row['right']-row['left']) ) # data_den. note: use cri timespan, not sample timespan

        trail_count = trail.shape[0]
        if not validateDataset(trail):
            features = features + ( [np.nan]*2 ) # data_den_trail, data_den_trail2
        else:
            features.append(trail_count / (trail.time.iloc[-1]-trail.time.iloc[0])) # data_den_trail
            features.append(trail_count / (15*60)) # data_den_trail2

        # TODO(gus): Ignoring this for waveforms for now
        if isWaveform(series):
            features.append(np.nan) # delta_t
        else:
            features.append(df.loc[data.index[0], 'time'] - df.loc[data.index[0]-1, 'time'] if data.index[0] > 0 else np.nan) # delta_t
        features.append(data.time.diff().max()) # max_gap

        # skipping for now... spec_ratio
        # skipping for now... max_spec
#         with localconverter(robjects.default_converter + pandas2ri.converter):
#             ret = getDictFromR(robjects.r['get_spectral_lomb'](data['time'], data['value']))
#         print('spec:', ret)

        printout and print("Retrieving dl4.")
        dl4 = df[(df['time'] >= row['left']-tref-4*60) & (df['time'] < row['left']-tref)]
        if isWaveform(series):
            dl4 = pd.DataFrame(dl4)
        printout and print("Retrieving dl4 complete.")

        if not validateDataset(dl4):

            features = features + ( [np.nan]*10 ) # delta_mean, delta_std, MWW-test statistic & p-value, KS-test statistic & p-value, t-test statistic & p-value, f-test statistic & p-value

        else:

            features.append(abs(dl4.value.mean()-_mean)) # delta_mean
            features.append(abs(dl4.value.std()-_std)) # delta_std

            #print(data['value'], dl4['value'])
            try:
                (mww_stat, mww_pval) = scipy.stats.mannwhitneyu(data['value'], dl4['value'])
                features.append(mww_stat)# MWW-test statistic
                features.append(mww_pval)# MWW-test p-value
            except:
                features = features + ( [np.nan]*2 )

            (ks_stat, ks_pval) = scipy.stats.ks_2samp(data['value'], dl4['value'])
            features.append(ks_stat) # KS-test statistic
            features.append(ks_pval) # KS-test p-value

            (tt_stat, tt_pval) = scipy.stats.ttest_ind(data['value'], dl4['value'])
            features.append(tt_stat) # t-test statistic
            features.append(tt_pval) # t-test p-value

            try:
                with localconverter(robjects.default_converter + pandas2ri.converter):
                    ret = getDictFromR(robjects.r['var.test'](data['value'], dl4['value']))
                features.append(ret['statistic'][0]) # f-test statistic
                features.append(ret['p.value'][0]) # f-test p-value
            except Exception:
                traceback.print_exc()
                print(row, data['value'], dl4['value'], sep="\n")
                features.append(np.nan)
                features.append(np.nan)

        # this section calculates f-test in python, but results often differ from r
#         x = data['value']
#         y = dl4['value']
#         f = x.var()/y.var()
#         p = 2*scipy.stats.f.cdf(f, x.shape[0]-1, y.shape[0]-1)
#         print(f, p)
#         print('p-python/p-r', vtr['p.value'][0]/p)

        slrm = LinearRegression()
        slrm.fit(data[['time']], data['value'])
        features.append(slrm.coef_[0]) # slope

        with localconverter(robjects.default_converter + pandas2ri.converter):
            ret = robjects.r['get_robust_slope'](data['time'], data['value'])
        features.append(ret[0]) # rslope

        with localconverter(robjects.default_converter + pandas2ri.converter):
            ret = robjects.r['get_slope_break'](data['time'], data['value'])
        features.append(ret[0][0]) # slope_before
        features.append(ret[1][0]) # slope_after
        features.append(ret[2][0]) # num_breakpoint

        data_dropduptimes = data.drop_duplicates(subset=['time'], keep='last')
        dydx = (data_dropduptimes.value.diff() / data_dropduptimes.time.diff())[1:]
        features.append(dydx.max()) # diff1_max
        features.append(dydx.min()) # diff1_min

        # TODO(gus): Ignoring this for waveforms for now
        if isWaveform(series):
            features = features + ([np.nan]*2)
        else:
            with localconverter(robjects.default_converter + pandas2ri.converter):
                ret = getDictFromR(robjects.r['comp_grads'](data['time'], data['value']))
            features.append(ret['max_grad'][0]) # max_grad
            features.append(ret['min_grad'][0]) # min_grad

        with localconverter(robjects.default_converter + pandas2ri.converter):
            ret = getDictFromR(robjects.r['get_rsq'](data['time'], data['value']))
        features.append(ret['rsq'][0]) # quad_rsq
        features.append(ret['coef1'][0]) # quad_coef1
        features.append(ret['coef2'][0]) # quad_coef2
        features.append(ret['var.res'][0]) # quad_resvar

        with localconverter(robjects.default_converter + pandas2ri.converter):
            ret = getDictFromR(robjects.r['get.osci.index'](data['time'], data['value']))
        features.append(ret['up'][0]) # osi_up
        features.append(ret['down'][0]) # osi_down
        features.append(ret['ratio'][0]) # osi_ratio
        
        if sk == 'pleth':
            
            data['time'] = pd.to_datetime(data['time']+tref, unit='s')
            data.set_index('time', inplace=True)
            features.append(data.resample('3s').agg(lambda x: ((x-x.mean()/x.std())**3).sum()/x.shape[0]).mean()[0]) # skew
            features.append(data.resample('3s').agg(lambda x: ((x-x.mean()/x.std())**4).sum()/x.shape[0]).mean()[0]) # kurtosis
            features.append(data.resample('3s').agg(lambda x: (x**2 * np.log(x**2)).sum()*-1).mean()[0]) # entropy

    print("File done.")
    printout and print()
    
    return features

def featurize_targeted(row):
    
    # Grab the row content
    row = row[1]
    
    # Open the file in both audata and h5py
    (f, h5f) = openFiles(row['filename'])
    
    # This will hold this alert's computed features
    features = [row['filename'], row['series'], row['left'], row['right'], row['label']]
    
    sk = 'ecgii'
    
    (series, tref, df, data) = getBestSeries(f, h5f, seriesOptions[sk], row['left'])
    
    # If we have no available series for this signal type, add empty values for the features and move on.
    if series is None:
        
        print(f"No data available for {sk}!")
        features = features + [np.nan]
        
    else:
        
        try:
        
            freq = int(round(1/(data.loc[1, 'time']-data.loc[0, 'time'])))

            time_r_peaks, amplitude_r_peaks = bsnb.detect_r_peaks(data['value'], freq, time_units=True, plot_result=False)

            vpp_signal_ecg = np.ptp(data['value'])

            # For this task, we will follow the same procedure as shown before, but store the values in a list, so that we can then calculate the mean value.
            vpp_noise_ecg = []
            for t in time_r_peaks:
                start = int((t + 0.5) * freq) # 0.5 - time between a peak and a flat 
                end = int((t + 0.65)* freq) # 0.65 time between a peak and the end of the flat
                interval = data['value'][start:end]
                if interval.shape[0] == 0:
                    print("Skipping len=0 noise segment.")
                    continue
                vpp = np.ptp(interval)
                vpp_noise_ecg.append(vpp)

            vpp_noise_ecg = np.mean(vpp_noise_ecg)

            features.append(vpp_signal_ecg/vpp_noise_ecg)
            
        except Exception as e:
            
            print(f"Exception occurred while computing SNR: {e}")
            features.append(np.nan)
        
    (seriesHR, trefHR, dfHR, dataHR) = getBestSeries(f, h5f, seriesOptions['hr'], row['left'])
    (seriesRR, trefRR, dfRR, dataRR) = getBestSeries(f, h5f, seriesOptions['rr'], row['left'])
    if seriesHR is None or seriesRR is None:
        print(f"No data available for hr+rr!")
        features = features + [np.nan]*10
        
    else:
        
        dataHR.set_index('time', inplace=True)
        dataRR.set_index('time', inplace=True)
        hrrr = pd.concat([dataHR, dataRR], axis=1).dropna()
        hrrr.columns = ['hr', 'rr']
        features = features + getBasicStats(hrrr['hr'] / hrrr['rr'])

    return features<jupyter_output><empty_output><jupyter_text>Action<jupyter_code>downsample_to_ppinnc1 = False
printout = False
minSamplesToBeValid = 5

# Standard list of features to compute for each series
featuresList = [
    'mean', 'sd', 'cv', 'mad', 'n', 'min', 'max', 'median', 'range', 'range_ratio',
    'data_den', 'data_den_trail', 'data_den_trail2', 'delta_t', 'max_gap',
    # 'spec_ratio', 'max_spec',
    'delta_mean', 'delta_sd', 'MWW_stat', 'MWW_pvalue', 'KS_stat', 'KS_pvalue', 't_stat', 't_pvalue', 'F_stat', 'F_pvalue',
    'slope', 'rslope', 'slope_before', 'slope_after', 'num_breakpoint', 'diff1_max', 'diff1_min', 'max_grad', 'min_grad', 'quad_rsq', 'quad_coef1', 'quad_coef2', 'quad_resvar', 'osi_up', 'osi_down', 'osi_ratio',
]

# List of series to compute features for, by signal type
seriesOptions = {
    'hr':    ['/data/numerics/HR.HR'],
    'rr':    ['/data/numerics/RR.RR'],
    'spo2':  ['/data/numerics/SpO₂.SpO₂', '/data/numerics/SpO₂T.SpO₂T'],
    'bpd':   ['/data/numerics/NBP.NBPd', '/data/numerics/NBP.NBP-D'],
    'bps':   ['/data/numerics/NBP.NBPs', '/data/numerics/NBP.NBP-S'],
}
if not downsample_to_ppinnc1:
    seriesOptions['abpd'] = ['/data/numerics/AR1-D', '/data/numerics/ART.Diastolic']
    seriesOptions['abps'] = ['/data/numerics/AR1-S', '/data/numerics/ART.Systolic']
    seriesOptions['ecgii'] = ['/data/waveforms/II']
    seriesOptions['pleth'] = ['/data/waveforms/Pleth', '/data/waveforms/PlethT']

# The columns end up being filename, series, left, right, label, and then one set of featuresList column names for each series.
columns = \
    ['filename', 'series', 'left', 'right', 'label']+ \
    [ser +'_' + feat for ser in ['hr', 'rr', 'spo2', 'bpd', 'bps'] for feat in featuresList]+ \
    ([ser +'_' + feat for ser in ['abpd', 'abps'] for feat in featuresList] if not downsample_to_ppinnc1 else [])+ \
    ['ecgii_' + feat for feat in featuresList]+ \
    ['pleth_' + feat for feat in featuresList+['skew', 'kurtosis', 'entropy']]
    
# Label set to featurize
# k=0; featurizeThese = labels.iloc[k:k+1]
# featurizeThese = labels[20:45]
# featurizeThese = labels.drop(labels.index[97])[95:]
featurizeThese = labels

# run(featurize, featurizeThese, columns)


# Label set to featurize
prevFeatureSet = pd.read_csv('data_full_pleth2.csv')[317:318]
columns = \
    ['filename', 'series', 'left', 'right', 'label', 'ecgii_snr']+ \
    ['hrrr_ratio_'+feat for feat in ['mean', 'sd', 'cv', 'mad', 'n', 'min', 'max', 'median', 'range', 'range_ratio']]
# featurizedData = run(featurize_targeted, prevFeatureSet, columns)
featurizedData = run_async(featurize_targeted, prevFeatureSet, columns, 18)




print("Run complete.")

with pd.option_context("display.max_columns", 1000):
    print("Featurized Data:")
    display(featurizedData)
labels
prevFeatureSet<jupyter_output><empty_output><jupyter_text>Export<jupyter_code>featurizedDataMissingOneSafe = featurizedDataMissingOne.copy()
# ### Export to model-specific datasets

# suf = '_pleth'

# ds = featurizedData_PLETH

# ds.to_csv(f'data_full{suf}.csv', index=False)
# ds.drop(columns=['series', 'left', 'right']).to_csv(f'data_withfn{suf}.csv', index=False)

# ds[ds['series'] == '/data/numerics/HR.HR:value'].drop(columns=['series', 'left', 'right']).to_csv(f'data_hr_withfn{suf}.csv', index=False)
# ds[ds['series'] == '/data/numerics/RR.RR:value'].drop(columns=['series', 'left', 'right']).to_csv(f'data_rr_withfn{suf}.csv', index=False)
# ds[ds['series'].isin(['/data/numerics/SpO₂.SpO₂:value', '/data/numerics/SpO₂T.SpO₂T:value'])].drop(columns=['series', 'left', 'right']).to_csv(f'data_spo2_withfn{suf}.csv', index=False)
featurizedData.sort_values(by=['filename', 'series', 'left'])
prevFeatureSet.sort_values(by=['filename', 'series', 'left'])
a=featurizedData.sort_values(by=['filename', 'series', 'left']).reset_index()
b=prevFeatureSet.sort_values(by=['filename', 'series', 'left']).reset_index()
a[~a.isin(b).all(1)]
a=a[['filename', 'series', 'left']]
b=b[['filename', 'series', 'left']]
a.values.difference(b.values)
ds1 = set([tuple(line) for line in a.values])
ds2 = set([tuple(line) for line in b.values])
ds2.difference(ds1)
pd.concat([featurizedDataMissingOne, featurizedData], ignore_index=True).copy()[['filename', 'series', 'left', 'right']].drop_duplicates()
for f in list(nd.columns):
    print(f)
### Combine datasets (for targeted featurization)

# new data
#nd = featurizedData.copy()
nd = pd.concat([featurizedDataMissingOne, featurizedData], ignore_index=True).copy()
nd.set_index(['filename', 'series', 'left', 'right'], inplace=True)
nd.sort_index(inplace=True)
nd.drop(columns=['label'], inplace=True)

# old data
#od = pd.read_csv('data_full_pleth15s.csv')
#od = prevFeatureSet.copy()
od = pd.read_csv('data_full_pleth2.csv')
od.set_index(['filename', 'series', 'left', 'right'], inplace=True)
od.sort_index(inplace=True)
#od.drop(columns=['pleth_skew', 'pleth_kurtosis', 'pleth_entropy'], inplace=True)

combined = pd.concat([od, nd], axis=1).reset_index()
with pd.option_context("display.max_columns", 1000):
    display(combined)
with pd.option_context("display.max_columns", 1000):
    #display(combined.drop(columns=[pre+'_'+feat for pre in ['pleth', 'ecgii'] for feat in ['mean', 'sd', 'cv', 'mad', 'n', 'min', 'max', 'median', 'range', 'range_ratio']]))
    display(combined.drop(columns=[pre+'_'+feat for pre in ['bps', 'bpd'] for feat in featuresList]+[pre+'_'+feat for pre in ['pleth', 'ecgii'] for feat in ['mean', 'sd', 'cv', 'mad', 'n', 'min', 'max', 'median', 'range', 'range_ratio']]))
### Export to model-specific datasets

suf = '_try2'

ds = combined.drop(columns=[pre+'_'+feat for pre in ['bps', 'bpd'] for feat in featuresList]+[pre+'_'+feat for pre in ['pleth', 'ecgii'] for feat in ['mean', 'sd', 'cv', 'mad', 'n', 'min', 'max', 'median', 'range', 'range_ratio']])

ds.to_csv(f'data_full{suf}.csv', index=False)
ds.drop(columns=['series', 'left', 'right']).to_csv(f'data_withfn{suf}.csv', index=False)

ds[ds['series'] == '/data/numerics/HR.HR:value'].drop(columns=['series', 'left', 'right']).to_csv(f'data_hr_withfn{suf}.csv', index=False)
ds[ds['series'] == '/data/numerics/RR.RR:value'].drop(columns=['series', 'left', 'right']).to_csv(f'data_rr_withfn{suf}.csv', index=False)
ds[ds['series'].isin(['/data/numerics/SpO₂.SpO₂:value', '/data/numerics/SpO₂T.SpO₂T:value'])].drop(columns=['series', 'left', 'right']).to_csv(f'data_spo2_withfn{suf}.csv', index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15728607862105357, 'fraction_numerical': 0.010069474980212823, 'mean_word_length': 2.922559503276992, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '421275', 'license': 'MIT', 'n_tokens_mistral': 7799, 'n_tokens_neox': 7582, 'n_words': 1640, 'path': 'docs/example_notebooks/PPINNC – Featurization.ipynb', 'repo_name': 'autonlab/auviewer', 'size': '22742'}