################################################################################
# Scripts for a grade and tonnage model
################################################################################

##
## Preparatory steps
##

# Load the package
library(MapMark4)

# Set the working directory
oldWd <- setwd("D:\\Datap-Qmra\\MapMark4\\Temp\\TestCalculationsGatm")

write.csv(ExampleGatm, file = "myGatm.csv", row.names = FALSE)
write.csv(ExampleDepEst1, file = "myDepEst.csv", row.names = FALSE)

##
## Probability calculations
##

# Generate the meta data
oMeta <- Metadata("PT001", "myGatm.csv", "myDepEst.csv",
                  seed = 7, tractName = "Lucky Strike",
                  depositModel = "Computer generated gatm with copper and gold",
                  personName = "Mary Doe",
                  otherInfo = "This is an example assessment.")
summary(oMeta)

# Generate the pdf that represents the ore tonnages and check it
gatm <- read.csv(getModelFilename(oMeta))
oTonPdf <- TonnagePdf(gatm, seed = getSeed(oMeta))
plot(oTonPdf)
summary(oTonPdf)

# Generate the pdf the represents the grades and check it
oGradePdf <- GradePdf(gatm, seed = getSeed(oMeta))
plot(oGradePdf)
summary(oGradePdf)

# Generate the pmf that represents the number of undiscovered deposits
# and check it
depEst <- read.csv(getDepEstFilename(oMeta))
oPmf <- NDepositsPmf("NegBinomial", list(nDepEst = depEst))
plot(oPmf)
summary(oPmf)

# Simulate undiscovered deposits in the permissive tract
oSimulation <- Simulation(oPmf, oTonPdf, oGradePdf, seed = getSeed(oMeta))

# Generate the pdf that represents the total ore and resource tonnages for all
# simulated undiscovered deposits in the permissive tract and check it
oTotalTonPdf <- TotalTonnagePdf(oSimulation, oPmf, oTonPdf, oGradePdf)
plot(oTotalTonPdf)
plotMarginals(oTotalTonPdf)
summary(oTotalTonPdf)

##
## Archive of the calculation results
##

# Summaries of the probability calculations
filename <- paste(getTractId(oMeta), "_summary.txt", sep = "")
oldWidth <- options(width = 120)
sink(filename)
summary(oMeta)
summary(oTonPdf)
summary(oGradePdf)
summary(oPmf)
summary(oTotalTonPdf)
sink()
options(width = oldWidth$width)

# Simulated undiscovered deposits, which are needed for
# subsequent assessment calculations that are not part of MapMark4
filename <- paste(getTractId(oMeta), "_Simulation.csv", sep = "")
print(oSimulation, filename)

# Instantiated classes
filename <- paste(getTractId(oMeta), ".dat", sep = "")
save(oMeta, oTonPdf, oGradePdf, oPmf, oSimulation, oTotalTonPdf, file = filename)

# Store the plots in PNG format
png(filename = paste(getTractId(oMeta), "_TonPdf.png", sep = ""), height = 240)
plot(oTonPdf)
dev.off()

png(filename = paste(getTractId(oMeta), "_GradePdf.png", sep = ""))
plot(oGradePdf)
dev.off()

png(filename = paste(getTractId(oMeta), "_Pmf.png", sep = ""), height = 240)
plot(oPmf)
dev.off()

png(filename = paste(getTractId(oMeta), "_TotalTonPdf.png", sep = ""), height = 240)
plot(oTotalTonPdf)
dev.off()

png(filename = paste(getTractId(oMeta), "_TotalTonPdf_Marg.png", sep = ""))
plotMarginals(oTotalTonPdf)
dev.off()


# Store the plots in postscript format
postscript(paste(getTractId(oMeta), "_TonPdf.ps", sep = ""), height = 3.5)
plot(oTonPdf)
dev.off()

postscript(paste(getTractId(oMeta), "_GradePdf.ps", sep = ""))
plot(oGradePdf)
dev.off()

postscript(paste(getTractId(oMeta), "_Pmf.ps", sep = ""), height = 3.5)
plot(oPmf)
dev.off()

postscript(paste(getTractId(oMeta), "_TotalTonPdf.ps", sep = ""), height = 3.5)
plot(oTotalTonPdf)
dev.off()

postscript(paste(getTractId(oMeta), "_TotalTonPdf_Marg.ps", sep = ""))
plotMarginals(oTotalTonPdf)
dev.off()

setwd(oldWd)