rely.R_tmpl

setwd("#<r_work>#");
source("include/util_plot.R");

strDataFile = "#<data_file>#";
strPlotFile = "#<plot_file>#";


boolEventHist = #<rely_event_hist>#;

listSeries1Val = #<series1_list>#;

listLegendUsr = #<legend>#;

listPlotDisp	= #<plot_disp>#;

listPlotOrderSeries	= #<order_series>#;

listColors = #<colors>#;

listPch = #<pch>#;

listType = #<type>#;

listLty = #<lty>#;

listLwd = #<lwd>#;

listPlotCI = #<plot_ci>#;

listStatisticsCurvs	= #<summary_curves>#;

addReferenceLine = #<add_reference_line>#;

addSkillLine = #<add_skill_line>#;

boolIndyStagger= #<indy1_stagger>#;

options(warn=1)

legend=c();

if(length(listSeries1Val) > 0){
  listSeriesVar = names(listSeries1Val);
  matPermVal = permute(listSeries1Val);
  intSeriesIndex=0;
  series = list();
}

# specify the output file
strPlotType = "#<plot_type>#";
if( "jpeg" == strPlotType ){
  jpeg(
    strPlotFile,			      # output file name, with appropriate extension
  		width=#<plot_width>#,	  # width of the display region, in specified units
  		height=#<plot_height>#,	# height of the display region, in specified units
  		res=#<plot_res>#,		    # resolution, in dots per inch
  		units="#<plot_units>#",	# units for height and width
  		pointsize = 12,
  		quality = 100
  );
}else{
  bitmap(
    strPlotFile,			# output file name, with appropriate extension
    type="#<plot_type>#",	# type of bitmap
    width=#<plot_width>#,	# width of the display region, in specified units
    height=#<plot_height>#,	# height of the display region, in specified units
    res=#<plot_res>#,		# resolution, in dots per inch
    units="#<plot_units>#",	# units for height and width
    taa=4,					# bits of anti-aliasing for graphics
    gaa=4
  );
}

# specify plot window formatting
par(
  mar=#<mar>#,			# lines of margin in form c(bottom, left, top, right)
	mgp=#<mgp>#,			# margin line for axis title, axis labels and axis line in form c(title, labels, line)
	cex=#<cex>#,			# magnification for plot text
	las=2 					# axis label orientation - 0: parallel to axis, 1: horiz, 2: perp to axis, 3: vert
);


initSeries <- function(dfPerm,envir, intPermVal=1) {
  # adjust the percentile thresholds to the midpoints
  if(nrow(dfPerm) > 0){
    dfPerm_obar = dfPerm[dfPerm$stat_name == 'PSTD_BASER',]
    dfPerm_calibration = dfPerm[dfPerm$stat_name == 'PSTD_CALIBRATION',]
    dfPerm_n_i = dfPerm[dfPerm$stat_name == 'PSTD_NI',]

    dfPerm = dfPerm_calibration;
    dfPerm$thresh_ii = dfPerm$thresh_i;
    dfPerm$n_i = dfPerm_n_i$stat_value
    for(i in 2:nrow(dfPerm)){
      dfPerm[i-1,]$thresh_i = 0.5*(dfPerm[i,]$thresh_ii + dfPerm[i-1,]$thresh_ii);
    }
    dfPerm[i,]$thresh_i = 0.5*(1 + dfPerm[i,]$thresh_ii);

    dfPerm$no_skill = 0.5*(dfPerm$thresh_i + dfPerm_obar$stat_value);
    
    # initialize the plot data
    listCal = dfPerm$stat_value;

    listThresh = dfPerm$thresh_i;
    listBcl = dfPerm$stat_bcl;
    listBcu = dfPerm$stat_bcu;
  }else{
    dfPerm=NA;
    dfPermHist=data.frame(thresh_i=envir$listX, stat_value=rep(0, length(envir$listX)) );
  }
  if( !envir$isPlotCreated & any( !is.na(dfPerm) ) ){
    # create a new plot with no axes
    plot(
      x=listThresh, # series x data
      y=listCal, # series y data
      type="n", # series line type - p: points, l: lines, o: overplotted points and lines,
      # b or c: points joined by lines, s or S: stair steps, h: histogram, n: nothing
      xaxt="n", # x-axis type - n: no axis, s: compatibility with l and t
      xaxs="r", # x-axis style - r: regular, i: internal
      yaxt="n", # y1-axis type - n: no axis, s: compatibility with l and t
      yaxs="r", # x-axis style - r: regular, i: internal
      xlim=c(0,1), # x-axis limits in form c(x1,x2); x1 > x2 reverses axis
      ylim=c(0,1), # y1-axis limits
      xlab="", # x-axis label
      ylab="", # y1-axis label
      fg="gray"
    );
    
    # specify the axis formatting
    par(
      usr=c(0, 1, 0, 1), # vector c(x1, x2, y1, y2) giving the extreme coordinates of the plotting region
      xpd=NA # plot clipping - FALSE: plot, TRUE: figure, NA: device
    );
 
  
    # plot title
    mtext(
      "#<plot_title>#",		# label
	    font=#<title_weight>#,	# 1: plain text, 2: bold, 3: italic, 4: bold italic and 5: symbol
	    las=0, 					# axis label orientation - 0: parallel to axis, 1: horiz, 2: perp to axis, 3: vert
	    side=3,					# axis of interest, 1: below, 2: left, 3: above and 4: right
	    cex=#<title_size>#,		# relative magnification
	    padj=#<title_offset>#,	# axis label perpendicular location adjustment
	    adj=#<title_align>#,	# axis label parallel location adjustment
	    line=1					# margin line, starting at 0 counting outwards
    );
  
    # x-axis label
    mtext(
      "#<x_label>#",			# label
	    font=#<xlab_weight>#,	# 1: plain text, 2: bold, 3: italic, 4: bold italic and 5: symbol
	    las=0, 					# axis label orientation - 0: parallel to axis, 1: horiz, 2: perp to axis, 3: vert
	    side=1,					# axis of interest, 1: below, 2: left, 3: above and 4: right
	    cex=#<xlab_size>#,		# relative axis label magnification
	    padj=#<xlab_offset>#,	# axis label perpendicular location adjustment
	    adj=#<xlab_align>#,		# axis label parallel location adjustment
	    line=1					# margin line, starting at 0 counting outwards
    );
  
    # specify the x-axis
    axis(
      1,						# axis of interest, 1: below, 2: left, 3: above and 4: right
      at=seq(0,1,.1),			# locations of tick marks
      tick=TRUE,				# tick marks on/off
      tck=.01, 				# length of tick marks as a fraction of plot window
      las=#<xtlab_orient>#,	# axis label orientation - 0: parallel to axis, 1: horiz, 2: perp to axis, 3: vert
      padj=#<xtlab_perp>#,	# axis label perpendicular location adjustment
      hadj=#<xtlab_horiz>#,	# axis label parallel location adjustment
      cex.axis=#<xtlab_size>#,# relative axis label magnification
      col="gray"				# axis color
    );
  
    # y-axis label
    mtext(
      "#<y1_label>#",			# label
      font=#<ylab_weight>#,	# 1: plain text, 2: bold, 3: italic, 4: bold italic and 5: symbol
      las=0, 					# axis label orientation - 0: parallel to axis, 1: horiz, 2: perp to axis, 3: vert
      side=2,					# axis of interest, 1: below, 2: left, 3: above and 4: right
      cex=#<ylab_size>#,		# relative axis label magnification
      padj=#<ylab_offset>#,	# axis label perpendicular location adjustment
      adj=#<ylab_align>#,		# axis label parallel location adjustment
      line=1					# margin line, starting at 0 counting outwards
    );
    # specify the y1-axis
    axis(
      2, # axis of interest, 1: below, 2: left, 3: above and 4: right
      at=seq(0,1,.1), # locations of tick marks
      labels=TRUE, # tick labels on/off
      tick=TRUE, # tick marks on/off
      tck=.01, # length of tick marks as a fraction of plot window
      las=#<ytlab_orient>#,	# axis label orientation - 0: parallel to axis, 1: horiz, 2: perp to axis, 3: vert
      padj=#<ytlab_perp>#,	# axis label perpendicular location adjustment
      hadj=#<ytlab_horiz>#,	# axis label parallel location adjustment
      cex.axis=#<ytlab_size>#,# relative axis label magnification
      col="gray" # axis color
    );


    # set the axis bounds for y1 plotting
    listY1Usr = par()$usr;
    par(new=TRUE, usr=listY1Usr, xaxs="r");
    envir$isPlotCreated=TRUE;
  }

  actualInd = listPlotOrderSeries[intPermVal];

  #if do not show or no data - continue to the next
  if( listPlotDisp[actualInd] & any( !is.na(dfPerm) ) ){
    # if the event histogram is activated, plot the number of events by threshold on the y2 axis
    if( envir$boolEventHist ){
      dfPermHist = subset(dfPerm, select=c(thresh_i, n_i));
  
      if( intPermVal == 1){
        envir$res = dfPermHist;
      }else{
        envir$res = merge(x=envir$res, y=dfPermHist, by="thresh_i", all=TRUE, suffixes=c("", paste(".",(envir$suffixesInd), sep="")));
        envir$suffixesInd=envir$suffixesInd+1;
      }
    }
  
    envir$o_bar[[actualInd]] = dfPerm_obar$stat_value[1];
    envir$listThresh[[actualInd]] = listThresh;
    envir$listCal[[actualInd]] = listCal;
    envir$listBcl[[actualInd]] = listBcl;
    envir$listBcu[[actualInd]] = listBcu;
    dblNoSkillM = (dfPerm[2,]$no_skill - dfPerm[1,]$no_skill) / (dfPerm[2,]$thresh_i - dfPerm[1,]$thresh_i);
    dblNoSkillB = dfPerm[1,]$no_skill - dblNoSkillM * dfPerm[1,]$thresh_i;
    envir$dblNoSkillB[[actualInd]] = dblNoSkillB;
    envir$dblNoSkillM[[actualInd]] = dblNoSkillM;
 
  }

}


dfPct = read.delim(strDataFile,header=TRUE,sep='');
env <- new.env(parent=emptyenv())
env$isPlotCreated=FALSE;
env$boolEventHist = boolEventHist;
env$suffixesInd = 1;
env$listX = unique(dfPct$thresh_i);
env$listThresh=list();
env$listCal=list();
env$listBcl=list();
env$listBcu=list();
env$dblNoSkillB=list();
env$o_bar=list();

intCurveNum = 0;

if(length(listSeries1Val) > 0){
  for(intPermVal in 1:nrow(matPermVal)){
    seriesLegend="";
    listPermVal = matPermVal[intPermVal,];
    dfPerm = dfPct;
    for(intVar in 1:length(listSeriesVar)){
      valPerm = listPermVal[intVar];
      if( grepl("^[0-9]+$", valPerm) ){ valPerm = as.integer(valPerm); }
      dfPerm = dfPerm[dfPerm[[ listSeriesVar[intVar] ]] == valPerm,];
      seriesLegend = paste(seriesLegend,valPerm);
    }
    
    seriesLegend = paste(seriesLegend,'Reliability Curve');
    legend=append(legend, seriesLegend);
    initSeries(dfPerm,envir=env, intPermVal);
    intCurveNum = intPermVal;
  }
}else{
  dfPerm = dfPct;
  legend=append(legend, 'Reliability Curve');
  initSeries(dfPerm,envir=env);
  intCurveNum = 1;
}

if(length(listStatisticsCurvs) > 0){
  intPermVal = intCurveNum + 1;
  for(index in 1:length(listStatisticsCurvs)){
    groupStats = aggregate(stat_value ~ thresh_i, data = dfPct[dfPct$stat_name == 'PSTD_CALIBRATION',], FUN = listStatisticsCurvs[index]);
    groupStats$thresh_ii = groupStats$thresh_i;
    for(i in 2:nrow(groupStats)){
      groupStats[i-1,]$thresh_i = 0.5*(groupStats[i,]$thresh_ii + groupStats[i-1,]$thresh_ii);
    }
    groupStats[i,]$thresh_i = 0.5*(1 + groupStats[i,]$thresh_ii);
    seriesLegend = paste(listStatisticsCurvs[index],'Reliability Curve');
    legend=append(legend, seriesLegend);

    env$listThresh[[length(env$listThresh)+1]] = groupStats$thresh_i;
    env$listCal[[length(env$listCal)+1]] = groupStats$stat_value;
  }
}

#replace with user's legends
if( 0 < length(listLegendUsr) ){
  for(index in 1:length(listLegendUsr)){
    if(nchar(listLegendUsr[index]) > 0){
      legend[index]=listLegendUsr[index];
    }
  }
}

listColorsOrdered=c();
listLegendOrdered=c();
listLtyOrdered=c();
listLwdOrdered=c();
listCiOrdered=c();
listPchOrdered=c();
listTypeOrdered=c();
listPlotDispOrdered=c();
for(index in 1:length(listPlotOrderSeries)){
  listLwdOrdered[listPlotOrderSeries[index]]=listLwd[index];
  listCiOrdered[listPlotOrderSeries[index]]=listPlotCI[index];
  listColorsOrdered[listPlotOrderSeries[index]]=listColors[index];
  listLegendOrdered[listPlotOrderSeries[index]]=legend[index];
  listLtyOrdered[listPlotOrderSeries[index]]=listLty[index];
  listPchOrdered[listPlotOrderSeries[index]]=listPch[index];
  listTypeOrdered[listPlotOrderSeries[index]]=listType[index];
  listPlotDispOrdered[listPlotOrderSeries[index]]=listPlotDisp[index];
}

listLwdLeg = listLwdOrdered;
listLwdLeg[listLwdLeg > 3] = 3;

# mask the hidden elements from the legend and color list
listLegendDisp = listLegendOrdered[listPlotDispOrdered];
listColorsDisp = listColorsOrdered[listPlotDispOrdered];
listLtyDisp    = listLtyOrdered[listPlotDispOrdered];
listPchDisp    = listPchOrdered[listPlotDispOrdered];
listTypeDisp    = listTypeOrdered[listPlotDispOrdered];



if(boolEventHist == TRUE){
  listHistData=vector();
  listVisibleLegend=vector();
  listVisibleColors=vector();
  for(n in 2:length(names(env$res))){
    actualInd = listPlotOrderSeries[n-1];
    if(listPlotDisp[actualInd]){
      listHistData=rbind(listHistData,env$res[[actualInd+1]]);
      listVisibleLegend=append(listVisibleLegend, legend[actualInd]);
      listVisibleColors=append(listVisibleColors, listColors[actualInd]);
    }
  }
  # y2 axis label
  mtext(
    "#<y2_label>#",			# label
		font=#<y2lab_weight>#,	# 1: plain text, 2: bold, 3: italic, 4: bold italic and 5: symbol
		las=0, 					# axis label orientation - 0: parallel to axis, 1: horiz, 2: perp to axis, 3: vert
		side=4,					# axis of interest, 1: below, 2: left, 3: above and 4: right
		cex=#<y2lab_size>#,		# relative axis label magnification
		padj=#<y2lab_offset>#,	# axis label perpendicular location adjustment
		adj=#<y2lab_align>#,	# axis label parallel location adjustment
		line=1					# margin line, starting at 0 counting outwards
  );
  # determine and set the limits for the y2 axis
  resTemp = env$res;
  resTemp$thresh_i = NULL;
  listY2Tick = pretty(range(resTemp,na.rm = TRUE), 4);
  listY2Usr = c(par()$usr[1:2], listY2Tick[1], listY2Tick[length(listY2Tick)]);
  par(new=TRUE,usr=listY2Usr, xaxs="r");
  
 
  # build the bar plot
  barplot(
    listHistData,			# bar data 
    width = 1,				# vector of bar widths, recycled if necessary
    space = NULL,			# space left before each bar
    #legend.text = listVisibleLegend, 	# vector of text used to construct a legend
    beside = TRUE,			# logical value indicating stacked or side-by-side bars - TRUE: side-by-side, FALSE: stacked
    horiz = FALSE,			# logical value indicating horizontal or vertical bars
    col = listVisibleColors,		# vector of box colors
    border = par("fg"),		# box border
    #tcl=.1,					# tick mark length
    #ylim=#<y1_lim>#,		# y axis limits
    axes=FALSE    
  );
  # add the y2 axis
  axis(
    4, # axis of interest, 1: below, 2: left, 3: above and 4: right
    labels=TRUE, # tick labels on/off
    tick=TRUE, # tick marks on/off
    tck=.01, # length of tick marks as a fraction of plot window
    las=#<y2tlab_orient>#,	# axis label orientation - 0: parallel to axis, 1: horiz, 2: perp to axis, 3: vert
	  padj=#<y2tlab_perp>#,	# axis label perpendicular location adjustment
	  hadj=#<y2tlab_horiz>#,	# axis label parallel location adjustment
	  cex.axis=#<y2tlab_size>#,# relative axis label magnification
    col="gray" # axis color
  );
}


# specify the axis formatting
par(
  usr=c(0, 1, 0, 1), # vector c(x1, x2, y1, y2) giving the extreme coordinates of the plotting region
  xpd=NA # plot clipping - FALSE: plot, TRUE: figure, NA: device
);

if(addSkillLine){
  # draw the max-skill line
  abline(
    a=0, # intercept
    b=1, # slope
    lty=1, # line type - 0: blank, 1: solid, 2: dashed, 3: dotted, 4: dotdash, 5: longdash, 6: twodash
    col="gray", # line color
    lend=1, # line end style - 0: rounded, 1: butt, 2: square
    cex=1, # line magnification
    lwd=1.5, # line width
    xpd=FALSE # plot clipping - FALSE: plot, TRUE: figure, NA: device
  );
}

if(addReferenceLine){
  # draw the no-resolution line
  for (i in 1:length(env$o_bar)){
    if(!is.na(env$o_bar[i][[1]])){
      abline(
        a=env$o_bar[i][[1]], # intercept
        b=0, # slope
        lty=2, # line type - 0: blank, 1: solid, 2: dashed, 3: dotted, 4: dotdash, 5: longdash, 6: twodash
        col=listColorsDisp[i], # line color
        lend=1, # line end style - 0: rounded, 1: butt, 2: square
        cex=1, # line magnification
        lwd=1.5, # line width
        xpd=FALSE # plot clipping - FALSE: plot, TRUE: figure, NA: device
      );
    }else{
      cat(" WARNING: reference  line can't be created for the " , i , "series \n");
    }
  }
}

# draw the no-skill line
for (i in 1:length(env$dblNoSkillB)){
  if(!is.na(env$dblNoSkillB[i][[1]]) && !is.na(env$dblNoSkillM[i][[1]]) ){
    abline(
      a=env$dblNoSkillB[i][[1]], # intercept
      b=env$dblNoSkillM[i][[1]], # slope
      lty=2, # line type - 0: blank, 1: solid, 2: dashed, 3: dotted, 4: dotdash, 5: longdash, 6: twodash
      col=listColorsDisp[i], # line color
      lend=1, # line end style - 0: rounded, 1: butt, 2: square
      cex=1, # line magnification
      lwd=1.5, # line width
      xpd=FALSE # plot clipping - FALSE: plot, TRUE: figure, NA: device
    );
  }else{
      cat("  WARNING: no-skill line can't be created for the " , i , "series \n");
  }
}

# build the x-adjustment increments, if requested
dblAdjScale = (tail(env$listThresh[[length(env$listThresh)]], n=1) - head(env$listThresh[[length(env$listThresh)]], n=1)) / 150;
# build the x-adjustment increments, if requested
intNumStag = length(env$listThresh);
listX1Stag = rep(0, intNumStag + 1);
if( TRUE == boolIndyStagger ){
  listX1Stag = seq( -(intNumStag/2)*dblAdjScale, (intNumStag/2)*dblAdjScale, length=intNumStag + 1 ) + ( dblAdjScale/2 );
}


# remove the NAs from the x and y lists
xValues  = env$listThresh[i][[1]]+ listX1Stag[i];
yValues  = env$listCal[i][[1]];
listDataInd  = !( is.na(env$listCal[i][[1]]) );
xValues  = xValues[listDataInd];
yValues  = yValues[listDataInd];


# plot the series line
for (i in 1:length(env$listThresh)){

  # remove the NAs from the x and y lists
  xValues  = env$listThresh[i][[1]]+ listX1Stag[i];
  yValues  = env$listCal[i][[1]];
  listDataInd  = !( is.na(env$listCal[i][[1]]) );
  xValues  = xValues[listDataInd];
  yValues  = yValues[listDataInd];

  lines(
    x=xValues, # x coordinates
    y=yValues, # y coordinates
    type=listTypeDisp[i], # series line type - p: points, l: lines, o: overplotted points and lines,
    # b or c: points joined by lines, s or S: stair steps, h: histogram, n: nothing
    pch=listPchDisp[i], # point symbol - see example(points)
    font=5, # 1: plain text, 2: bold, 3: italic, 4: bold italic and 5: symbol
    lty=listLtyDisp[i], # line type - 0: blank, 1: solid, 2: dashed, 3: dotted, 4: dotdash, 5:
# longdash, 6: twodash
    col=listColorsDisp[i], # line color
    lend=1, # line end style - 0: rounded, 1: butt, 2: square
    cex=1, # line magnification
    lwd=listLwdLeg[i] # line width
  );
  # add arrows to the line, if requested
  if( "none" != listPlotCI[i] ) {
    listXFr  = env$listThresh[i][[1]] + listX1Stag[i];
    listYFr  = env$listBcl[i][[1]];
    listYTo  = env$listBcu[i][[1]];

    listXTo  = listXFr;

    for(ind in 1: length( listXFr) ){

      arrows(
        listXFr[ind],                      # x coordinate from which to draw
        listYFr[ind],                      # y coordinate from which to draw
        listXTo[ind],                      # x coordinate to which to draw
        listYTo[ind],                      # y coordinate to which to draw
        length=.02,                    # length of the edges of the arrow head (in inches)
        angle=90,                      # angle from the shaft of the arrow to the edge of the arrow head
        code=3,                        # arrow type - 1: up, 2: down and 3: both
        lwd=listLwd[i],              # line width
        col=listColorsDisp[i],      # line color
        ljoin=1,                      # line join style - 0: round, 1: mitre and 2: bevel
        lmitre=1,                      # line mitre limit, range 1-10
        lend=2                        # line end style - 0: round, 1: butt and 2: square
      );
    }
  }

}


# add the legend
legend(
  x="bottom", # x coordinate of legend
  listLegendDisp, # vector of legend labels
  col=listColorsDisp, # vector of point and line colors
  lty=listLtyDisp, # vector of line types
  lwd=listLwdLeg, # vector of line widths
  pch=listPchDisp, # vector of point symbols
  merge=FALSE, # merge points and lines
  cex=#<legend_size>#, # relative magnification
  bty="#<legend_box>#", # legend box type - o: box, n: none
  adj=0, # string adjustment for text
  xpd=TRUE, # clipping - FALSE: plot region, TRUE: figure region, NA: device region
  ncol=#<legend_ncol>#, # number of columns in legend
  inset=#<legend_inset>#, # inset distance from margin
  x.intersp=1, # horizontal character spacing factor
  y.intersp=.8 # vertical character spacing factor
);


# plot caption
mtext(
"#<plot_caption>#",				# caption
	font=#<caption_weight>#,		# 1: plain text, 2: bold, 3: italic, 4: bold italic and 5: symbol
	col="#<caption_col>#",			# font color
	las=0, 							# axis label orientation - 0: parallel to axis, 1: horiz, 2: perp to axis, 3: vert
	side=1,							# axis of interest, 1: below, 2: left, 3: above and 4: right
	cex=#<caption_size>#,			# relative magnification
	padj=#<caption_offset>#,		# axis perpendicular location adjustment
	adj=#<caption_align>#,			# axis parallel location adjustment
	line=1	
);  

#<plot_cmd>#

invisible(dev.off());