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

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


###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')

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]


###Using weighted jarcard distance to generate the chemical matrix
def weighted_jaccard_distance(x, y):
    arr = np.array([x, y])
    return 1 - arr.min(axis=0).sum() / arr.max(axis=0).sum()

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

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

#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)

#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')


###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)

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

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]

#iteration starts here
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

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 

#save the list of novel wells
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')

#save well clusters using the full shale sample (without iteration)
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)

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')


###repeat the above process for all shale samples

