In [1]:

    import csv
    import pandas as pd
    import os
    import glob
    from numpy import savetxt
    from scipy.spatial.distance import jaccard, pdist, squareform
    import numpy as np
    from sklearn.cluster import DBSCAN
    import matplotlib.pyplot as plt
    from sklearn import metrics
    from sklearn.preprocessing import StandardScaler
    from tqdm import tqdm

In [2]:

    os.getcwd()
    os.chdir(r'C:\Users\anson\Dropbox\Dissertation\Fracfocus\FracFocusCSV')

In [3]:

    ###find minpts and epsilon
    #the online appendix uses the austin chalk shale to illustrate how to find epsilon and minPts
    df = pd.read_csv('playDBSCAN_austin chalk.csv')

In [4]:

    df['date'] = pd.to_datetime(df['date'])
    df['Percentsum'] = df['PercentHFJob'].groupby(df['UploadKey']).transform('sum')
    df = df[df['PercentHFJob'] > 0]
    df_clean = df[df['Percentsum']!=0]
    df_clean = df_clean.sort_values(['date','UploadKey'])
    df_clean = df_clean.reset_index(drop=True)
    df_clean['gid'] = df_clean.groupby('UploadKey', sort=False).ngroup()
    last_gid = df_clean["gid"].iloc[-1]

In [5]:

    ###Using weighted jarcard distance to generate the chemical matrix

In [6]:

    def weighted_jaccard_distance(x, y):
        arr = np.array([x, y])
        return 1 - arr.min(axis=0).sum() / arr.max(axis=0).sum()

In [7]:

    def input_numpy(dbscan1000):
        p = dbscan1000.pivot_table(index='gid', columns='CAS_cleaned', values='PercentHFJob').fillna(0)
        sim_weighted = pd.DataFrame(data=squareform(pdist(p, weighted_jaccard_distance)), index=p.index, columns=p.index)
        sim_numpy = sim_weighted.to_numpy()
        return sim_numpy

In [8]:

    #minpts (for each shale play, use the average number of chemicals used per well)
    wellCAScount = df_clean.groupby('UploadKey').agg({'CAS_cleaned': ['count']})

In [9]:

    #epsilon
    sim_numpy = input_numpy(df_clean)
    sim_numpy_k = [row[25] for row in sim_numpy]
    sim_numpy_k = np.sort(sim_numpy_k)

In [10]:

    #reproduce online appendix figure b.2
    plt.figure(figsize=(15,10))
    plt.plot(sim_numpy_k)

    plt.xlabel('Number of Data Points sorted by distance',fontsize=18)
    plt.ylabel('Jaccard Distance',fontsize=18)
    plt.show()
    plt.savefig('knn.jpeg')

    <Figure size 640x480 with 0 Axes>

In [11]:

    ###iterate through each shale sample and use DBSCAN to obtain a list of novel (noise) wells
    #note that this method only identifies the first well using the novel chemical choice
    #i used the Austin Chalk shale sample above to reproduce online Appendix Figure B.2
    #it took approximately 17 hours to iterate through the Austin Chalk shale sample
    #to save time for iteration, I switched to a smaller shale sample below (the Granite Wash shale; epsilon=0.2, minPts=15)

In [12]:

    df = pd.read_csv('playDBSCAN_granite wash.csv')

In [13]:

    df['date'] = pd.to_datetime(df['date'])
    df['Percentsum'] = df['PercentHFJob'].groupby(df['UploadKey']).transform('sum')
    df = df[df['PercentHFJob'] > 0]
    df_clean = df[df['Percentsum']!=0]
    df_clean = df_clean.sort_values(['date','UploadKey'])
    df_clean = df_clean.reset_index(drop=True)
    df_clean['gid'] = df_clean.groupby('UploadKey', sort=False).ngroup()
    last_gid = df_clean["gid"].iloc[-1]

In [14]:

    #iteration starts here

In [15]:

    def dbscan_noisecheck (sim_numpy):
        db = DBSCAN(eps=0.2, min_samples=15, metric='precomputed')
        db.fit(sim_numpy)
        labels = db.labels_
        noise_index = list(np.nonzero(labels == -1)[0])

        if 1000 in noise_index:
            well_noise = True
        else:
            well_noise = False

        return well_noise

In [16]:

    uploadkley_noise = []
    for gid in tqdm(range(1000, last_gid+1)):
        dbscan1000 = df_clean[df_clean['gid'].isin(range(gid-1000,gid+1))]
        sim_numpy = input_numpy(dbscan1000)
        well_noise = dbscan_noisecheck(sim_numpy)

        if well_noise is True:
            noisekey = dbscan1000["UploadKey"].iloc[-1]
            uploadkley_noise.append(noisekey)
        else:
            continue

    100%|██████████| 418/418 [40:59<00:00,  5.88s/it]

In [17]:

    #save the list of novel wells

In [18]:

    df_uploadkeynoise = pd.DataFrame(uploadkley_noise, columns=['Uploadkey_noise'])
    df_uploadkeynoise.to_csv(r'playTS\noisekey\playDBSCAN_granite wash.csv', index=False, encoding='utf-8-sig')

In [19]:

    #save well clusters using the full shale sample (without iteration)

In [20]:

    db = DBSCAN(eps= 0.2, min_samples =15, metric='precomputed')
    sim_numpy = input_numpy(df_clean)
    db.fit(sim_numpy)
    labels = db.labels_
    set(labels)

Out[20]:

    {-1, 0, 1, 2, 3}

In [21]:

    cluster1_id = list(np.nonzero(labels == 0)[0])
    df_cluster1 = df_clean[df_clean['gid'].isin(cluster1_id)]
    df_cluster1.to_csv(r'playTS\cluster_1k\playDBSCAN_granite wash_cluster1.csv')

    cluster2_id = list(np.nonzero(labels == 1)[0])
    df_cluster2 = df_clean[df_clean['gid'].isin(cluster2_id)]
    df_cluster2.to_csv(r'playTS\cluster_1k\playDBSCAN_granite wash_cluster2.csv')

    cluster3_id = list(np.nonzero(labels == 2)[0])
    df_cluster3 = df_clean[df_clean['gid'].isin(cluster3_id)]
    df_cluster3.to_csv(r'playTS\cluster_1k\playDBSCAN_granite wash_cluster3.csv')

    cluster4_id = list(np.nonzero(labels == 3)[0])
    df_cluster4 = df_clean[df_clean['gid'].isin(cluster4_id)]
    df_cluster4.to_csv(r'playTS\cluster_1k\playDBSCAN_granite wash_cluster4.csv')

In [22]:

    ###repeat the above process for all shale samples