############################################################################### 
# Master Recession Curve (MRC) Plotting Program MRCplot
# John Nimmo
# Started 2/9/17
# This R program plots parametric Master Recession Curves (MRCs).
############################################################################### 

#### User inputs ############################################################## 

# Folders for input data and output:
inputlocation = "P:/current_hydrograph-sites/SH_CZO/outputs/"
outputlocation = "P:/current_hydrograph-sites/SH_CZO/outputs/"

# File names for fitted MRC parameters and points, for fits 1, 2, 3, . . . 
inputfiles = c("SH-Q.2009.7mo.MRCfit1a.out.txt", 
               "SH-Q.2009.7mo.MRCfit2a.out.txt") 
inputpaths = as.vector(paste0(inputlocation, inputfiles))
nfits = length(inputfiles)

# If no output file is desired, use an empty string ("") for the file name
outputtextfile = "MRC.SH-Q.5.txt"

# Graphics specifications ###############

# Axis limits (set = for autogenerated limits)
xlimsl = c(0, 0)
ylimsl =  c(0, 0)

# Put specifications on MRC plot?)
annotate_specs = TRUE

#### End user inputs ########################################################## 
 
timerun = substr(Sys.time(), 1, 19) 
outfilepath = paste0(outputlocation, outputtextfile)

# coefs is a matrix whose columns contain the parameters of a fitted segment.
# Each row of this MRC coefficient matrix represents a different part (portion 
# of range) of the MRC, with the intercept in the first column and the slope in 
# the second.
coefs = matrix(nrow=nfits, ncol=2) 

# Read fitted MRC parameters and points 

linesf = lapply(inputpaths, readLines) 
timecalc = vapply(1:nfits, function(iarg) {paste0("Fit ", iarg, ":  ", 
                      linesf[[iarg]][3])}, character(1)) 
allpts = list() 
for (i1 in 1:nfits) {
  icoef = which(substring(linesf[[i1]], 1, 8)=="mrc_coef") 
  strcoef = substr(linesf[[i1]][icoef], 11, 80) 
  coefs[i1,] = eval(parse(text=paste0("coefs[i1,] =", strcoef))) 
  allpts[[i1]] = read.table(inputpaths[i1], header=TRUE, skip=icoef+1)
}
# Find intersection points to use in evaluating MRC function: 
h_intersect = vector()
s_intersect = vector()
for (i1 in 1:(nfits-1)) { 
  i1p1 = i1 + 1
  h_intersect[i1] = (coefs[i1,1]-coefs[i1p1,1]) / (coefs[i1p1,2]-coefs[i1,2]) 
  s_intersect[i1] = coefs[i1,1] + coefs[i1,2] * h_intersect 
}
# Define MRC function: 
MRCfn = function(rsp) {           # rsp is dummy variable for response 
  ipart = vector(length = length(rsp)) 
  hintvec = c(-1e20, h_intersect)
  for (i in 1:length(rsp)) {ipart[i] = max(which(hintvec < rsp[i]))}
  cat(rsp, "\n", ipart, "\n") 
  MRCslope = coefs[ipart,1] + coefs[ipart,2] * rsp 
  return(MRCslope)
}
allptsslope = allpts[[1]][,2]
allptsslope = c(allptsslope, allpts[[2]][,2]) ## change later to allow nfits>3
allptsresp = allpts[[1]][,1]
allptsresp = c(allptsresp, allpts[[2]][,1]) ## change later to allow nfits>3 

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

# Show results on screen and save to file
dataname = substr(inputfiles[1], 1, 8)    ## change later, to be more general 
title = paste0("Master Recession Curve for ", dataname)
linesout = vector()
linesout[1:2] = c(paste("***** ", title, " *****"), 
          paste0("\n\nPlotted at ", timerun, "\n"))
for (i1 in 1:nfits) {linesout[i1+2] = paste0("coefs[", i1, "] = c(", 
                                      coefs[i1,1], ", ", coefs[i1,2], ")\n")}
linesout[3+nfits] = "\nMRC intersection point (H-value slope-value):\n"
for (i1 in 1:(nfits-1)) {
  linesout[i1+3+nfits] = paste("         ", h_intersect[i1], s_intersect[i1])} 
nl = length(linesout) + 1
linesout[nl] = "\n\n" 
for (i in 1:nfits) linesout[i+nl] = paste0(timecalc[i], "\n") 
cat("\n", linesout, "\n")
cat(file=outfilepath, linesout, "\n")

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

# Plot computed points and fitted MRC:			
dev.new(width=15,height=10)
if(xlimsl[1] == xlimsl[2]) {
  xlimsl = range(allptsslope) 
  leftedge = xlimsl[1]*1.03   ## unneeded?
} else {
  leftedge = xlimsl[1]
}
rightedge = xlimsl[2]
if(ylimsl[1] == ylimsl[2]) {ylimsl = range(allptsresp)} 
mrcverticesy = c(ylimsl[1], h_intersect[1:(nfits-1)], ylimsl[2])
mrcverticesx = MRCfn(mrcverticesy)
plot(mrcverticesx, mrcverticesy, type="l", xlim=xlimsl, ylim=ylimsl, xaxs="r",
     xlab="dH/dt", ylab="Hydrograph Response H", pch=5, 
     cex=1., col="black",	lwd=3, main=title) 

for (i1 in 1:nfits) {points(allpts[[i1]][,2], allpts[[i1]][,1], pch=16, 
                            cex=1., col=i1+1)} 
legend("bottomleft", "Fitted MRC", col="black", lty=1, lwd=3, bty="n")
if(annotate_specs) {for (i in seq(nfits, 1, -1)) {
  mtext(timecalc[i], cex=.8, side=3, adj=1, at=rightedge, line=.9*(i-nfits-1))} 
}
###############################################################################