<jupyter_start><jupyter_text>To start, we import all of the packages that we will need in this notebook.<jupyter_code>import glob
import pandas as pd
import numpy as np
from scipy.integrate import quad
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import seaborn as sns<jupyter_output><empty_output><jupyter_text>Stage 1: Preprocess The data  Adding the features we will need to optimize to a DataFrame and then exporting to a csv Here we read in the well production.csv file<jupyter_code>well_productions = pd.read_csv("well productions/well production.csv")<jupyter_output><empty_output><jupyter_text>Here we read in the csv file for each well<jupyter_code>import csv
datas = []
for file in glob.glob("well productions/*"):
    if "well production.csv" not in file:
        frame = pd.read_csv(file)
        # Strip off the extra things on the end
        frame["Name"] = file[17:-4]
        datas.append(frame)<jupyter_output><empty_output><jupyter_text>The following function calculates the well length for each well and adds the well length of each well to the DataFrame<jupyter_code>def well_length(dataframe: pd.DataFrame):
    dataframe["well length"] = dataframe["easting"].iloc[-1] - dataframe["easting"][0]<jupyter_output><empty_output><jupyter_text>The following function calculates the number of frac stages for each well and adds the calculated value to the DataFrame<jupyter_code>def frac_stages(dataframe : pd.DataFrame):
    dataframe["frac stages"] = dataframe[dataframe["proppant weight (lbs)"].isna() == False].shape[0]<jupyter_output><empty_output><jupyter_text>The following function sets proppant per stage as "ppf." It will default to the maximum weight.<jupyter_code>"""
Sets a proppant per stage as "ppf" defaults to the maximum weight
method can be ["min, "max", "avg"]
"""
def proppant_per_stage(dataframe: pd.DataFrame, method="avg"):
    if method == "min":
        val = min(dataframe["proppant weight (lbs)"])
    elif method == "max":
        val = max(dataframe["proppant weight (lbs)"])
    elif method == "avg":
        val = dataframe["proppant weight (lbs)"].describe()["mean"]
    else:
        raise Exception(f"No method found for {method}")
    dataframe["ppf"] = val<jupyter_output><empty_output><jupyter_text>The following function sets pump rate as "pr." It will default to the maximum rate.<jupyter_code>"""
Sets a pump rate as "pr" defaults to the maximum weight
method can be ["min, "max", "avg"]
"""
def pump_rate(dataframe: pd.DataFrame, method="avg"):
    if method == "min":
        val = min(dataframe["pump rate (cubic feet/min)"])
    elif method == "max":
        val = max(dataframe["pump rate (cubic feet/min)"])
    elif method == "avg":
        val = dataframe["pump rate (cubic feet/min)"].describe()["mean"]
    else:
        raise Exception(f"No method found for {method}")
    dataframe["pr"] = val<jupyter_output><empty_output><jupyter_text>The below function finds how many years the well will economically produce oil, under the assumption that up to 93 barrels of oil is economical. To arrive at this value we **NEED TO RATIONALIZE THIS VALUE**<jupyter_code>def life_of_res(qi, D):
    lifetime = 1/D*np.log(qi/93)
    if lifetime < 0:
        return 0
    else: 
        return lifetime<jupyter_output><empty_output><jupyter_text>exponential loss curve$ q(t) = q_ie^{-Dt} $Where $D$ is the decline rate and $q_i$ is the initial rate of productionExponential decline curve equation    Arguments:        t:  Number of months the well has been producing oil         qi: Float. Initial production rate when well first came online.        di: Float. Nominal decline rate (constant)            Output:         Returns q, or the expected production rate at time t. Float.<jupyter_code>def exponential_loss(t, qi, D): 
    return qi*np.exp(-D*t)<jupyter_output><empty_output><jupyter_text>$ q(t) = \frac{q_i}{(1+bDt)^{1.0/b}} $  Hyperbolic decline curve equation     Arguments:        t:  Number of months the well has been producing oil        qi: Float. Initial production rate when well first came online.        b:  Float. Hyperbolic decline constant        di: Float. Nominal decline rate at time t=0    Output:         Returns q, or the expected production rate at time t. Float.<jupyter_code>def hyperbolic_loss(t, qi, b, di):
    return qi/((1.0+b*di*t)**(1.0/b))<jupyter_output><empty_output><jupyter_text>This function allows you to look at the first X months of production, and selects the highest production month as max initial production. It returns the max initial production in the first X months along with a series that contains the values of oil production for the first 12 months.         Arguments:        number_first_months: int. Number of months from the point the well comes online        to compare to get the max initial production rate qi (this looks at multiple months        in case there is a production ramp-up)                well_name: String. name of the well where we're attempting to get        the max volume from.<jupyter_code>def get_max_initial_production(number_first_months, well_name):
    
    row = well_productions.loc[well_productions["well name"] == well_name]
    val=0      
    row = row.filter(regex='oil')
    row=row.T.squeeze()
    
    for i in range(number_first_months):
        val = max(row[i], val)

    return val, row<jupyter_output><empty_output><jupyter_text>The below two functions find the defined integral of the fitted exponential and hyperbolic equations between 0 and the calculated life of the reservoir. Our team found that was easier to use an integral then implement the equation fround in the "Hinge Basin" notebook. $ N_p = \int_0^{L_t} q(t)dt $ Where $N_p$ is total production and $L_t$ is the lifetime of the well<jupyter_code>def get_cumulative_exponential(qi, D):
    return quad(exponential_loss, 0, life_of_res(qi, D), args=(qi,D))

def get_cumulative_hyperbolic(qi, b, di):
    return quad(hyperbolic_loss, 0, life_of_res(qi, di), args=(qi, b, di))<jupyter_output><empty_output><jupyter_text>$N_{p} = (\frac{q_{i}}{D})ln(\frac{q_{i}}{q})$<jupyter_code>def get_cumulative_other(qi, D, qt):
    t = life_of_res(qi, D)
    q = qt(t, qi, D) # for the exponential equation
    return (qi/D) * np.log(qi/q)<jupyter_output><empty_output><jupyter_text>The following function will add the expected lifetime of the well when using the exponential_loss function and the hyperbolic function then add the respective values to the dataframe. Following that, it will find the cumulative production with the least amount of error out of the two functions, and add that value to the DataFrame.<jupyter_code>time_series=pd.Series(list(range(12)))
def cum_production(dataframe: pd.DataFrame, find_min=False):
    #find the name of the current well being examined
    name = dataframe["Name"][0]
    #Set qi equal to the maximum production in the first 5 months of this well's lifespan, and row equal to the series of the 
    #first 12 months' production values
    qi, row = get_max_initial_production(5, name)
    
    #Use the scipy curve_fit function to get the best possible exponential and hyperbolic curves
    popt_exp, pcov_exp = curve_fit(exponential_loss, time_series, row, bounds=(0, [qi,20]))
    popt_hyp, pcov_hyp=curve_fit(hyperbolic_loss, time_series, row,bounds=(0, [qi,2,20]))
    
    #Get the defined integral's value and error for both the exponential and hyperbolic curves
    cp_exp = get_cumulative_exponential(*popt_exp)
    cp_hyp = get_cumulative_hyperbolic(*popt_hyp)

    #Set cumulative production to whichever defined integral has the least error
    if cp_exp[1] < cp_hyp[1]: 
        cum_production = cp_exp[0]
    else: 
        cum_production = cp_hyp[0]

    #Add the respected calculated lifetimes for the exponential and hyperbolic curves to the dataframe then add the 
    #calculated cumulative production. 
    dataframe["lifetime_exp"] = life_of_res(*popt_exp)
    dataframe["lifetime_hyp"] = life_of_res(popt_hyp[0], popt_hyp[2])
    dataframe["cum_production"] = cum_production
    
    
    dataframe["cum_production_other"] = get_cumulative_other(popt_exp[0], popt_exp[1], exponential_loss)<jupyter_output><empty_output><jupyter_text>Here we call all of the functions we wrote that add the values we will need to optimize to the DataFrame<jupyter_code>%%capture
list(map(cum_production, datas))
list(map(well_length, datas))
list(map(frac_stages, datas))
list(map(proppant_per_stage, datas))
list(map(pump_rate, datas))<jupyter_output><empty_output><jupyter_text>Here we drop the duplicate wells from the dataframe<jupyter_code>%%capture
big_df = pd.concat(datas)
## Removing duplicates makes graphing well locations significantly harder. Remove duplicates when regressing
# big_df.drop_duplicates(subset=['Name'], inplace=True)<jupyter_output><empty_output><jupyter_text>Here we export our DataFrame to a new csv file so we no longer need to run preprocessing cells above<jupyter_code>big_df.to_csv("bigPoppa.csv", index=False)<jupyter_output><empty_output><jupyter_text>Stage2: Creating a model to predict cumulative output  Creating a model that uses the features we added to the dataframe to predict cumulative output **We should finally be able to get to the fun stuff now 😁**<jupyter_code>big_df = pd.read_csv("bigPoppa.csv")
data = big_df.drop_duplicates(subset=['Name'])
data["lifetime_exp"].describe() # should be measured in years<jupyter_output><empty_output><jupyter_text>In order to apply any sort of regression to our data, we must first clean our data and remove all non-numerical values that will not add any specific input into our final predictions. For this dataset, we see that the Names column along with any rows that have NaN/Infinity/Empty data should not be included in our model. Therefore, we will remove these anomalies from our data.<jupyter_code>def clean_dataset(df):
    df = df.dropna()
    indices_to_keep = ~df.isin([np.nan, np.inf, -np.inf]).any(1)
    return df[indices_to_keep].astype(np.float64)
from scipy import stats<jupyter_output><empty_output><jupyter_text>Remove outliers on the right...right skewed graph<jupyter_code>clean_data = clean_dataset(data.drop(["Name"], axis=1))
clean_data = clean_data[np.abs(stats.zscore(clean_data) < 2.5).all(axis=1)]<jupyter_output><empty_output><jupyter_text>4027 After cleaning our data, we will use all but the cumulative production data in our dataframe to predict our cumulative production for each well.  Find the Important Features<jupyter_code>def feature_importances(model, features_df, thresh=0.1):
    print("  Feature".ljust(36) + "  Score")
    print("-"*60)
    for i, feature in enumerate(model.feature_importances_):
        if abs(feature) >= thresh:
            print("| ", end="")
            print(f"{features_df.columns[i]}".ljust(36), end="|")
            print(f"{feature}".ljust(20), end="")
            print(" |")
            print("-"*60)
def graph_acuracy(x, y, data, model, feature):
    plt.scatter(x[feature], rf.predict(x), label="Prediction")
    plt.scatter(x[feature], y, label="Acutal")
    plt.xlabel(f"{feature} value")
    plt.ylabel("Cumulative production")
    plt.legend()
def graph_accuracy_importances(x, y, model, data, importances, thresh=0.1):
    for i, feature in enumerate(importances):
        if abs(feature) >= thresh:
            feat_name = x.columns[i]
            graph_acuracy(x, y, data, model, feat_name)
            plt.title(f"Acuracy of {feat_name} - Importance {feature}")
            plt.show()
from yellowbrick.regressor import ResidualsPlot

def show_residusal(model, train_tup, test_tup):
    resPlot = ResidualsPlot(model)
    resPlot.fit(*train_tup)
    resPlot.score(*test_tup)
    resPlot.show()
!pip install yellowbrick
data.head()
drop_columns = ['cum_production', 'cum_production_other', 'lifetime_hyp', 'lifetime_exp','frac stages', 'northing']
# Changed this to just look at the values we are supposed to 
x = clean_data.drop(drop_columns, axis = 1)
#x = clean_data[["ppf", "pr", "well length", "frac stages"]]
y = np.log(clean_data['cum_production']) # predict on the ln(cum) hopefully will linearlize it
from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x,y, test_size=0.25, random_state=86)
train, test = (x_train, y_train), (x_test, y_test)
from sklearn.ensemble import RandomForestRegressor
rf = RandomForestRegressor(n_jobs=11,  random_state=86)
show_residusal(rf, train, test)
feature_importances(rf, x, 0.1)<jupyter_output>Feature                             Score
------------------------------------------------------------
| easting                             |0.12935111374237     |
------------------------------------------------------------
| permeability                        |0.12003789461556819  |
------------------------------------------------------------
| ppf                                 |0.4760935726042412   |
------------------------------------------------------------
| pr                                  |0.12156238855280192  |
------------------------------------------------------------<jupyter_text>Lets see how these features correlate to our production<jupyter_code>graph_accuracy_importances(x, y, rf, data, rf.feature_importances_)<jupyter_output><empty_output><jupyter_text>Below we will see the differences between the regressor's predictions and the actual cumulative production<jupyter_code>from sklearn.ensemble import ExtraTreesRegressor
etr = ExtraTreesRegressor(random_state=86)
show_residusal(etr, train, test)
feature_importances(etr, x)<jupyter_output>Feature                             Score
------------------------------------------------------------
| easting                             |0.15796966128138012  |
------------------------------------------------------------
| permeability                        |0.10344420140720699  |
------------------------------------------------------------
| ppf                                 |0.22447668418572217  |
------------------------------------------------------------
| pr                                  |0.22172454753237059  |
------------------------------------------------------------<jupyter_text>Predict the following values:- porosity- perm- Young -<jupyter_code>def create_regressor(data:pd.DataFrame, x_names: list, y_name:str):
    ds = clean_dataset(data)
    x = ds[x_names]
    y = ds[y_name]
    
    rf = RandomForestRegressor(random_state=86)
    rf.fit(x, y)
    return rf
vals = ["easting","ppf", "pr", "pump rate (cubic feet/min)",'pump rate (cubic feet/min)', 'proppant weight (lbs)', 'well length']
dicto = {}
optim = list(x.drop(vals, axis=1).keys())
for var in optim:
    reg = create_regressor(big_df.drop("Name", axis=1), ["easting", "northing"], var)
    dicto[var] = reg<jupyter_output><empty_output><jupyter_text>The data here seems to be easily clusterable<jupyter_code>sns.scatterplot(data=big_df, x="easting", y="northing", hue="permeability")<jupyter_output><empty_output><jupyter_text>Stage 3: Choosing Parameters We want to look at the- [ ] Well Length- [ ] Frac stage count- [ ] Proppant per stage - Variable- [ ] Pump rate per stage - Variable<jupyter_code>!pip install tqdm
def iterate_over(df, var_name, step=100):
    min_ = min(df[var_name])
    max_ = max(df[var_name])
    for i in np.arange(min_, max_, step):
        yield i
def iter_hypers():
    values = {}
    for north in iterate_over(data, "northing", step=1000):
        values["northing"] = north
        for east in iterate_over(data, "easting", step=1000):
            values["easting"] = north
            # Predict 
            for key in dicto.keys():
                arr = np.array([north, east]).reshape(1, -1)
                values[key] = dicto[key].predict( arr )
            for well in iterate_over(data, "well length"):
                values["well length"] = well
                for frac in iterate_over(data, "frac stages"):
                    values["frac stages"] = frac
                    for prop in iterate_over(data, "ppf"):
                        values["ppf"] = prop
                        for pump in iterate_over(data, "pr"):
                            values["pr"] = pump
                            yield pd.DataFrame(values)
even_larger_daddy = pd.concat(iter_hypers())
even_large_daddy.to_csv("biggerPoppa.csv", index=False)<jupyter_output><empty_output><jupyter_text>Stage 4: Well location Graph all well locations<jupyter_code>sns.scatterplot(data=big_df, x="easting", y="northing", hue="easting")
plt.ylabel("Northing")
plt.xlabel("Easting")
plt.title("Location of Horizontal wells, starting position")
plt.legend('', frameon=False) # remove the legend
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13673530289465832, 'fraction_numerical': 0.014920919128618322, 'mean_word_length': 3.3979002624671915, 'pattern_counts': {'":': 6, '<': 110, '<?xml version=': 0, '>': 109, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '27', 'dir': 'jupyter-structured-clean-dedup', 'id': '134528', 'license': 'MIT', 'n_tokens_mistral': 5390, 'n_tokens_neox': 5046, 'n_words': 1521, 'path': 'Tools.ipynb', 'repo_name': 'fastfists/project', 'size': '16755'}