r_plots_part5_LTA_sensitivity_analysis

r_plots_part5_LTA_sensitivity_analysis.R

# Plots for PArt 5 LTA sensitivity analysis
library(sm)
library(pheatmap)
library(ggplot2)

# original Grid 
M1 = c(0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300)
M2 = c(0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.100, 0.220, 0.280, 0.250, 0.240, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300)
M3 = c(0.300, 0.300, 0.300, 0.300, 0.300, 0.215, 0.255, 0.215, 0.145, 0.275, 0.170, 0.240, 0.250, 0.250, 0.280, 0.290, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300)
M4 = c(0.300, 0.300, 0.300, 0.300, 0.300, 0.170, 0.270, 0.190, 0.190, 0.285, 0.250, 0.225, 0.275, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300)
M5 = c(0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300)
M6 = c(0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300)
inspection_matrix = t(data.frame(M1,M2,M3,M4,M5,M6))


# Random Grid 
random_set = seq(.200,.3,.005)
sample(random_set,1)
M1 = c(0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300)
M2 = c(0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.100, sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300)
M3 = c(0.300, 0.300, 0.300, 0.300, 0.300, sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), 0.300, 0.300, 0.300, 0.300, 0.300, 0.300)
M4 = c(0.300, 0.300, 0.300, 0.300, 0.300, sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), sample(random_set,1), 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300)
M5 = c(0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300)
M6 = c(0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300)
inspection_matrix = t(data.frame(M1,M2,M3,M4,M5,M6))


require(RColorBrewer)
#heatmap(inspection_matrix, Rowv=NA, Colv=NA,col = brewer.pal(9,"RdYlGn"),scale="none",xlab = "Longitudinal Plane", ylab =  "Circumfrential Plane", main ="Inspection Data")

#install.packages("pheatmap")
require(pheatmap)

cols = colorRampPalette(colors = c("blue", "white"))
pheatmap(inspection_matrix,          color = brewer.pal(9,"RdYlGn"),
         #breaks = c(0.3, .2, .1), 
         display_numbers = T,cluster_rows=F,cluster_cols=F, main = "Random Grid Structure #3 +0.100")

# Determine CTP 
CTPc = c(min(M1),min(M2),min(M3),min(M4),min(M5),min(M6))
CTPl = c(min(inspection_matrix[,1]),min(inspection_matrix[,2]),min(inspection_matrix[,3]),min(inspection_matrix[,4]),min(inspection_matrix[,5]),min(inspection_matrix[,6]),min(inspection_matrix[,7]),
         min(inspection_matrix[,8]),min(inspection_matrix[,9]),min(inspection_matrix[,10]),min(inspection_matrix[,11]),min(inspection_matrix[,12]),min(inspection_matrix[,13]),min(inspection_matrix[,14]),min(inspection_matrix[,15]),
         min(inspection_matrix[,16]),min(inspection_matrix[,17]),min(inspection_matrix[,18]),min(inspection_matrix[,19]),min(inspection_matrix[,20]),min(inspection_matrix[,21]),
         min(inspection_matrix[,22]))

# Plot Pit Profiles 
max = max(CTPl, na.rm=TRUE)
plot(x = 1:length(CTPl),CTPl, type = "l", pch = 19, col = "red", xlab = "No.", ylab = "t", main = "Longitudinal Critical Thickness Profile - Random Grid Structure #3 +0.100", cex.main=1.0,axes = FALSE)
axis(side = 1, at = 1:length(CTPl),labels = 1:length(CTPl), las =1,cex.axis=.7,tck = 1, lty = 2, col = "grey")
axis(side = 2, at = seq(0,max+.1,.1),labels = seq(0,max+.1,.1), las =0,cex.axis=.7,tck = 1, lty = 2, col = "grey")
box()


data_dist_1 = read.csv("C:/Users/Andrew.Bannerman/OneDrive - Shell/Documents/API CODES/level_2_LTA_grid_assumption_MAWPr_out_case_1.csv", header=T, stringsAsFactors = F)
data_dist_1$label_1 = "MAWPr_dist_1"
mean_dist_1 = mean(data_dist_1$x1)
median_dist_1 = median(data_dist_1$x1)
data_dist_2 = read.csv("C:/Users/Andrew.Bannerman/OneDrive - Shell/Documents/API CODES/level_2_LTA_grid_assumption_MAWPr_out_case_2.csv", header=T, stringsAsFactors = F)
data_dist_2$label_1 = "MAWPr_dist_2"
mean_dist_2 = mean(data_dist_2$x1)
median_dist_2 = median(data_dist_2$x1)
data_dist_3 = read.csv("C:/Users/Andrew.Bannerman/OneDrive - Shell/Documents/API CODES/level_2_LTA_grid_assumption_MAWPr_out_case_3.csv", header=T, stringsAsFactors = F)
data_dist_3$label_1 = "MAWPr_dist_3"
mean_dist_3 = mean(data_dist_3$x1)
median_dist_3 = median(data_dist_3$x1)

all_out = rbind(data_dist_1,data_dist_2,data_dist_3)

data = data_dist_3
  
require(ggplot2)
ggplot(data, aes(x=x1)) + 
  geom_histogram( binwidth=30, fill="#69b3a2", color="#e9ecef", alpha=0.9)+
  ggtitle("LTA Level 2 Grid Assumption - 10,000 Iteration Output Distribution")+
  scale_x_continuous(breaks=seq(round(min(data$x1),0),round(max(data$x1),0),50))+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  xlab("MAWPr")+
  ylab("Count")+
  geom_vline(xintercept=mean(data$x1),colour="#BB0000")+
  annotate("text", x = mean(data$x1)+100, y = 1800, label = paste("Mean MAWPr",round(mean(data$x1),2)),colour="#BB0000")+
  geom_vline(xintercept=3008.19)+
  annotate("text", x = 3008.19-100, y = 1800, label = "Actual Known Grid MAWPr\n3008.19")

# 1x plot
plot(data$x1, dnorm(data$x1, 2, 0.5), type = "l")
d <- density(data$x1) # returns the density data
plot(d) # plots the results
polygon(d, col="red", border="blue")

#+++++++++++
# https://towardsdatascience.com/how-to-find-probability-from-probability-density-plots-7c392b218bab
# density(x, ...)
## Default S3 method:
density(data$x1, bw = "nrd0", adjust = 1,
        kernel = c("gaussian", "epanechnikov", "rectangular",
                   "triangular", "biweight",
                   "cosine", "optcosine"),
        weights = NULL, window = kernel, width,
        give.Rkern = FALSE,
        n = 512, from, to, cut = 3, na.rm = FALSE, ...)



# Compare MAWPr distributions

# create value labels
value_labels <- factor(all_out$label_1, levels= c(1,2,3),
                labels = c("MAWPr_dist_1", "MAWPr_dist_2","MAWPr_dist_3"))

all_out$label_1 = factor(all_out$label_1)

sm.density.compare(all_out$x1, all_out$label_1, xlab="MAWPr")
title(main="MAWPr distribution (retain max loss - 30,000 random grid samples)\n Random grid samples selected from 3x loss to remaining ratio ranges (2,1,0.5)\n Distribution comparisons highlight the effect on MAWPr when there is thicker regions\n around local minimums - local reinforcement effects increase MAWPr",cex.main=.85)
text(mean_dist_1-100, 0.0025, labels = "loss/remaining ratio <=2", col = "red")
text(mean_dist_2-102, 0.0035, labels = "loss/remaining ratio <=1", col = "darkgreen")
text(mean_dist_3-100, 0.0053, labels = "loss/remaining ratio <=0.5", col = "blue")
mtext(side=3, line=2, at=0.07, adj=0, cex=0.7, "mysubtitle")
axis(side = 1, at = seq(1500,3500,250))