#################################################################
##### Title:   MapMark4GUI: A R Script for Running MapMark4GUI
##### Author:  Jason Shapiro
##### Date:    03/22/2018
#################################################################
#################################################################
#############To run, click Edit, Run All in the R Console 
#################################################################
#################################################################
###### References
###################################################################
##Ellefsen, K.J., 2017b, User’s guide for MapMark4—An R package for the probability calculations in three-part mineral resource assessments: Techniques and Methods.
##Singer, D., and Menzie, W.D., 2010, Quantitative Mineral Resource Assessments: An Integrated Approach: Oxford University Press.

######################################################
##### Reference usage comments 
########################################################
## This GUI is a tool to support implementation of the MapMark4 R package, many of the process scripts and resulting summaries are from or based off of Ellefsen's MapMark4 Package.
## MARK 3 PMF option section script is from Singer and Menzie, 2010.
## Bottom of this script is the PMF process script from the MapMark4 Package, to enable usage of the process with the PMF option. 


#################################################################
##Uploads required R packages 
#################################################################
library(MapMark4)
#install.packages('gWidgets')
library(gWidgets)
#install.packages('gWidgetstcltk')
library(gWidgetstcltk)
#install.packages("dplyr")
library(dplyr)

#################################################################
## Launches the GUI dialog
#################################################################win <- gwindow("MapMark Inputs")
win <- gwindow("MapMark Inputs")
grp_name <- ggroup(horizontal= FALSE, container = win)
grp2 <- ggroup(horizontal= FALSE, container = win)

#################################################################
## Sets frames in the GUI dialog
#################################################################
tmp <- gframe("Working Directory Input", container=grp_name)
info <- gframe("Background Information", container=grp_name)
input <- gframe("Input Files",horizontal= FALSE, container=grp_name)
sf <- gframe("Seed Information",container=grp_name)
modelI <- gframe("Model Information",horizontal= FALSE, container=grp_name)
runP <- gframe("Run Processes",container=grp_name)
ds <- gframe("Download Information",container=grp_name)
ps <- gframe("Plots",container=grp_name)
 
#################################################################
### Dialog development with buttons 
#################################################################

###################
##Working directory 
##################

### Creation of label for working directory button    
lbl_data_frame_name <- glabel(
  "Working Directory: ",
  container = tmp
)
### Creation of browse folder button
obj <- gfilebrowse("Select folder ",type = "selectdir", cont=tmp)
addhandlerchanged(obj, handler=function(h,...) 
	wdir1 <<- svalue(h$obj))

### sets button to change the working directory 
obj <- gbutton(
	text   = "Set Working Directory",
      container=tmp,
      handler = function(h,...)
	{
	setwd(wdir1)
	WD <<- getwd()
	cat (WD)
	})

###################
##Tract ID
##################
lbl_data_frame_name <- glabel(
  "Tract ID: ",
  container = info
)

obj1 <- gedit("Enter Tract ID here", container = info)
addhandlerchanged(obj, handler=function(h,...) 
	)

lbl_data_frame_name <- glabel(
  "Run ID: ",
  container = info
)

###################
##Run ID
##################
obj2 <- gedit("Enter Run ID here", container = info)
addhandlerchanged(obj, handler=function(h,...) 
	)
	

###################
##GTModel input
##################

lbl_data_frame_name <- glabel(
  "GTM Input: (GT deposits >= 20)",
  container = input
)
obj3 <- gfilebrowse(text = "Select a file...", type = "open", container = input)
addhandlerchanged(obj, handler=function(h,...) 
	)
	

###################
##NB estimates input
##################

lbl_data_frame_name <- glabel(
  "Estimates file for NegBinomial Type: ",
  container = input
)
obj4 <- gfilebrowse(text = "Select a file...", type = "open", container = input)
addhandlerchanged(obj, handler=function(h,...) 
	)
	#print out


###################
##Random seed
##################

obj <- gbutton(
	text   = "Create Random Seed",
      container=sf,
      handler = function(h,...)
	{
	seednum <<- sample(1:100, 1, replace=TRUE)
	obj9 <<- seednum
	lbl_data_frame_name <- glabel(
 	 seednum,
 	 container = sf
)
	})

###################
##Set specific seed
##################
obj <- gbutton(
	text   = "Set Specific Seed",
      container=sf,
      handler = function(h,...)
	{
	obj9 <<- gedit("",width = 4,container = sf)
	addhandlerchanged(obj, handler=function(h,...) {}
	)
	})

###################
##PMF Model
##################

lbl_data_frame_name <- glabel(
  "PMF model: ",
  container = modelI
)

obj5 <- gdroplist(c("Select a model", "NegBinomial", "MARK3"), container = modelI)
addhandlerchanged(obj, handler= function(h,...)
	)
##################	
## MARK3 estimates input	
##################
lbl_data_frame_name <- glabel(
  "Enter the N90, N50, N10, N05, and N01 percentiles [MARK3 Option Only] - \n If no estimates were made at N05 and N01, use the value from N10 ",
  container = modelI
)

obj6 <- gedit("NA", container = modelI)
addhandlerchanged(obj, handler=function(h,...) 
	)
	
##################	
## Truncate?	
##################

lbl_data_frame_name <- glabel(
  "Truncate?: ",
  container = modelI
)

obj7 <- gdroplist(c("Select true or false","TRUE", "FALSE"), container = modelI)
addhandlerchanged(obj, handler= function(h,...)
	)
	
##################	
## Distribution type	
##################
	
lbl_data_frame_name <- glabel(
  "Distribution type: (kde should be only used if the GTM has > 50 deposits)",
  container = modelI
)

obj8 <- gdroplist(c("Select a type","normal", "kde"), container = modelI)
addhandlerchanged(obj, handler= function(h,...)
	)

##################	
## Run PMF
##################


obj <- gbutton(
	text   = "1. Run Pmf",
      container=runP,
      handler = function(h,...)
	{
	


#################################################################
## PMF model run
#################################################################
	
#####################
## Sets PMF Parameters
######################
	tractname1 <<- svalue(obj1)
	testname1 <<- svalue(obj2)
	GTM <<- svalue(obj3)
	estimates <<- svalue(obj4)
	type <<- svalue(obj5)
	#type <<- strtoi(type , base = 0L)
	idp <<-c(svalue(obj6))
	trunc1 <<- svalue(obj7)
	#trunc1 <<- strtoi(type , base = 0L)
	typed <<- svalue(obj8)
	seednum<<- svalue(obj9)

################################################################
## MARK3 PMF Option Run- Based off of Singer, D., and Menzie, W.D., 2010, Quantitative Mineral Resource Assessments: An Integrated Approach: Oxford University Press.
#################################################################
	
	if (type == "MARK3") 
		{
		type <<- "Custom"
		idp5 <<- strsplit(idp, " ")
		idp6 <<- as.numeric(unlist(idp5))
		##adding a way to fill the rest of the line
		N9 <<- idp6[1]
		N5 <<- idp6[2]
		N1 <<- idp6[3]
		N05 <<- idp6[4]
		N01 <<- idp6[5]

 max1 <<- max(idp6)
nlist <<- c(0:max1)

P = c(0)

##########################################
##S Calculation 
##########################################
##calculation 

S1 <<- (1 + 2 * N9)

if (N9 == N5) 
{
   S2 <<- 1 
	
} else 
{
	S2 <<- 2 + 2 * (N5 - N9 - 1)
}

if (N5 == N1) 
{
   S3 <<- 1 
} else 
{
	S3 <- ( 2 + 2 * (N1 - N5 - 1))
}

if (N1 == N05) 
{
	S4 <<- 1 
} else 
{
	S4 <<- ( 2 + 2 * (N05 - N1 - 1))
}

if (N01 == N05) 
{
	S5 <<- 1 
}else 
{
	S5 <<- ( 2 + 2 * (N01 - N05 - 1))
}

#######################
##Setting Variables
######################
sum <- 0.0
expect <-0.0

###############################################
##For statements - general MapMark loop run
##################################################
for(n in 1:(N01 + 1))
{

	if (n < (N9 + 1)) 
	{
		
		P[n] <- 0.2/S1    ## first P, n = 1  [0.067]
	} else 
	{
		if (n == N9 + 1) 
		{
			P[n] = 0.1/S1 + 0.4/S2  #second P, n= 2 [0.233]
			if (N9 == N5) 
			{
				P[n] = P[n] + 0.4/S3
				if (N5 == N1) 
				{
					P[n] = P[n] + 0.05/S4
				}
				if (N5 == N05) 
				{
					P[n] <- P[n] + 0.04/S5
				}
			}	
		} else
		{
			if (n < N5 + 1)    
			{
				P[n] <- 0.8/S2  
				if (N9 == N5) 
				{
					P[n] <- 0.8/S3
				}
			} else 
			{
				if (n == N5 + 1) 
				{
					P[n] = 0.4/S2 + 0.4/S3
					if (N5 == N1) 
					{
						P[n] = P[n] + 0.05/S4
						if (N1 == N05) 
						{
							P[n] <- P[n] + 0.04/S5
						}
					}
				} else
				{
					if (n < N1 + 1) 
					{
						P[n] = 0.8/S3
						if (N5 == N1) 
						{
							P[n] <- 0.0
						}
					} else
					{
						if (n == N1 + 1) 
						{
							P[n] = 0.4/S3 + 0.05/S4
							if (N1 == N05) 
							{
								P[n] <- (P[n] + 0.04/S5 )    ## last P , n =5 [0.19]
							}
						} else
						{
							if (n < N05 + 1) 
							{
								P[n] <- 0.1/S4
							} else
							{
								if (n == N05 + 1) 
								{
									P[n] = 0.05/S4 + 0.04/S5
								} else
								{
									if (n < N01 + 1) 
									{
										P[n] <- 0.08/S5
									} else
									{
										if (n == N01 + 1) 
										{
											P[n] = 0.08/S5 + 0.01
											if (N01 == N05) 
											{
												P[n] <- P[n] + 0.025
											}
										}
									}
				
								}
							}
						}
					}
				}
			}
		}
	}
}

for(n in 1:(N01 + 1))
{
sum <- sum + P[n]
}

## changed the count to plus 1 to follow the others mabd added the last p , as the probability plus the rest of the sum
P[N01 + 1] <- P[N01 + 1] + (1.0 - sum)

#######################
## Setting Expect Variable
#######################
for(n in 1:(N01+1))
{
expect = expect + P[n] * n
}

rel6 <<- P
############################################################
### Runs PMF from MapMark4 package using the MARK 3 Option
#############################################################
		oPmf <<- NDepositsPmf2( type, list(nDeposits=nlist, relProbabilities=rel6))
##################################################
## Plots PMF model result and saves to .eps and .jpg files
##################################################	


		filename <- paste(testname1,"pmfPlot",".jpg", sep = "")
		jpeg(file=filename)
		plot(oPmf)
		dev.off() 


		filename <- paste(testname1,"pmfPlot",".eps", sep = "")
		plot(oPmf)

		postscript(file=filename)
		plot(oPmf)
		dev.off() 
		}
#################################################################
### Running PMF model with the MapMark4 Negative Binomial Option
#################################################################
		else if(type == "NegBinomial")
			{
			depEst <- read.csv(estimates)
###############################
## Runs PMF using MapMark4 and NB option
###################################
			oPmf <<- NDepositsPmf("NegBinomial", list(nDepEst = depEst))
			plot(oPmf)
############################################
##Plots PMF model and saves to .eps and .jpg files
############################################

			filename <- paste(testname1,"pmfPlot",".jpg", sep = "")
			jpeg(file=filename)
			plot(oPmf)
			dev.off() 


			filename <- paste(testname1,"pmfPlot",".eps", sep = "")
			postscript(file=filename)
			plot(oPmf)
			dev.off() 

			}

		
	}
 )

#########################
### Tonnage PDF model
#########################
obj <- gbutton(
	text   = "2. Run TonPdf",
      container=runP,
      handler = function(h,...)
	{
	gatm <- read.csv(GTM)
#################################################################
## Running PDF of the tonnage - using MapMark 4 package 
#################################################################
	oTonPdf <<- TonnagePdf(gatm, seed = seednum, pdfType = "normal", isTruncated = trunc1)
########################################
## Plots Ton PDF and saves to .eps and .jpg files
#######################################
	plot (oTonPdf)
	filename <- paste(testname1,"TonPdfPlot",".eps", sep = "")
	postscript(file=filename)
	plot(oTonPdf)
	dev.off() 


	filename <- paste(testname1,"TonPdfPlot",".jpg", sep = "")
	jpeg(file=filename)
	plot(oTonPdf)
	dev.off() 

	})
#########################
### Grade PDF model
#########################

obj <- gbutton(
	text   = "3. Run GradePdf",
      container=runP,
      handler = function(h,...)
	{
	gatm <- read.csv(GTM)
#################################################################
## Running PDF of the Grades - Using MapMark4 Package 
#################################################################	
gradePdf <<- GradePdf(gatm, seed = seednum, pdfType = "normal", isTruncated = trunc1)
########################################
## Plots Grades PDF and saves to eps and .jpg files
#######################################
	plot (gradePdf)
	filename <- paste(testname1,"GradPdfPlot",".eps", sep = "")
	postscript(file=filename)
	plot(gradePdf)
	dev.off() 

	filename <- paste(testname1,"GradPdfPlot",".jpg", sep = "")
	jpeg(file=filename)
	plot(gradePdf)
	dev.off() 

	})

####################################################
### Monte Carlo Simulation of undiscovered deposits 
####################################################
obj <- gbutton(
	text   = "4. Run Simulation",
      container=runP,
      handler = function(h,...)
	{
	gatm <- read.csv(GTM)
	oTonPdf <- TonnagePdf(gatm, seed = seednum, pdfType = typed, isTruncated = trunc1)
	gradePdf <- GradePdf(gatm, seed = seednum, pdfType = typed, isTruncated = trunc1)
#################################################################
## MARK3 PMF Option Run - for simulation- Based off of Singer, D., and Menzie, W.D., 2010, Quantitative Mineral Resource Assessments: An Integrated Approach: Oxford University Press.
#################################################################

	if (type == "MARK3") 
		{
		type <<- "Custom"
		idp5 <<- strsplit(idp, " ")
		idp6 <<- as.numeric(unlist(idp5))
		##adding a way to fill the rest of the line
		N9 <<- idp6[1]
		N5 <<- idp6[2]
		N1 <<- idp6[3]
		N05 <<- idp6[4]
		N01 <<- idp6[5]

 max1 <<- max(idp6)
nlist <<- c(0:max1)

P = c(0)


############
##calculation 
###############

S1 <<- (1 + 2 * N9)


if (N9 == N5) 
{
   S2 <<- 1 
	
} else 
{
	S2 <<- 2 + 2 * (N5 - N9 - 1)
}

if (N5 == N1) 
{
   S3 <<- 1 
} else 
{
	S3 <- ( 2 + 2 * (N1 - N5 - 1))
}

if (N1 == N05) 
{
	S4 <<- 1 
} else 
{
	S4 <<- ( 2 + 2 * (N05 - N1 - 1))
}

if (N01 == N05) 
{
	S5 <<- 1 
}else 
{
	S5 <<- ( 2 + 2 * (N01 - N05 - 1))
}


#####################
#Setting variables
####################
sum <- 0.0
expect <-0.0


#####################################
### General loop for the MARK3 run
#####################################  

for(n in 1:(N01 + 1))
{

	if (n < (N9 + 1)) 
	{
		
		P[n] <- 0.2/S1    ## first P, n = 1  [0.067]
	} else 
	{
		if (n == N9 + 1) 
		{
			P[n] = 0.1/S1 + 0.4/S2  #second P, n= 2 [0.233]
			if (N9 == N5) 
			{
				P[n] = P[n] + 0.4/S3
				if (N5 == N1) 
				{
					P[n] = P[n] + 0.05/S4
				}
				if (N5 == N05) 
				{
					P[n] <- P[n] + 0.04/S5
				}
			}	
		} else
		{
			if (n < N5 + 1)    
			{
				P[n] <- 0.8/S2  
				if (N9 == N5) 
				{
					P[n] <- 0.8/S3
				}
			} else 
			{
				if (n == N5 + 1) 
				{
					P[n] = 0.4/S2 + 0.4/S3
					if (N5 == N1) 
					{
						P[n] = P[n] + 0.05/S4
						if (N1 == N05) 
						{
							P[n] <- P[n] + 0.04/S5
						}
					}
				} else
				{
					if (n < N1 + 1) 
					{
						P[n] = 0.8/S3
						if (N5 == N1) 
						{
							P[n] <- 0.0
						}
					} else
					{
						if (n == N1 + 1) 
						{
							P[n] = 0.4/S3 + 0.05/S4
							if (N1 == N05) 
							{
								P[n] <- (P[n] + 0.04/S5 )    ## last P , n =5 [0.19]
							}
						} else
						{
							if (n < N05 + 1) 
							{
								P[n] <- 0.1/S4
							} else
							{
								if (n == N05 + 1) 
								{
									P[n] = 0.05/S4 + 0.04/S5
								} else
								{
									if (n < N01 + 1) 
									{
										P[n] <- 0.08/S5
									} else
									{
										if (n == N01 + 1) 
										{
											P[n] = 0.08/S5 + 0.01
											if (N01 == N05) 
											{
												P[n] <- P[n] + 0.025
											}
										}
									}
				
								}
							}
						}
					}
				}
			}
		}
	}
}

for(n in 1:(N01 + 1))
{
sum <- sum + P[n]
}

## changed the count to plus 1 to follow the others and added the last p , as the probability plus the rest of the sum
P[N01 + 1] <- P[N01 + 1] + (1.0 - sum)
#####################
### Setting expect variable
########################
for(n in 1:(N01+1))
{
expect = expect + P[n] * n
}

rel6 <<- P

###########################################################
#### PMF run for run simulation using MapMark4 package with the MARK3 option
#########################################################
		oPmf <<- NDepositsPmf2( type, list(nDeposits=nlist, relProbabilities=rel6))		}

###########################################################
#### PMF run for run simulation using MapMark4 package with the Neg B option
#########################################################
	else if(type == "NegBinomial")
		{
		depEst <- read.csv(estimates)
		oPmf <- NDepositsPmf("NegBinomial", list(nDepEst = depEst))
		}


####################################################
### Runs Simulation process using MapMark4 package
##################################################
	oSimulation <<- Simulation(oPmf, oTonPdf, gradePdf, seed = seednum)

###################################################
### Runs Total Tonnage PDF using MapMark4 Package 
#####################################################
	oTotalTonPdf <<- TotalTonnagePdf(oSimulation, oPmf, oTonPdf, gradePdf)

####################################################
### Plots Total Tonnage PDF and saves it as eps and .jpg files
####################################################
	plot (oTotalTonPdf)
	filename <- paste(testname1,"SimTotalPlot",".eps", sep = "")
	postscript(file=filename)
	plot(oTotalTonPdf)
	dev.off() 

	filename <- paste(testname1,"SimTotalPlot",".jpg", sep = "")
	jpeg(file=filename)
	plot(oTotalTonPdf)
	dev.off() 





#################################################################
#### Records simulation data file and temporarily saves it as a csv file.
#################################################################
 	filename2 <<- paste(testname1,"Datafile",".csv", sep = "")
	print(oSimulation,filename2)

#################################################################
## Develops the Simulation result EF file, file EF_5
#################################################################
dat = read.csv(filename2, header = TRUE)
LD <- length(names(dat))   # length of the columns in the data

#### Saving variables for calcualtion of EF 05 file from the above data file 
SimI <- dat[1]
NumD <- dat[2]
SimDI <- dat[3]
Ore<- dat[4]
Gran<-dat[LD]
G1<- dat[5]
Grades0<- G1
Tons0<- (Ore*dat[5]/100)
g <- 6
TonsList<-c("Tons5")
GradesList<-c("G5")

for (nam in 6:LD)
{
print (nam)
var<- paste("G",nam,sep="")
print (var)
assign(var,dat[g])


######################################
####calculate the ton for each mineral 
#######################################
varTon<- paste("Ton",nam,sep="")
assign(varTon,((Ore*(get(var)))/100))
##combine the grades and tons 
TonsList<-c(TonsList,varTon)

#########################
### Saves the grades
#########################
GradesList<-c(GradesList,var)
Grades0 <- cbind(Grades0,get(var))
names(Grades0)<-GradesList
Tons0<- cbind(Tons0,get(varTon))
names(Tons0)<-TonsList
g <- g + 1
}

TbNames <- names(dat)
TbLen <- length(TbNames)
MinStop <- TbLen -1 
minerals<- TbNames[5:MinStop]
OreN<- TbNames[4]
NamesMins<- sub('.grade', '', minerals)

NamesBegin <- TbNames[1:3]
Gangue <- TbNames[TbLen]
Gan<- sub('.grade', '', Gangue)


#######################
##Creates a joint table 
#########################
Cont <- cbind(SimI, NumD, SimDI, Ore,Grades0,Gran,Tons0) 
lenCont1 <- length(Cont)

lenMins<- length(NamesMins)
ContMath <- 5 + lenMins - 1
ContEndM <- ContMath + 2
ContBegin <- Cont[1:4]
ContMins <- Cont[5: ContMath]
ContMath1 <- ContMath + lenMins
ContEnd <- Cont [ContEndM: lenCont1]
NewCont <- cbind(ContBegin,ContMins,ContEnd)
MinTons<- paste(NamesMins,'_MetricTons')
OreN<- sub(".tonnage", "_MetricTons",OreN)
NamesMins<- paste(NamesMins,"_pct")
NameList12<- c(NamesBegin,OreN,NamesMins,Gan,MinTons ) 
lenCont2 <- length(NewCont)
con9 <- lenCont2 - 1
colnames(NewCont) <- NameList12 
NewCont <- NewCont[1:con9 ]

#############################################################
## Saving the simulation results - 05 EF file to a csv file.
#############################################################
filename5 <<- paste(testname1,"_05_SIM_EF",".csv", sep = "")
write.csv(NewCont, file = filename5)

#################################################################
#### Aggregation pivot calculcaiton, based on number of grades-   Contained Totals 
#################################################################
## make a pivot table 
cols<- names(Cont)

newtab<- Cont[1]  #creaitng a newtable so it can use consistent variable name
xy<- 1
LCont <- length(Cont)


NewT <-Cont[1:4]

#####################################################
## If statements based on number of grades- columns
####################################################

if (LCont == 7) 
{
NewT <- cbind(NewT,Cont[7])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1))
}
if (LCont == 9)
{
NewT <- cbind(NewT,Cont[8:9])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits), Ore = sum(Ore),Tons1 = sum(Tons1), Tons2 = sum(Tons2))
}
if (LCont == 11)
{
NewT <- cbind(NewT,Cont[9:11])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2","Tons3")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1), Tons2 = sum(Tons2), Tons3= sum(Tons3))
}
if (LCont == 13)
{
NewT <- cbind(NewT,Cont[10:13])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2","Tons3","Tons4")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1), Tons2 = sum(Tons2), Tons3= sum(Tons3), Tons4= sum(Tons4))
}
if (LCont == 15)
{
NewT <- cbind(NewT,Cont[11:15])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2","Tons3","Tons4","Tons5")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1), Tons2 = sum(Tons2), Tons3= sum(Tons3), Tons4= sum(Tons4), Tons5 = sum(Tons5))
}
if (LCont == 17)
{
NewT <- cbind(NewT,Cont[12:17])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2","Tons3","Tons4","Tons5","Tons6")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1), Tons2 = sum(Tons2), Tons3= sum(Tons3), Tons4= sum(Tons4), Tons5 = sum(Tons5), Tons6 = sum(Tons6))
}

if (LCont == 19)
{
NewT <- cbind(NewT,Cont[13:19])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2","Tons3","Tons4","Tons5","Tons6","Tons7")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1), Tons2 = sum(Tons2), Tons3= sum(Tons3), Tons4= sum(Tons4), Tons5 = sum(Tons5), Tons6 = sum(Tons6), Tons7 = sum(Tons7))
}

if (LCont == 21)
{
NewT <- cbind(NewT,Cont[14:21])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2","Tons3","Tons4","Tons5","Tons6","Tons7","Tons8")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1), Tons2 = sum(Tons2), Tons3= sum(Tons3), Tons4= sum(Tons4), Tons5 = sum(Tons5), Tons6 = sum(Tons6), Tons7 = sum(Tons7), Tons8 = sum(Tons8))
}

if (LCont == 23)
{
NewT <- cbind(NewT,Cont[15:23])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2","Tons3","Tons4","Tons5","Tons6","Tons7","Tons8","Tons9")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1), Tons2 = sum(Tons2), Tons3= sum(Tons3), Tons4= sum(Tons4), Tons5 = sum(Tons5), Tons6 = sum(Tons6), Tons7 = sum(Tons7), Tons8 = sum(Tons8), Tons9 = sum(Tons9))
}

if (LCont == 25)
{
NewT <- cbind(NewT,Cont[16:25])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2","Tons3","Tons4","Tons5","Tons6","Tons7","Tons8","Tons9","Tons10")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1), Tons2 = sum(Tons2), Tons3= sum(Tons3), Tons4= sum(Tons4), Tons5 = sum(Tons5), Tons6 = sum(Tons6), Tons7 = sum(Tons7), Tons8 = sum(Tons8), Tons9 = sum(Tons9), Tons10 = sum(Tons10))
}

if (LCont == 27)
{
NewT <- cbind(NewT,Cont[17:27])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2","Tons3","Tons4","Tons5","Tons6","Tons7","Tons8","Tons9","Tons10", "Tons11")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1), Tons2 = sum(Tons2), Tons3= sum(Tons3), Tons4= sum(Tons4), Tons5 = sum(Tons5), Tons6 = sum(Tons6), Tons7 = sum(Tons7), Tons8 = sum(Tons8), Tons9 = sum(Tons9), Tons10 = sum(Tons10), Tons11 = sum(Tons11))
}

if (LCont == 29)
{
NewT <- cbind(NewT,Cont[18:29])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2","Tons3","Tons4","Tons5","Tons6","Tons7","Tons8","Tons9","Tons10", "Tons11", "Tons12")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1), Tons2 = sum(Tons2), Tons3= sum(Tons3), Tons4= sum(Tons4), Tons5 = sum(Tons5), Tons6 = sum(Tons6), Tons7 = sum(Tons7), Tons8 = sum(Tons8), Tons9 = sum(Tons9), Tons10 = sum(Tons10), Tons11 = sum(Tons11), Tons12 = sum(Tons12))
}

if (LCont == 31)
{
NewT <- cbind(NewT,Cont[19:31])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2","Tons3","Tons4","Tons5","Tons6","Tons7","Tons8","Tons9","Tons10", "Tons11", "Tons12", "Tons13")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1), Tons2 = sum(Tons2), Tons3= sum(Tons3), Tons4= sum(Tons4), Tons5 = sum(Tons5), Tons6 = sum(Tons6), Tons7 = sum(Tons7), Tons8 = sum(Tons8), Tons9 = sum(Tons9), Tons10 = sum(Tons10), Tons11 = sum(Tons11), Tons12 = sum(Tons12), Tons13 = sum(Tons13))
}

if (LCont == 33)
{
NewT <- cbind(NewT,Cont[20:33])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2","Tons3","Tons4","Tons5","Tons6","Tons7","Tons8","Tons9","Tons10", "Tons11", "Tons12", "Tons13", "Tons14")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1), Tons2 = sum(Tons2), Tons3= sum(Tons3), Tons4= sum(Tons4), Tons5 = sum(Tons5), Tons6 = sum(Tons6), Tons7 = sum(Tons7), Tons8 = sum(Tons8), Tons9 = sum(Tons9), Tons10 = sum(Tons10), Tons11 = sum(Tons11), Tons12 = sum(Tons12), Tons13 = sum(Tons13), Tons14 = sum(Tons14))
}

if (LCont == 35)
{
NewT <- cbind(NewT,Cont[21:35])
colnames(NewT)<-c("SimIndex","NumDeposits","SimDepIndex","Ore","Tons1","Tons2","Tons3","Tons4","Tons5","Tons6","Tons7","Tons8","Tons9","Tons10", "Tons11", "Tons12", "Tons13", "Tons14", "Tons15")
Tb <- summarise(group_by(NewT,SimIndex),NumDep = mean(NumDeposits),Ore = sum(Ore), Tons1 = sum(Tons1), Tons2 = sum(Tons2), Tons3= sum(Tons3), Tons4= sum(Tons4), Tons5 = sum(Tons5), Tons6 = sum(Tons6), Tons7 = sum(Tons7), Tons8 = sum(Tons8), Tons9 = sum(Tons9), Tons10 = sum(Tons10), Tons11 = sum(Tons11), Tons12 = sum(Tons12), Tons13 = sum(Tons13), Tons14 = sum(Tons14), Tons15 = sum(Tons15))
}


nrc <- length(Tb) 
nc <- nrc -1
Tb <- Tb[1:nc]

##############################
## Setting table column names
#############################
TbNames <- names(dat)
TbLen <- length(TbNames)
TbLen<- TbLen -1
TbStart<- TbNames[1:2]
TbEnd<-TbNames[4:TbLen]
TbNames <- c(TbStart,TbEnd)


countNNN <- length(TbNames)
#Names123<- namelist4[3:countNNN]
NamesTT<- TbNames[3]
Names1 <- TbNames[1:2]
Names3 <- TbNames[4:countNNN]
Names3<- sub('.grade', '_mT', Names3)
NamesTT<- sub('.tonnage', '_mT', NamesTT)

NamesNew <- c(Names1,NamesTT,Names3)
colnames(Tb) <- NamesNew 

####################################################
#### Writes aggregated sim contained totals csv file
########################################################

filename6 <<- paste(testname1,"_07_SIM_Contained_Totals",".csv", sep = "")
write.csv(Tb, file = filename6)

})
############################################################
#### Sim Matrix
############################################################

obj <- gbutton(
	text   = "5. Run Sim Matrix",
      container=runP,
      handler = function(h,...)
	{
	plotMarginals(oTotalTonPdf)
	filename <- paste(testname1,"SimMatrixPlot",".eps", sep = "")
	postscript(file=filename)
	plotMarginals(oTotalTonPdf)
	dev.off() 

	filename <- paste(testname1,"SimMatrixPlot",".jpg", sep = "")
	jpeg(file=filename)
	plotMarginals(oTotalTonPdf)
	dev.off() 
	})

#################################################################
#### Review plots buttons
#################################################################

###################
## Plot PMF
###################
obj <- gbutton(
	text   = "Plot PMF ",
      container=ps,
      handler = function(h,...)
	{
plot (oPmf)
	})

###################
## Plot tonnage PDF
###################
obj <- gbutton(
	text   = "Plot Ton PDF ",
      container=ps,
      handler = function(h,...)
	{
plot (oTonPdf)
	})

###################
## Plot grade PDF
###################
obj <- gbutton(
	text   = "Plot Grade PDF ",
      container=ps,
      handler = function(h,...)
	{
	plot(gradePdf)
	})
###################
## Plot simulation results- total ton PDF
###################
obj <- gbutton(
	text   = "Plot Simulation Results",
      container=ps,
      handler = function(h,...)
	{
	plot(oTotalTonPdf)


	})
###################
## Plot sim matrix
###################
obj <- gbutton(
	text   = "Plot Sim Matrix",
      container=ps,
      handler = function(h,...)
	{
	plotMarginals(oTotalTonPdf) 
	})

####################################################################  Download information section
#################################################################


#########################################
### Sets download parmmeters and summary
#########################################

obj <- gbutton(
	text   = "Download Parameters & Summary",
      container=ds,
      handler = function(h,...)
	{
	date1 <<- Sys.Date()
	date2 <<- format(date1,"%a %b %d %Y")
	time1 <<- Sys.time()
	time2 <<- format(time1, "%X ")
	if (estimates == "Select a file..." ) { e2 <<- "None" }
	if (estimates != "Select a file..." ) { e2 <<- estimates }
	RS1 <- rbind (date2,time2, WD, tractname1, testname1, GTM, e2 , seednum, type, idp, trunc1, typed)
	rownames(RS1) <- c("Date", "Time", "Working Directory", "Tract ID", "Run ID", "GT Model Input", "NegBi Estimates Input", "Seed number", "PMF model", "N[9] N[5] N[1] N[05] N[01] -MARK3", "Truncate?", "Distribution Type")
	
#####################################
## Saving the parameters to a csv file
#####################################	
	RS2 <- paste(testname1,"_01_InputParameters.csv", sep = "")
	write.csv(RS1, file = RS2 , row.names=TRUE)

###############################################
### Creates and saves the MapMark4 summary file
##############################################

filename <- paste(testname1,"_Summary.txt", sep ="")

cat("Run Parameters:\n", file=filename )
capture.output( RS1, file=filename , append=TRUE)
cat("\n \n", file=filename , append=TRUE )
capture.output( summary(oTonPdf), file=filename , append=TRUE)
cat("\n \n", file=filename  , append=TRUE)
capture.output( summary(gradePdf), file=filename , append=TRUE)
cat("\n \n", file=filename , append=TRUE )
capture.output( summary(oPmf), file=filename , append=TRUE)
cat("\n \n", file=filename , append=TRUE )
capture.output( summary(oTotalTonPdf), file=filename , append=TRUE)
})

#########################################
### Sets download PMF stats
#########################################
obj <- gbutton(
	text   = "Download PMF Stats",
      container=ds,
      handler = function(h,...)
	{
###############
### PMF mean
################
	meanopmf <<- oPmf[7]
###############
### PMF variance 
################
	varopmf <<- oPmf[8]
	statsopmf <<- cbind(meanopmf,varopmf)
	if (type == "NegBinomial"){probs<<-oPmf[5]}
	if (type == "Custom") {probs <<- oPmf[2]}	
	
	

#################################
### Records the PMF probs file
##################################
	opfile2 <- paste(testname1,"_03_PMF_Probs.csv", sep = "")
	write.csv(probs, file = opfile2 , row.names=TRUE)

###################################
## If PMF option is Negative Binominal - pmf file saving
##################################
	if (type == "NegBinomial"){         
	dat1 <<- read.csv(opfile2 , header = TRUE)
	dat1[1] <- dat1[1] - 1
	colnames(dat1) <- c("NDeposits", "RelProbs") 
	write.csv(dat1, file = opfile2 , row.names=TRUE)			
	}
###################################
## If PMF option is MARK3 - pmf file saving
##################################
	if (type == "Custom"){         
	dat1 <<- read.csv(opfile2 , header = TRUE)
	dat1 <<- dat1[2:3]
	colnames(dat1) <- c("NDeposits", "RelProbs") 
	write.csv(dat1, file = opfile2 , row.names=TRUE)			
	}


##################################
##  Establishes entropy information
###################################
enf12 <- c()

opfile25 <- paste(testname1,"_03_PMF_Probs.csv", sep = "")
bb <- read.csv(opfile25)
for (hh in bb$RelProbs)
{


if (hh == 0)
{
print (hh)
en1 <- 0
enf12 <- c(enf12,en1)
}

if (hh > 0)
{
print (hh)
ln1 <- log(hh)
print (ln1)
en1 <<- (-hh * ln1)
print (en1)
enf12 <- c(enf12,en1)
}
}
InformationEntropy <<- sum(enf12)
print (InformationEntropy )
print (opfile25)
statsopmf <<- cbind(meanopmf,varopmf,InformationEntropy )

###############################
###Records PMF stats
###############################
	opfile <- paste(testname1,"_04__PMF_Stats.csv", sep = "")
	rownames(statsopmf) <- "stats"
	write.csv(statsopmf, file = opfile , row.names=TRUE)


})




#################################################################
### Sets Download Simulation and GTM stats button
#################################################################
obj <- gbutton(
	text   = "Download Sim and GTM Stats",
      container=ds,
      handler = function(h,...)
	{
#################################################################
## Simulation EF Stats 
#################################################################
## sim model Stats 
Rsim<- read.csv(filename5) ### reads in the sim data results- ef 5 file
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

#####################################
## Create means for each variable 
##################################
OMeans <- mean(v5)

for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- mean(v5)

OMeans<- cbind(OMeans,OM)
yy<- yy + 1
}
n <- length(OMeans)
#n <- n-1
OMeans <- OMeans[6:n]  # final means for each mineral and ore
#################################################################
## sim model max 
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)




####################################
## Create max for each variable 
####################################
OMaxs <- max(v5)

for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- max(v5)

OMaxs<- cbind(OMaxs,OM)
yy<- yy + 1
}
OMaxs<- OMaxs[6:n] # final maximums for each mineral and ore

#################################################################
## sim model min 
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###############################
## Create min for each variable 
############################## 
OMins <- min(v5)

for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- min(v5)

OMins<- cbind(OMins,OM)
yy<- yy + 1
}
OMins<- OMins[6:n]  # final minimums for each mineral and ore


#################################################################
## sim model median 
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###########################################
## Create median for each variable 
########################################
OMeds <- median(v5)

for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- median(v5)

OMeds<- cbind(OMeds,OM)
yy<- yy + 1
}
OMeds<- OMeds[6:n]   # final medians for each mineral and ore

#################################################################

## sim model standard deviations 
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

####################################
## Create standard deviation for each variable 
#############################################
OSds <- sd(v5)

for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- sd(v5)

OSds <- cbind(OSds,OM)
yy<- yy + 1
}
OSds<- OSds[6:n] # final standard deviations for each mineral and ore

#################################################################
## sim model Q99s
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###########################################
## 99 percentile for each variable 
##########################################
Q99s <- quantile(v5, c(.99))
for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.99))

Q99s <- cbind(Q99s,OM)
yy<- yy + 1
}
Q99s <- Q99s[6:n]  # final 99 percentiles for each mineral and ore

#################################################################
## sim model Q90s
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)


###########################################
## 90 percentile for each variable 
##########################################
Q90s <- quantile(v5, c(.90))
for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.90))

Q90s <- cbind(Q90s,OM)
yy<- yy + 1
}
Q90s <- Q90s[6:n]  # final 90 percentiles for each mineral and ore

#################################################################
## sim model Q80s
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###########################################
## 80 percentile for each variable 
##########################################  
Q80s <- quantile(v5, c(.80))
for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.80))

Q80s <- cbind(Q80s,OM)
yy<- yy + 1
}
Q80s <- Q80s[6:n]  # final 80 percentiles for each mineral and ore

#################################################################

## sim model Q70s
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###########################################
## 70 percentile for each variable 
##########################################
Q70s <- quantile(v5, c(.70))
for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.70))

Q70s <- cbind(Q70s,OM)
yy<- yy + 1
}
Q70s <- Q70s[6:n] # final 70 percentiles for each mineral and ore

#################################################################

## sim model Q60s
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###########################################
## 60 percentile for each variable 
##########################################
Q60s <- quantile(v5, c(.60))
for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.60))

Q60s <- cbind(Q60s,OM)
yy<- yy + 1
}
Q60s <- Q60s[6:n]  # final 60 percentiles for each mineral and ore

#################################################################
## sim model Q50s
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)



###########################################
## 50 percentile for each variable 
##########################################
Q50s <- quantile(v5, c(.50))
for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.50))

Q50s <- cbind(Q50s,OM)
yy<- yy + 1
}
Q50s <- Q50s[6:n] # final 50 percentiles for each mineral and ore

#################################################################

## sim model Q40s
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###########################################
## 40 percentile for each variable 
##########################################
Q40s <- quantile(v5, c(.40))
for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.40))

Q40s <- cbind(Q40s,OM)
yy<- yy + 1
}
Q40s <- Q40s[6:n] # final 40 percentiles for each mineral and ore

#################################################################

## sim model Q30s
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###########################################
## 30 percentile for each variable 
##########################################
Q30s <- quantile(v5, c(.30))
for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.30))

Q30s <- cbind(Q30s,OM)
yy<- yy + 1
}
Q30s <- Q30s[6:n]  # final 30 percentiles for each mineral and ore

#################################################################
## sim model Q20s
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###########################################
## 20 percentile for each variable 
########################################## 
Q20s <- quantile(v5, c(.20))
for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.20))

Q20s <- cbind(Q20s,OM)
yy<- yy + 1
}
Q20s <- Q20s[6:n]  # final 20 percentiles for each mineral and ore


#################################################################
## sim  model Q10s
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###########################################
## 10 percentile for each variable 
##########################################
Q10s <- quantile(v5, c(.10))
for (g in namelist4){

print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.10))

Q10s <- cbind(Q10s,OM)
yy<- yy + 1
}
Q10s <- Q10s[6:n]  # final 10 percentiles for each mineral and ore
#################################################################
## Sim model Q01s
Rsim<- read.csv(filename5)
namelist4 <- names(Rsim)
yy <- 1
v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###########################################
## 01 percentile for each variable 
##########################################
Q01s <- quantile(v5, c(.01))
for (g in namelist4){
print (g)
print (yy)
v4NA <- na.omit(Rsim[yy])
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.01))

Q01s <- cbind(Q01s,OM)
yy<- yy + 1
}
Q01s <- Q01s[6:n] # final 01 percentiles for each mineral and ore
#################################################################
##Create stats list
#################################################################
StatsList <- cbind(OMeans,OMaxs,OMins,OMeds,OSds,Q01s, Q10s, Q20s, Q30s, Q40s, Q50s, Q60s, Q70s, Q80s, Q90s, Q99s)
colnames(StatsList) <- c("Means", "Max", "Min", "Median", "STD", "P99", "P90", "P80", "P70", "P60", "P50", "P40", "P30", "P20", "P10", "P1")

#Rsim<- read.csv(GTM)
namelist5 <- names(Rsim)

xx <<- (n- 1)
rownames(StatsList) <- namelist5[5:xx] 


##############################################
##Downlaod Sim stats to csv file
############################################
Stats1 <<- paste(testname1,"_06_SIM_EF_Stats",".csv", sep = "")
write.csv(StatsList, file = Stats1, row.names=TRUE)

#################################################################### contained total stats 
#################################################################

Rsim<- read.csv(filename6)
n8 <<- ncol(Rsim)
n8 <<-( n8 + 1)
namelist4 <- names(Rsim)
yy <- 1

Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[1]
#v4NA <- na.omit(Rsim[1])
v5<- unlist(v4NA)
v5<- as.numeric(v5)

##############################################
## Create means for each variable 
##############################################
OMeans <- mean(v5)

for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- mean(v5)

OMeans<- cbind(OMeans,OM)
yy<- yy + 1
}
OMeans<- OMeans[5:n8]  # final means for each mineral and ore
#################################################################

## contained model max 
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
#################################
## Create max for each variable 
###################################
OMaxs <- max(v5)

for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- max(v5)

OMaxs<- cbind(OMaxs,OM)
yy<- yy + 1
}
OMaxs<- OMaxs[5:n8]   # final maxinums for each mineral and ore

#################################################################

## contained model min 
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###################################
## Create min for each variable
################################## 
OMins <- min(v5)

for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- min(v5)

OMins<- cbind(OMins,OM)
yy<- yy + 1
}
OMins<- OMins[5:n8]  # final mininums for each mineral and ore

#################################################################

## contained model med 
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

##############################################
## Create median for each variable 
#############################################
OMeds <- median(v5)

for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- median(v5)

OMeds<- cbind(OMeds,OM)
yy<- yy + 1
}
OMeds<- OMeds[5:n8] # final medians for each mineral and ore


#################################################################
## contained model STD Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

##################################
## Create standard deviations for each variable 
#####################################
OSds <- sd(v5)

for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- sd(v5)

OSds <- cbind(OSds,OM)
yy<- yy + 1
}
OSds<- OSds[5:n8]  # final standard deviations for each mineral and ore

#################################################################

## contained model Q99s
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

#################################
## Create 99 percentile for each variable 
#################################
Q99s <- quantile(v5, c(.99))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.99))

Q99s <- cbind(Q99s,OM)
yy<- yy + 1
}
Q99s <- Q99s[5:n8]   # final 99 percentiles for each mineral and ore

#################################################################
## Contained model Q90s
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

#########################################
## Create 90 percentiles for each variable 
########################################## 
Q90s <- quantile(v5, c(.90))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.90))

Q90s <- cbind(Q90s,OM)
yy<- yy + 1
}
Q90s <- Q90s[5:n8]   # final 90 percentiles for each mineral and ore

#################################################################

## Contained model Q80s
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

####################################
## Create 80 percentiles for each variable 
#######################################
Q80s <- quantile(v5, c(.80))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.80))

Q80s <- cbind(Q80s,OM)
yy<- yy + 1
}
Q80s <- Q80s[5:n8]  # final 80 percentiles for each mineral and ore

#################################################################

## Contained model Q70s
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
###################################
## Create 70 percentiles for each variable
################################### 

Q70s <- quantile(v5, c(.70))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.70))

Q70s <- cbind(Q70s,OM)
yy<- yy + 1
}
Q70s <- Q70s[5:n8] # final 70 percentiles for each mineral and ore

#################################################################

## Contained model Q60s
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

######################################
## Create 60 percentiles for each variable 
#####################################
Q60s <- quantile(v5, c(.60))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.60))

Q60s <- cbind(Q60s,OM)
yy<- yy + 1
}
Q60s <- Q60s[5:n8] # final 60 percentiles for each mineral and ore


#################################################################

## contained model Q50s
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

##################################################
## Create 50 percentiles for each variable 
####################################################
Q50s <- quantile(v5, c(.50))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.50))

Q50s <- cbind(Q50s,OM)
yy<- yy + 1
}
Q50s <- Q50s[5:n8]  # final 50 percentiles for each mineral and ore

#################################################################
## Contained  model Q40s
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

##########################################
## Create 40 percentiles for each variable 
##########################################
Q40s <- quantile(v5, c(.40))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.40))

Q40s <- cbind(Q40s,OM)
yy<- yy + 1
}
Q40s <- Q40s[5:n8]  # final 40 percentiles for each mineral and ore

#################################################################

## Contained model Q30s
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###########################################
## Create 30 percentiles for each variable 
#############################################
Q30s <- quantile(v5, c(.30))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.30))

Q30s <- cbind(Q30s,OM)
yy<- yy + 1
}
Q30s <- Q30s[5:n8]  # final 30 percentiles for each mineral and ore

#################################################################

## Contained model Q20s
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

##########################################
## Create 20 percentiles for each variable 
#########################################
Q20s <- quantile(v5, c(.20))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.20))

Q20s <- cbind(Q20s,OM)
yy<- yy + 1
}
Q20s <- Q20s[5:n8] # final 20 percentiles for each mineral and ore

#################################################################

## Contained model Q10s
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

#########################################
## Create 10 percentiles for each variable 
###########################################
Q10s <- quantile(v5, c(.10))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.10))

Q10s <- cbind(Q10s,OM)
yy<- yy + 1
}
Q10s <- Q10s[5:n8]  # final 10 percentiles for each mineral and ore

#################################################################

## Contained model Q01s
Rsim<- read.csv(filename6)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###############################################
## Create 01 percentiles for each variable 
###################################################
Q01s <- quantile(v5, c(.01))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.01))

Q01s <- cbind(Q01s,OM)
yy<- yy + 1
}
Q01s <- Q01s[5:n8]  # final 01 percentiles for each mineral and ore


#################################################################

##############################
## Calculating percent zero
###############################
Rsim<- read.csv(filename6)
#zero prob and mean
zero <-  length(which(Rsim[3]==0))
total <- length(which(Rsim[3]> -1))
PZero <- (zero/total)

######################################
## Calculating greater than or equal to mean
############################################
namelist4 <- names(Rsim)
yy <- 1

## create > mean  for each variable 
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

GOL <- length(which(v5 > mean(v5))) 
GEMS <- (GOL/total)



for (g in namelist4)
{
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

OM <- length(which(v5 >= mean(v5))) 
OM2 <- (OM/total)

GEMS <- cbind(GEMS,OM2)
yy<- yy + 1
}
GEMS <- GEMS[5:n8]  # final means for each mineral and ore


#################################################################
##Create stats list
StatsList <- cbind(OMeans,OMaxs,OMins,OMeds,OSds,Q01s, Q10s, Q20s, Q30s, Q40s, Q50s, Q60s, Q70s, Q80s, Q90s, Q99s, PZero,GEMS)
colnames(StatsList) <- c("Means", "Max", "Min", "Median", "STD", "P99", "P90", "P80", "P70", "P60", "P50", "P40", "P30", "P20", "P10", "P1","Prob of Zero", "Prob >= Mean" )

##########################################
### Setting row names for stats table
##########################################

Rsim<- read.csv(GTM)
namelist5 <- names(Rsim)
countNNN <- length(namelist5)
Names123<- namelist5[3:countNNN]
Names123<- paste(Names123,'_MetricTons')
rownames(StatsList) <- Names123

#######################################
##Downloading stats table to csv file
######################################
Stats1 <<- paste(testname1,"_08_SIM_Contained_Stats",".csv", sep = "")
write.csv(StatsList, file = Stats1, row.names=TRUE)



###################################################################GTM stats
#################################################################

## GTM model Stats 
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
n99 <<- ncol(Rsim)
n99 <<- (n99 + 1)

#####################################################
## Create means for each variable 
#######################################################
OMeans <- mean(v5)

for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- mean(v5)

OMeans<- cbind(OMeans,OM)
yy<- yy + 1
}
OMeans<- OMeans[4:n99]  # final means for each mineral and ore
#################################################################

## GTM model max 
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

#################################################################
## Create maxinums for each variable 
#################################################################
OMaxs <- max(v5)

for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- max(v5)

OMaxs<- cbind(OMaxs,OM)
yy<- yy + 1
}
OMaxs<- OMaxs[4:n99]  # final means for each mineral and ore

#################################################################

## GTM model min 
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)





#################################################################
## Create mininums for each variable 
#################################################################
OMins <- min(v5)

for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- min(v5)

OMins<- cbind(OMins,OM)
yy<- yy + 1
}
OMins<- OMins[4:n99]  # final mininums for each mineral and ore


#################################################################
## GT model med 
#################################################################
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

#################################################################
## Create median for each variable 
#################################################################
OMeds <- median(v5)

for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- median(v5)

OMeds<- cbind(OMeds,OM)
yy<- yy + 1
}
OMeds<- OMeds[4:n99]  # final medians for each mineral and ore




#################################################################
## GT model STD 
#################################################################
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

######################################
## Create standard deviations for each variable 
######################################
OSds <- sd(v5)

for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- sd(v5)

OSds <- cbind(OSds,OM)
yy<- yy + 1
}
OSds<- OSds[4:n99]  # final std for each mineral and ore

################################################################

## GT model Q99s
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

#############################################################
## Create 99 percentile for each variable 
##################################################
Q99s <- quantile(v5, c(.99))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.99))

Q99s <- cbind(Q99s,OM)
yy<- yy + 1
}
Q99s <- Q99s[4:n99]   # final 99 percentile for each mineral and ore

#################################################################

## GT model Q90s
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

#############################################
## Create 90 percentile for each variable 
##############################################
Q90s <- quantile(v5, c(.90))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.90))

Q90s <- cbind(Q90s,OM)
yy<- yy + 1
}
Q90s <- Q90s[4:n99]  # final 90 percentile for each mineral and ore

#################################################################

## GT model Q80s
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###############################################
## Create 80 percentile for each variable 
############################################### 
Q80s <- quantile(v5, c(.80))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.80))

Q80s <- cbind(Q80s,OM)
yy<- yy + 1
}
Q80s <- Q80s[4:n99]    # final 80 percentile for each mineral and ore

#################################################################

## GT model Q70s
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

############################################
## Create 70 percentile for each variable 
############################################
Q70s <- quantile(v5, c(.70))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.70))

Q70s <- cbind(Q70s,OM)
yy<- yy + 1
}
Q70s <- Q70s[4:n99]  # final 70 percentile for each mineral and ore

#################################################################

## GT model Q60s
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

###############################################
## Create 60 percentile for each variable 
###############################################
Q60s <- quantile(v5, c(.60))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.60))

Q60s <- cbind(Q60s,OM)
yy<- yy + 1
}
Q60s <- Q60s[4:n99]   # final 60 percentile for each mineral and ore

#################################################################

## GT model Q50s
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

######################################################
## Create 50 percentile for each variable 
######################################################
Q50s <- quantile(v5, c(.50))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.50))

Q50s <- cbind(Q50s,OM)
yy<- yy + 1
}
Q50s <- Q50s[4:n99]   # final 50 percentile for each mineral and ore

#################################################################

## GT model Q40s
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

#####################################################
## Create 40 percentile for each variable 
####################################################
Q40s <- quantile(v5, c(.40))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.40))

Q40s <- cbind(Q40s,OM)
yy<- yy + 1
}
Q40s <- Q40s[4:n99]  # final 40 percentile for each mineral and ore

#################################################################

## GT model Q30s
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

##########################################################
## Create 30 percentile for each variable 
#######################################################
Q30s <- quantile(v5, c(.30))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.30))

Q30s <- cbind(Q30s,OM)
yy<- yy + 1
}
Q30s <- Q30s[4:n99]   # final 30 percentile for each mineral and ore

#################################################################
## GT model Q20s
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

##################################################
## Create 20 percentile for each variable 
##################################################
Q20s <- quantile(v5, c(.20))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.20))

Q20s <- cbind(Q20s,OM)
yy<- yy + 1
}
Q20s <- Q20s[4:n99]   # final 20 percentile for each mineral and ore

#################################################################

## GT model Q10s
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

##################################################
## Create 10 percentile for each variable 
##################################################
Q10s <- quantile(v5, c(.10))
for (g in namelist4){

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.10))

Q10s <- cbind(Q10s,OM)
yy<- yy + 1
}
Q10s <- Q10s[4:n99]   # final 10 percetile for each mineral and ore

#################################################################

## GT model Q01s
Rsim<- read.csv(GTM)
namelist4 <- names(Rsim)
yy <- 1
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)

##############################################
## Create 01 percentile for each variable 
################################################
Q01s <- quantile(v5, c(.01))
for (g in namelist4)
{

print (g)
print (yy)
Rsim[is.na(Rsim)] <- 0
v4NA <- Rsim[yy]
v5<- unlist(v4NA)
v5<- as.numeric(v5)
OM <- quantile(v5, c(.01))

Q01s <- cbind(Q01s,OM)
yy<- yy + 1
}
Q01s <- Q01s[4:n99]  # final 01 percentile for each mineral and ore


#################################################################
##Create stats list
StatsList2 <- cbind(OMeans,OMaxs,OMins,OMeds,OSds,Q01s, Q10s, Q20s, Q30s, Q40s, Q50s, Q60s, Q70s, Q80s, Q90s, Q99s)
colnames(StatsList2) <- c("Means", "Max", "Min", "Median", "STD", "P99", "P90", "P80", "P70", "P60", "P50", "P40", "P30", "P20", "P10", "P1" )

#########################################
## Creates names for the GTM stats table
#########################################


Rsim<- read.csv(GTM)
namelist5 <- names(Rsim)
countNNN <- length(namelist4)
NamesGrades<- namelist5[4:countNNN]
OreG <-  namelist5[3]
NamesGrades<- paste(NamesGrades,'_pct')
OreT1<- paste(OreG ,'_MetricTons')
Names123 <- c(OreT1,NamesGrades)
rownames(StatsList2) <- Names123


#################################################################
## Saves the GTM stats to a csv file
################################################################

Stats2 <<- paste(testname1,"_02_GTM_Stats",".csv", sep = "")
write.csv(StatsList2, file = Stats2, row.names=TRUE)	


#####################################################
### Removes unneeded intermediate data file to support organization
###############################################################

if (file.exists(filename2))
 {file.remove(filename2)}


	})



#################################################################
###  MapMark4 PMF function process code - NDepositsPmf.R from Ellesfen's MapMark4 code- added to this script - to allow MARK3 PMF option to be used with it
#################################################################



CalcNegBinomialPmf <- function( nDepEst,
                                nGridPoints = 40,
                                power = 1.0,
                                isOptimized = TRUE,
                                probRightTail = 0.001 ){

  CheckInput <- function(df){

    colNames <- c("Name", "Weight", "N90", "N50", "N10")
    if(any(colnames(df) != colNames)){
      stop( sprintf( "Dataframe with the number of deposit estimates has\n"),
            sprintf( "one or more incorrect column names.\n"),
            sprintf( "The column names must be \"Name Weight N90 N50 N10\".\n"),
            call. = FALSE )
    }

    if(any(is.na(df))){
      stop( sprintf( "Dataframe with the number of deposit estimates has\n"),
            sprintf( "one or more missing entries.\n"),
            call. = FALSE )
    }
    if(any(df[, -1] < 0)){
      stop( sprintf( "Dataframe with the number of deposit estimates has\n"),
            sprintf( "one or more negative values.\n"),
            call. = FALSE )
    }
    for(i in 1:nrow(df)){
      if(!(df[i, "N90"] <= df[i, "N50"] && df[i, "N50"] <= df[i, "N10"])){
        stop( sprintf( "Row %4d of dataframe with the number of deposit\n", i),
              sprintf( "estimates is mis-ordered.\n"),
              call. = FALSE )
      }
    }

  }

  CalcParms <- function(nDepEst, nGridPoints, power, isOptimized){

    fn <- function(params, nDepEst, power){

      prevOp <- options(warn = -1)

      theProb <- params[1]
      theSize <- params[2]
      nDeposits <- qnbinom(c(0.90, 0.50, 0.10), size = theSize,
                           prob = theProb, lower.tail = FALSE)
      cost <- 0.0
      for(k in 1:nrow(nDepEst)){
        absDiff <- abs(nDeposits - nDepEst[k, c("N90", "N50", "N10")])
        currentCosts <- nDepEst[k, "Weight"] * (absDiff)^power
        cost <- cost + sum(currentCosts)
      }

      options(prevOp)

      return(cost)
    }

    lowBnd <- sum(nDepEst$N90 * nDepEst$Weight) / sum(nDepEst$Weight)
    upBnd <- sum(nDepEst$N10 * nDepEst$Weight) / sum(nDepEst$Weight)

    mean_seq <- seq(from = lowBnd + 0.01, to = upBnd, length.out = nGridPoints)

    sdMax <- upBnd - lowBnd
    varMax <- sdMax^2

    cost <- matrix(0, nrow = nGridPoints, ncol = nGridPoints)
    theProb <- matrix(NA_real_, nrow = nGridPoints, ncol = nGridPoints)
    theSize <- matrix(NA_real_, nrow = nGridPoints, ncol = nGridPoints)

    for(i in 1:nGridPoints){

      var_seq <- seq(from = mean_seq[i] + 0.01,
                     to = varMax,
                     length.out = nGridPoints)

      for(j in 1:nGridPoints){

        theProb[i,j] = mean_seq[i] / var_seq[j]
        theSize[i,j] <- mean_seq[i] * theProb[i,j] / (1 - theProb[i,j])

        cost[i,j] <- fn(c(theProb[i,j], theSize[i,j]), nDepEst, power)
      }
    }

    index <- which.min(cost)
    j <- index %/% nGridPoints + 1
    i <- index - (j-1) * nGridPoints

    if(isOptimized){
      tmp <- optim(c(theProb[i,j], theSize[i,j]), fn, NULL, nDepEst, power)
      if(tmp$convergence > 0){
        stop( sprintf( "Optimization failed\n"),
              call. = FALSE )
      } else{
        result <- list(theProb = tmp$par[1], theSize = tmp$par[2])
      }
    } else {
      result <- list(theProb = theProb[i,j], theSize = theSize[i,j])
    }

    return(result)
  }


  CheckInput(nDepEst)

  nDepEst$Name <- NULL


  nDepEst <- nDepEst[nDepEst$Weight > 0, ]

  params <- CalcParms(nDepEst, nGridPoints, power, isOptimized)


  x <- qnbinom( 1-probRightTail, size = params$theSize, prob = params$theProb )


  nDeposits <- 0:x
  pmf <- dnbinom( nDeposits, size = params$theSize, prob = params$theProb )

  areTooSmall <- (pmf < pmf[length(pmf)])
  nDeposits <- nDeposits[!areTooSmall]
  pmf <- pmf[!areTooSmall]


  pmf <- pmf/sum(pmf)

  return( list( nDeposits = nDeposits, probs = pmf ) )
}


CalcCustomPmf2 <- function( nDeposits, relProbabilities ) {

  if( length( nDeposits ) != length( relProbabilities ) ) {
    stop( sprintf( "Function CalcCustomPmf\n" ),
          sprintf( "The length of nDeposits must equal the length of relProbabilities.\n" ),
          call. = FALSE )
  }

  if( any( nDeposits < 0 ) ) {
    stop( sprintf( "Function CalcCustomPmf\n" ),
          sprintf( "All elements of nDeposits must be >= 0.\n" ),
          call. = FALSE )
  }
###########################################
### Change from orginal MapMark4 code - deletion 
####################################
 # if( any( diff( nDeposits ) <= 0 ) ) {
 #   stop( sprintf( "Function CalcCustomPmf\n" ),
 #         sprintf( "All elements of nDeposits must be in strictly ascending order.\n" ),
 #         call. = FALSE )
 # }

  if( any( relProbabilities < 0.0 ) ) {
    stop( sprintf( "Function CalcCustomPmf\n" ),
          sprintf( "All elements of relProbabilities must be >= 0.\n" ),
          call. = FALSE )
  }

  probs <- relProbabilities / sum( relProbabilities )

  return( list( nDeposits=nDeposits, probs=probs ) )

}

plot.NDepositsPmf <- function( object,
                               isMeanPlotted = FALSE,
                               areLinesAdded = TRUE,
                               isUsgsStyle = TRUE) {


  df <- data.frame(nDeposits = object$nDeposits,
                   probs = object$probs)

  p <- ggplot2::ggplot(df)

  if(isMeanPlotted == TRUE) {
    p <- p + ggplot2::geom_vline(xintercept = object$theMean, colour = "red")
  }

  if(areLinesAdded){
    for(i in 1:length(object$nDeposits)){
      df1 <- data.frame( x = rep.int(object$nDeposits[i], 2),
                        y = c(0, object$probs[i]))
      p <- p + ggplot2::geom_line(ggplot2::aes(x = x, y = y), data = df1)
    }
  }

  p <- p + ggplot2::geom_point(ggplot2::aes(x = nDeposits, y = probs)) +
    ggplot2::scale_x_continuous(name = "Number of undiscovered deposits",
                                limits = range(df$nDeposits) + c(-1,1)) +
    ggplot2::scale_y_continuous(name = "Probability",
                                limits = c(0, max(df$probs)))

  if(isUsgsStyle)
    p <- p + ggplot2::theme_bw()

  if(object$type == "Custom") {
    print(p)
  } else {

    df$accdf <- c(1, 1-cumsum(df$probs)[-length(df$probs)])

    a <- table(object$pmf.args$nDepEst$N90)
    b <- table(object$pmf.args$nDepEst$N50)
    c <- table(object$pmf.args$nDepEst$N10)

    df2 <- data.frame(nDeposits = c(as.integer(names(a)),
                                    as.integer(names(b)),
                                    as.integer(names(c))),
                      Percentiles = c(rep.int(90, length(a)),
                                      rep.int(50, length(b)),
                                      rep.int(10, length(c))),
                      Sizes = c(as.vector(a),
                                as.vector(b),
                                as.vector(c)))

    q <- ggplot2::ggplot() +
      ggplot2::geom_point(ggplot2::aes(x = nDeposits,
                                       y = Percentiles,
                                       size = Sizes),
                          data = df2, colour = "red", shape = 1, stroke = 1.1) +
      ggplot2::geom_point(ggplot2::aes(x = nDeposits,
                                       y = accdf * 100),
                          data = df) +
      ggplot2::scale_y_continuous(name = "Elicitation percentiles", limits = c(0,100),
                                  breaks = c(0, 10, 50, 90, 100)) +
      ggplot2::scale_x_continuous(name = "Number of undiscovered deposits") +
      ggplot2::scale_size(name = "Count", breaks = min(df2$Sizes):max(df2$Sizes))

    if(isUsgsStyle)
      q <- q + ggplot2::theme_bw()

    q <- q + ggplot2::theme(legend.position = c(1,1), legend.justification = c(1.1,1.1))

    if(isUsgsStyle) {
      p <- p + ggplot2::ggtitle("A.") +
        ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0,
                                                          face = "italic"))
      q <- q + ggplot2::ggtitle("B.") +
        ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0,
                                                          face = "italic"))
    } else {
      p <- p + ggplot2::ggtitle("(a)") +
        ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0))
      q <- q + ggplot2::ggtitle("(b)") +
        ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0))
    }

    grid::grid.newpage()
    grid::pushViewport(grid::viewport(layout=grid::grid.layout(1,2)))
    plot(p, vp=grid::viewport(layout.pos.row=1, layout.pos.col=1))
    plot(q, vp=grid::viewport(layout.pos.row=1, layout.pos.col=2))

  }

}


summary.NDepositsPmf <- function( object) {

  cat(sprintf("Summary of the pmf for the number of undiscovered deposits\n"))
  cat(sprintf("within the permissive tract.\n"))
  cat(sprintf("------------------------------------------------------------\n"))
  cat( sprintf( "Type: %s\n", object$type ))
  cat( sprintf( "Description: %s\n", object$description ))
  cat( sprintf( "Mean: %g\n", object$theMean ))
  cat( sprintf( "Variance: %g\n", object$theVar ))
  cat( sprintf( "Standard deviation: %g\n", sqrt(object$theVar) ))
  cat( sprintf( "Mode: %d\n", object$nDeposits[which.max(object$probs)] ))
  cat( sprintf( "Smallest number of deposits in the pmf: %d\n",
                min(object$nDeposits) ))
  cat( sprintf( "Largest number of deposits in the pmf: %d\n",
                max(object$nDeposits) ))
  cat( sprintf( "Information entropy: %g\n",
                max(object$entropy) ))

  cat( sprintf( "\n###############################################################\n"))
  cat(sprintf("\n\n\n\n"))

}

########################################
###  New function name NDepositsPmf2
########################################


NDepositsPmf2 <- function( type, pmf.args, description="" ) {


  CalcSummaryStats <- function( probs, nDeposits, base=exp(1) )
  {

    theMean <- sum( probs * nDeposits )
    theVar <- sum( probs * nDeposits^2 ) - theMean^2
    theEntropy <- - sum( probs * log( probs, base=base ) )

    return( list( theMean=theMean, theVar=theVar, theEntropy=theEntropy ) )

  }

  if(!any(type == c("NegBinomial", "Custom"))) {
    stop( sprintf( "Function NDepositsPmf\n" ),
          sprintf( "Argument type must be either NegBinomial or Custom.\n"),
          sprintf( "It is specified as %s\n", type),
          call. = FALSE )
  }

  pmf <- switch( type,
                 NegBinomial = do.call( CalcNegBinomialPmf, pmf.args ),
                 Custom = do.call( CalcCustomPmf2, pmf.args ) )
  nTotalDraws <- 20000

  stats <- CalcSummaryStats( pmf$probs, pmf$nDeposits )

  rval <- list( type = type,
                pmf.args = pmf.args,
                description = description,
                call = sys.call(),
                probs = pmf$probs,
                nDeposits = pmf$nDeposits,
                theMean = stats$theMean,
                theVar = stats$theVar,
                entropy = stats$theEntropy )

  class(rval) <- "NDepositsPmf"
  return(rval)
}