# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Programming For Data Analytics
#
# # Project
#
# # GMIT GRAD. DIP. DATA ANALYTICS
#
#
# ## PROBLEM STATEMENT
#
# Fictitious Pharmaceutical Company ‘Ollopa Pharmaceuticals’ manufactures the Active Pharmaceutical Ingredient (API) for its small molecule blockbuster drug ‘Asclepius’ by stepwise organic synthesis in a processing plant. The Manufacturing Science Group have been asked to optimise the production process to increase the yield of the active ingredient by monitoring the reaction endpoint using process analytical technology (such as 'at-line' or 'in-line' FTIR spectroscopy). The effect of changing process conditions and their relationship with reaction kinetics need to be monitored to achieve this. The Manufacturing Science Group have identified 4 Process Variables which they would like to investigate to see if there is any correlation between them and the reaction endpoint. The importance of identifying the reaction endpoint is critical for reactions where over reaction may lead to impurities being formed and the underreaction may lead to leftover starting materials whcih may not be purged from the process. 
#
# ## 1.CREATE THE VARIABLES
#
# **Critical Process Parameters**
#
# |RAP(%)     |Reaction Temp      |  pH                   | CO2                     |
# |:-------:  |:-----------------:|:--------------------: |:-----------------------:|
# | 0.2-99.9% |          <25C     | pH = 4.5-7.5          | 0-0.5%                  |
#
# #### **Variables Explained**
#
# **RAP** *Relative area %*: The % of Active ingredient produced during the course of the Reaction relative to the starting material (as the level of active increases, the level of strting matrial decreases. *We hope to determine the effects of the following variables on its rate*
#
# * **Reaction Temp:** Temp the reaction reaches (it is an exothermic reaction, i.e. one in which heat is produced)
# * **Pre crystallisation pH:** The pH of the batch during the reaction
# * **CO2 (%):** CO2 Levels produced during the reaction as a by-product
#
# The In-Process variables are the parameters which influence the yield and quality of the active ingredient. A manufacturing process which deviates from these parameters may have low yields and/or high impurity levels.
#
# For the purposes of creating a dataset we will investigate 100 data points of variables i.e. we will monitor a reaction over 10hrs, with the process variables being measured ten times per hour using in-line or at-line process analytical techniques.
#
# The effect or contribution (if any) of these variables on the reaction will be investigated.
#
#
#

# ***Loose* inspiration for the ideas behind the project came from papers such as the following:**

# [Model-based NIR spectroscopy implementation for in-line assay monitoring during a pharmaceutical suspension manufacturing process](https://www.sciencedirect.com/science/article/pii/S0378517318303521)
#
# [Monitoring of antibody-drug conjugation reactions with UV/Visspectroscopy](https://reader.elsevier.com/reader/sd/pii/S0168165618306552?token=6A2BBDF6C1063761A0D69C123693DF9BC4885E7D7FE6DA67488F45D0DAB6BF78C30729E380A16466D5ACDFD4F4D8A835)
#
# [Fiber-optic UV visible composition monitoring for process
# control of batch reactions](https://ac.els-cdn.com/S0169743998001051/1-s2.0-S0169743998001051-main.pdf?_tid=01e6088d-cc43-4b3f-8a07-4b35e9cfcbff&acdnat=1544731837_93d5b5e1af1778030dc063132dcfc019)

# For many years the US FDA (Food and Drug Administration) whichregulates the manufacture and sale of human medicines in the US and those manufactured abroad for the US market, has promoted the use of Process Analytical technologies for the monitoring of pharmaceutical processes https://www.fda.gov/downloads/drugs/guidances/ucm070305.pdf.
#
# *"This guidance is intended to describe a regulatory framework (Process Analytical Technology,
# PAT) that will encourage the voluntary development and implementation of innovative
# pharmaceutical development, manufacturing, and quality assurance. Working with existing
# regulations, the Agency has developed an innovative approach for helping the pharmaceutical
# industry address anticipated technical and regulatory issues and questions."*
#
# The industry itself has been slow to use this technology in commercial processes, mainly utilising it in scale-up and process development in small molecule manufacturing.
#
# With the emergence and growth of Biologics manufacturing, PAT use has increased in the industry

# **CREATION OF RANDOM VARIABLES**
#
# Import of python packages i.e. Numpy, Pandas, Scipy, odeint (a sub package of scipy for differential equations), Seaborn. The number of observations was set at 100 and a random seed generated. It was necessary to set a random seed in order to preserve the consistency of the data - since random variables were created each time the Jupyter Notebook was ran.
#
# The Time variable was created using the ```np.arange``` function which returns evenly spaced values in a given interval.
#
# The Relative Area Percent variable was created by trial and error using the ```odeint``` function. dy/dt vs time could then be plotted
#
# pH and CO2 Variables were created using a random normal function rounded to 2 decimal places
#
#
#

# +
import numpy as np
import pandas as pd
import scipy as sp
from scipy.integrate import odeint #odeint is a library which allows us to create and solve differential equations in python
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
# %matplotlib inline  

#number of observations

size = 100 

#Random seed set to maintain repeatability of the dataset
np.random.seed(876556788)

#List of variables
Time = np.arange(1,101) # Number of data points (with step size of 1) from start of reaction 

#RAP Will be calculated as a model to mimic a first order kinetic curve

#RxnTemp data will be created to correlate with the change in RAP over time.

pH = np.round(np.random.normal(loc=7.0, scale=0.25, size=size),2)

CO2 = np.round(np.random.normal(loc=0.25, scale=0.025, size=size),2)


# -

# **We will now generate a function to create a differential curve to model the increase in RAP during the reaction**
#
# The 'model' is a function which returns values of dy/dt for y and t. Initial value fro y was set at 0 and ```np.linspace``` returns evenly spaced numbers over a specified number of samples. ```odeint``` then calculates y(t). The resulting variables are then plotted as y V t defining y as RAP for the purposes of our dataset.

# +
#Credit: Adapted from (https://apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations)
#Placeholder for some simple code for graphing a differential curve to be adapted
#in order to mimic first order reaction  kinetics

# function that returns dy/dt
def model(y, t):
    dydt = -y + 100.0
    return dydt

# initial condition for y
y0 = 0

# time points
t = np.linspace(1, 10, 100)# set of data simulating time (100 data points over 10 hours)

# solve ODE
y = (odeint(model,y0,t)) #uses odeint function to model the relationship between dy/dt and t

# plot results
plt.plot(t,y)
plt.xlabel('time')
plt.ylabel('RAP(%)')
plt.show()
# -

# We now have a set of random variables of y which we can use to simulate the increase in Relative Area Percent (RAP) of product relative to starting material as the reaction progresses

#Output of y as a multidimensional array.
y

# 'y' Will need to be flattened into a single dimension array in order to be concatenated with the other variables using ```np.flatten()```
#
# We will now call this set of variables 'RAP'

RAP = y.flatten() #Converts multidimensional array to single dimensional array, allowing it be be concatenated with the 
# other random arrays in our dataset
RAP


# Another dy/dt v t model was created in order to model a change in reaction temperature over time. A lot of experimentation was used to find the right function for dy/dt. Eventually the ```np.cos```[](https://docs.scipy.org/doc/numpy/reference/generated/numpy.cos.html) command proved useful: 
#
# Since it is a simulated dataset, the objective is to create at least two variables which are correlated and RAP and Temperature were the two variables chosen, hence the need for them to be simulated using differential equations

# We will call Reaction temp 'RxnTemp' for the purposes of naming axes in the plot

# +
#Some experimentation to try to model Temp data to fit with RAP data
#Lots of random experimentation with help from (https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.math.html)

def model(z, t):
    dydt = (np.cos(1/t) - 1)*(-40)
    return dydt

# initial condition for y
z0 = 0

# time points
t = np.linspace(1, 10, 100) 

# solve ODE
z = (odeint(model,z0,t)) #uses odeint function to model the relationship between dy/dt and t

#RAP = pd.DataFrame(y)
# plot results
plt.plot(t,z)
plt.xlabel('time(hrs)')
plt.ylabel('RxnTemp')
plt.show()
# -

RxnTemp = z.flatten() #Converts mulitdimensional array to single dimensional array, allowing it be be concatenated with the 
# other random arrays in our dataset
RxnTemp

# All the variables were output as NP arrays for inspection together

#Output the Data as NUMPY arrays
print(('Time(hrs)=', Time), ('RAP =', RAP), ('pH =', pH), ('Reaction Temp(C) =', RxnTemp), ('CO2 =', CO2))

# ## 2. CREATE A DATAFRAME OF ALL VARIABLES
#
# **Each of the 4 numpy arrays will now be concatenated to create 1 dataframe**
#
# [1][Column names](https://docs.python.org/3.4/tutorial/datastructures.html#dictionaries)
#
# All of variables were concatenated to form one dataframe. Easily done using Pandas dataframe command.Names added to each column also. The purpose of this was to assemble a dataframe of simulated variables representing the hypothetical data gathered for the reaction of one batch of active ingredient. 100 data points are present, representing 100 measurements of process variables during the reaction
#

# +

columns = {"Time(hrs)": Time, "RAP": RAP, "pH": pH, "Reaction Temp(C)": RxnTemp,  "CO2": CO2} #Add labels to the columns

Batchdf = pd.DataFrame(columns, 
                  columns=["Time(hrs)", "RAP", "pH", "Reaction Temp(C)", "CO2"]) #Create the dataframe

Batchdf
# -

# **The batches.describe() function tabulates the summary statistics of the dataset and gives an overview of what the data looks like and how the Process Variables relate to eachother**
#
# Summary statistics give us a tabular view of the data and enable us to see how the variables might relate to eachother using the ```np.describe()``` command

Batchdf.describe() #Table of summary statistics

# ## 3. CREATE A SCATTERPLOT MATRIX OF THE  DATAFRAME
#

#[Credit] (https://github.com/spyder-ide/spyder/issues/7312)
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning) #Workaround to remove 'FutureWarnings' message which appears to be a bug. 

# **The following function ```sns.pairplot``` using the Seaborn package, creates a scatterplot matrix of variables and includes a density estimate (in the diagonal axis). It's a useful package to visualise data with multiple variables and get a quick idea of the relationships between variables. [11]**
#
# The ```sns.pairplot``` function is an excellent way of starting the process of Multivariate statistical analysis

# The Seaborn pairplot function enables a broad overview of how each variable relates to one another, showing scatter plots, linear relationships and Kernel density estimates. Kernel density estimates provide an estimate of how the data in each variable is distributed which may be useful to ascertain the types of subsequent statistical tests to be performed on the data

# %matplotlib inline
sns.pairplot(Batchdf, diag_kind='kde', plot_kws={'alpha':1.0});#scatterplots and kde matrices of all variables


# Simple scatter plots give an immediate indication if there is a linear relationship between variables, before correlation test are ran

# ## 4. EXAMINE THE  RELATIONSHIPS BETWEEN VARIABLES

# **Example scatter plot of Reaction Temp V RAP**
#
# From the seaborn scatter plot matrix above it can be seen that there appears to be a positive linear relationship between Reaction Temperature and RAP. The extent of this relationship may be further investigated by calculating the correlation coefficient

# +
g = sns.FacetGrid(Batchdf, height=10, )
g.map(plt.scatter, "Reaction Temp(C)", "RAP", s=50, linewidth=.5, edgecolor="red")
plt.rcParams.update({'axes.titlesize' : 24, 'axes.labelsize' : 20, 'lines.linewidth' : 3, 'lines.markersize' : 3, 'xtick.labelsize' : 16, 'ytick.labelsize' : 16, 'lines.linewidth' : 1, 'lines.markersize' : 10}) #Sets the plot font sizes etc
g.add_legend();


# -

# **Plot The linear regression models for each variable against RAP in order to determine if any significant linear relationships exist between variables**
#
# Using the Seaborn package to plot variables and output the regression equation. Of particular value is the R^2 value which tells us if there is correlation between variables.

# #### Linear Regression equation: 
#
# The relationship between variables in a linear relationship can be described by 
#
# \$Y_i = \beta_0 +\beta_1 X_i + \epsilon_i\$
#
# $Where:  \beta = slope, \epsilon$ = random error
#
# assuming y is a linear function of x but where for a fixed value of x the actual value of y is calculated by the linear model with a random error $\epsilon$  [13]

# +
# Get coefficients of linear fit
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning) #Workaround to remove 'FutureWarnings' message which appears to be a bug. 
# Linear regression for pH V RAP
slope, intercept, r_value, p_value, std_err = stats.linregress(Batchdf['pH'],Batchdf['RAP'])

# Set line label for legend
ax = sns.lmplot(x='pH', y='RAP', data=Batchdf, height=10)

plt.title("pH V RAP {}: y={:.1f}x+{:.1f} (p: {:.1f}, R^2: {:.2f})".format(ax, slope, intercept,  p_value, r_value)) #Adds title and linear equation
plt.show() # use of Seaborn to plot the relationships between data and generate the regression equation


# +
# get coefficients of linear fit

# Linear regression for Reaction Temp V RAP
slope, intercept, r_value, p_value, std_err = stats.linregress(Batchdf['Reaction Temp(C)'],Batchdf['RAP'])

# use line_kws to set line label for legend
ax = sns.lmplot(x='Reaction Temp(C)', y='RAP', data=Batchdf, height = 10)

plt.title("RAP V Reaction Temp(C) {}: y={:.1f}x+{:.1f} (p: {:.1f}, R^2: {:.2f})".format(ax, slope, intercept,  p_value, r_value))
plt.show()

# +
# get coefficients of linear fit

# Linear regression for IPC Result(%) V RAP
slope, intercept, r_value, p_value, std_err = stats.linregress(Batchdf['CO2'],Batchdf['RAP'])

# use line_kws to set line label for legend
ax = sns.lmplot(x='CO2', y='RAP', data=Batchdf, height = 10)
plt.title("RAP V CO2 {}: y={:.1f}x+{:.1f} (p: {:.1f}, R^2: {:.2f})".format(ax, slope, intercept,  p_value, r_value))
plt.show()
# -

# ## 5. DISCUSSION OF LINEAR REGRESSION

# It appears that the linear regression models indicate that the only relationship between variables which is linear is RAP V Reaction temp with a R^2 Value of 0.98, i.e. as the reaction progresses the temp increases. Display of the Correlation coefficients for each variable and it's relationship with RAP will help us determine the next steps in the analysis of the data.
#
# Uses the ```sp.stats.pearsonr()``` function
#
# The correlation coefficient is a measure of how closely related two variables are to eachother
# [Pearson Correlation](http://www.statsoft.com/Textbook/Statistics-Glossary/P/button/p#Pearson%20Correlation)

# +
 #Calculates the Pearson correlation coefficient for all variables V RAP and the *p* value for non-correlation
pHVRAP = sp.stats.pearsonr(pH,RAP)
RxnTempVRAP = sp.stats.pearsonr(RxnTemp,RAP)
CO2VRAP = sp.stats.pearsonr(CO2,RAP)

print ('pH V RAP: Correlation Coefficient', pHVRAP)
print ('Reaction Temperature V RAP: Correlation Coefficient', RxnTempVRAP,)
print ('CO2 V RAP: Correlation Coefficient', CO2VRAP)


# -

# The correlation coefficient for Reaction Temperature V RAP, calculated here using the 'scipy.stats.pearson' function shows a positive correlation between the two variables. This combined with the low *p* value confirms the hypothesis that the two variables are correlated.

# ## 6. CONCLUSION

# The relationship between temperature and RAP shows us that as the reaction progresses the temperature rises reaching a maximum as the reaction is complete. Knowing how as the reaction progresses, the temperature increases and levels off when complete will enable us to use the temperature as a guide when the reaction is complete and therefore when to sample for analysis to be confirmed offline in the laboratory.
#
# [9] Further analysis of the data using multivariate statistical techniques may enable statistical process monitoring. Such techniques can predict and determine the cause of out of control processes in many industries. The Hotellings T^2 statistic (the multivariate version of the students t-test statistic) and the Q statistic (test for randomness) are two methods which are used for the identification of out of control processes. [10] These two statistics are determined on data on which principal component analysis has been carried out. Principal Component analysis is a statistical technique as described by Jackson, 1991, where correlated variables are linearly transformed into 'principal components' which themselves are not correlated to eachother. PCA is in effect a technique for differentiating variables from eachother and in the process *reducing dimensions*. While there are only 4 variables in our data set and two which are correlated, our first principal component would be RAP/Temp as both are correlated variables. 
#
# PCA may be carried out using the ```sklearn.decomposition``` package and ```PCA``` class in Scikitlearn.
#
#

# ## 7. KEY LEARNINGS
#
# * Creation of random datasets, particularly ones which have distinct distributions enables the simulation of real world data heavy scenarios. This allows the testing of statistical tests and machine learning algorithms to model various outcomes.
# * Pandas Dataframes are a versatile data wrangling package for the manipulation and interrogation of datasets.
# * The python ecosystem is a rich environment with many and varied open source tools for Data Analytics.
# * Large datasets with multiple variables may be effectively analysed for latent realtionships using PCA and multivariate techniques.

# ## 8. **REFERENCES**
#
# 1. [Linspace](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linspace.html)
# 2. [Numpy](http://www.numpy.org/)
# 3. [Pandas](https://pandas.pydata.org/)
# 4. [Scipy](https://scipy.org/)
# 5. [scipy.integrate.odeint](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html)
# 6. [MatplotLib](https://matplotlib.org/index.html)
# 7. [Seaborn](https://seaborn.pydata.org)
# 8. [Github](https://help.github.com/articles/basic-writing-and-formatting-syntax/#links)
# 9. [A User’s Guide To Principal Components, <NAME>, 1991](https://pdfs.semanticscholar.org/70a9/769e8027b48bce5a9ce3a82d0ce88e0f0730.pdf) 
# 10. [Survey on data-driven industrial process monitoring and diagnosis](https://www.sciencedirect.com/science/article/pii/S1367578812000399)
# 11. [Python for Data Analysis; <NAME>.; O'Reilly; p.282](https://doc.lagout.org/programmation/python/Python%20for%20Data%20Analysis_%20Agile%20Tools%20for%20Real-World%20Data%20%5BMcKinney%202012-11-01%5D.pdf)
# 12. [sklearn.decomposition](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html)
# 13. [Applied Statistics and Probability for Engineers](http://www.um.edu.ar/math/montgomery.pdf) 
# 14. [Latex Styling](https://www.latex-project.org/)

StarcoderdataJupyterScriptsDedupFiltered