########################################################################
##
##   Name: mode_summary.R
##
##   Description:
##      This R script summarizes the contents of one or more MODE
##      object statistics file.
##
##   Usage:
##      Rscript mode_summary.R mode_obj_file_list
##
##   Arguments:
##      "mode_obj_file_list" is a space-seperated list of the MODE
##      object statistics file names to be used.
##
##   Details:
##      This R script summarizes the contents of one or more MODE
##      object statistics files.  It reads in the MODE file(s) passed
##      on the command line, counts the numbers/areas of the
##      matched/unmatched forecast/observation objects and prints them
##      to the screen.  It also computes the maximum interest value for
##      each object and prints out the median of them for the forecast
##      objects, observation objects, and both.
##
##      This R script is meant as an example.  Users are welcome and
##      encouraged to adapt it to perform the type of analysis they
##      need.
##
##   Examples:
##      Rscript mode_summary.R METv1.1/out/mode/mode_APCP_12_SFC_vs_APCP_12_SFC_120000L_20050807_120000V_120000A_obj.txt
##
##   Author:
##      John Halley Gotway, MET Development Team, 10/22/2008
##      johnhg@ucar.edu
##
########################################################################

########################################################################
#
# Handle the arguments.
#
########################################################################

# Retreive the arguments
args <- commandArgs(TRUE)

# Check the number of arguments
if(length(args) < 1) {
   cat("Usage: mode_summary.R mode_obj_file_list\n")
   cat("   where mode_obj_file_list is a space-separated list of MODE object statistic files.\n")
   quit()
}

# Constants
i_start<-1

# Variables to store the object counts
sum_n_fcst<-0
sum_n_fcst_match<-0
sum_n_fcst_unmatch<-0
sum_n_obs<-0
sum_n_obs_match<-0
sum_n_obs_unmatch<-0
sum_n_fcst_comp<-0
sum_n_obs_comp<-0

# Variables to store the object areas
sum_a_fcst<-0
sum_a_fcst_match<-0
sum_a_fcst_unmatch<-0
sum_a_obs<-0
sum_a_obs_match<-0
sum_a_obs_unmatch<-0
sum_a_fcst_comp<-0
sum_a_obs_comp<-0

# Handle each of the MODE files one at a time
for(i in i_start:length(args)) {

   # Store the current MODE file
   mode_file<-args[i]

   # Read the MODE file passed as an argument
   mode<-read.table(mode_file, header=TRUE)

   cat(paste("Read", dim(mode)[1], "lines from MODE file:", mode_file, "\n"))

   # Check for an empty file containing only the header row
   if(dim(mode)[1] > 0) {

      # Replace instances of -9999 with NA
      mode[mode==-9999]=NA

      ########################################################################
      #
      # Compute indicators for each of the line types.
      #
      ########################################################################

      ind_simp_fcst <-regexpr("^F[0-9][0-9][0-9]$", mode$OBJECT_ID) == 1
      ind_simp_obs  <-regexpr("^O[0-9][0-9][0-9]$", mode$OBJECT_ID) == 1
      ind_simp_pair <-regexpr("^F[0-9][0-9][0-9]_O[0-9][0-9][0-9]$", mode$OBJECT_ID) == 1

      ind_comp_fcst <-regexpr("^CF[0-9][0-9][0-9]$", mode$OBJECT_ID) == 1
      ind_comp_obs  <-regexpr("^CO[0-9][0-9][0-9]$", mode$OBJECT_ID) == 1
      ind_comp_pair <-regexpr("^CF[0-9][0-9][0-9]_CO[0-9][0-9][0-9]$", mode$OBJECT_ID) == 1

      ind_match     <-(regexpr("^C[F,O]000$", mode$OBJECT_CAT) == -1) *
                      (regexpr("^C[F,O][0-9][0-9][0-9]$", mode$OBJECT_CAT) == 1)
      ind_unmatch   <-regexpr("^C[F,O]000$", mode$OBJECT_CAT) == 1

      ########################################################################
      #
      # Keep a running sum of object counts and areas.
      #
      ########################################################################

      # Keep a running sum of the object counts
      sum_n_fcst         <-sum_n_fcst+sum(ind_simp_fcst)
      sum_n_fcst_match   <-sum_n_fcst_match+sum(ind_simp_fcst*ind_match > 0)
      sum_n_fcst_unmatch <-sum_n_fcst_unmatch+sum(ind_simp_fcst*ind_unmatch > 0)

      sum_n_obs          <-sum_n_obs+sum(ind_simp_obs)
      sum_n_obs_match    <-sum_n_obs_match+sum(ind_simp_obs*ind_match > 0)
      sum_n_obs_unmatch  <-sum_n_obs_unmatch+sum(ind_simp_obs*ind_unmatch > 0)

      sum_n_fcst_comp    <-sum_n_fcst_comp+sum(ind_comp_fcst)
      sum_n_obs_comp     <-sum_n_obs_comp+sum(ind_comp_obs)

      # Keep a running sum of the object areas
      sum_a_fcst         <-sum_a_fcst+sum(mode$AREA[ind_simp_fcst])
      sum_a_fcst_match   <-sum_a_fcst_match+sum(mode$AREA[ind_simp_fcst*ind_match > 0])
      sum_a_fcst_unmatch <-sum_a_fcst_unmatch+sum(mode$AREA[ind_simp_fcst*ind_unmatch > 0])

      sum_a_obs          <-sum_a_obs+sum(mode$AREA[ind_simp_obs])
      sum_a_obs_match    <-sum_a_obs_match+sum(mode$AREA[ind_simp_obs*ind_match > 0])
      sum_a_obs_unmatch  <-sum_a_obs_unmatch+sum(mode$AREA[ind_simp_obs*ind_unmatch > 0])

      sum_a_fcst_comp    <-sum_a_fcst_comp+sum(mode$AREA[ind_comp_fcst])
      sum_a_obs_comp     <-sum_a_obs_comp+sum(mode$AREA[ind_comp_obs])

      ########################################################################
      #
      # Keep track of maximum interest values for each object.
      #
      ########################################################################

      # Use only the simple pair lines.
      ind<-ind_simp_pair

      # Find the maximum interest for each forecast object.
      mif <-aggregate(mode$INTEREST[ind], by=list(obj=substr(mode$OBJECT_ID[ind], 1, 4)), FUN=max)

      # Store the area for each of the forecast objects.
      ind_area <-charmatch(mif$obj, mode$OBJECT_ID)
      mif$area <-mode$AREA[ind_area]

      # Find the maximum interest for each observation object.
      mio <-aggregate(mode$INTEREST[ind], by=list(obj=substr(mode$OBJECT_ID[ind], 6, 9)), FUN=max)

      # Store the area for each of the observation objects.
      ind_area <-charmatch(mio$obj, mode$OBJECT_ID)
      mio$area <-mode$AREA[ind_area]

      # Append this array to the previous ones
      if(exists("list_mif")) list_mif <-rbind(list_mif, mif)
      else                   list_mif <-mif

      if(exists("list_mio")) list_mio  <-rbind(list_mio, mio)
      else                   list_mio  <-mio
   } # end if statement

} # end for loop

# Combine the maximum interests by fcst and obs objects.
list_mia <-rbind(list_mif, list_mio)

########################################################################
#
# Print the results to the screen.
#
########################################################################

# Print out the counts and areas of the single objects
cat(paste("Total Number of Files Processed =",
          length(args), "\n"))
cat(paste("Total Number of Single Objects =",
          sum_n_fcst+sum_n_obs, "\n"))

cat(paste("Number of Single Fcst Objects =",
          sum_n_fcst, "\n"))
cat(paste("Number of Matched Single Fcst Objects =",
          sum_n_fcst_match, "\n"))
cat(paste("Number of Unmatched Single Fcst Objects =",
          sum_n_fcst_unmatch, "\n"))

cat(paste("Number of Single Obs Objects =",
          sum_n_obs, "\n"))
cat(paste("Number of Matched Single Obs Objects =",
          sum_n_obs_match, "\n"))
cat(paste("Number of Unmatched Single Obs Objects =",
          sum_n_obs_unmatch, "\n"))

cat(paste("Total Number of Composite Objects =",
          sum_n_fcst_comp+sum_n_obs_comp, "\n"))
cat(paste("Total Number of Fcst Composite Objects =",
          sum_n_fcst_comp, "\n"))
cat(paste("Total Number of Obs Composite Objects =",
          sum_n_obs_comp, "\n"))

cat(paste("Total Area of Objects =",
          sum_a_fcst+sum_a_obs, "\n"))

cat(paste("Area of Single Fcst Objects =",
          sum_a_fcst, "\n"))
cat(paste("Area of Matched Single Fcst Objects =",
          sum_a_fcst_match, "\n"))
cat(paste("Area of Unmatched Single Fcst Objects =",
          sum_a_fcst_unmatch, "\n"))

cat(paste("Area of Single Obs Objects =",
          sum_a_obs, "\n"))
cat(paste("Area of Matched Single Obs Objects =",
          sum_a_obs_match, "\n"))
cat(paste("Area of Unmatched Single Obs Objects =",
          sum_a_obs_unmatch, "\n"))

# Print interest measures for fcst objects
cat(paste("Median of Max Interest by Fcst Object =",
          median(list_mif$x), "\n"))
cat(paste("Mean of Max Interest by Fcst Object =",
          mean(list_mif$x), "\n"))
cat(paste("Area-Weighted Mean of Max Interest by Fcst Object =",
          weighted.mean(list_mif$x, list_mif$area), "\n"))

# Print interest measures for obs objects
cat(paste("Median of Max Interest by Obs Object =",
          median(list_mio$x), "\n"))
cat(paste("Mean of Max Interest by Obs Object =",
          mean(list_mio$x), "\n"))
cat(paste("Area-Weighted Mean of Max Interest by Obs Object =",
          weighted.mean(list_mio$x, list_mio$area), "\n"))

# Print interest measures for fcst + obs objects
cat(paste("Median of Max Interest by Fcst+Obs Objects =",
          median(list_mia$x), "\n"))
cat(paste("Mean of Max Interest by Fcst+Obs Object =",
          mean(list_mia$x), "\n"))
cat(paste("Area-Weighted Mean of Max Interest by Fcst+Obs Object =",
          weighted.mean(list_mia$x, list_mia$area), "\n"))

# Optionally, save the output
# save.image()