docs: def multicategory_scatterplot(output_directory, file_prefix, df,
                                x_series_index, y_series_index, category_series_index,
                                series_color, plot_title = '',
                                x_axis_label = '', y_axis_label = '',
                                min_predicted_ddg = None, max_predicted_ddg = None, min_experimental_ddg = None, max_experimental_ddg = None):
      '''This function was adapted from the covariation benchmark.'''
  
      # todo: Abstract this graph from the current usage (DDG measurements).
      # todo: make the capped value for unquantified but classified measurements (e.g. DDG > 7 kcal/mol) parameterizable
      # todo: add an option to identify outliers by standard deviations (over the set of errors |x - y|) rather than by fixed value
      # todo: add an option to use geom_text_repel to avoid/reduce overlapping text
      # todo: allow users to provide colors for the facets / categories
  
      # Changeset
      # todo: Change it to take in a pandas dataframe instead of the data_table_headers + data_table parameters.
      # todo: Add exception if number of cases > 2 so the general case can be implemented once we have test data.
      # todo: use one column as the category e.g. "PDB". assert that there is a maximum number of categories. Test with > 2 categories
      # todo: remove all references to SNX27 and NHERF1 below and loop over the set of categories instead
  
      #print(df[facet_index])
      color_map = {}
      categories = list(df.ix[:, category_series_index].unique())
      print(type(categories))
      num_categories = len(categories)
      category_colors = get_spaced_plot_colors(num_categories)
      for x in xrange(num_categories):
          color_map[categories[x]] = '#' + category_colors[x]
  
      df['CategorizationColor'] = df.apply(lambda r: color_map[r[category_series_index]], axis = 1)
      categorization_color_index = len(df.columns.values) - 1
  
      # Monday: continue here
      print(df)
      sys.exit(0)
      try: os.mkdir(output_directory)
      except: pass
      assert(os.path.exists(output_directory))
  
  
      df['Categorization'] = df.apply(lambda r: _determine_fraction_correct_class(r[x_series_index], r[y_series_index])[0], axis = 1)
      categorization_index = len(df.columns.values) - 1
      df['CategorizationShape'] = df.apply(lambda r: _determine_fraction_correct_class(r[x_series_index], r[y_series_index])[1], axis = 1)
      categorization_shape_index = len(df.columns.values) - 1
  
  
  
  
      # Create the R script
      boxplot_r_script = '''
  library(ggplot2)
  library(gridExtra)
  library(scales)
  library(qualV)
  
  # PNG generation
  png('%(file_prefix)s.png', width=2560, height=2048, bg="white", res=600)
  txtalpha <- 0.6
  redtxtalpha <- 0.6
  
  %(png_plot_commands)s
      '''
  
      xy_table_filename = '{0}.txt'.format(file_prefix)
      xy_table_filepath = os.path.join(output_directory, xy_table_filename)
      write_file(xy_table_filepath, '\n'.join(','.join(map(str, line)) for line in [data_table_headers] + data_table))
  
      single_plot_commands = '''
  # Set the margins
  par(mar=c(5, 5, 1, 1))
  
  xy_data <- read.csv('%(xy_table_filename)s', header=T)
  
  names(xy_data)[%(x_series_index)d + 1] <- "xvalues"
  names(xy_data)[%(y_series_index)d + 1] <- "yvalues"
  
  # coefs contains two values: (Intercept) and yvalues
  coefs <- coef(lm(xvalues~yvalues, data = xy_data))
  fitcoefs = coef(lm(xvalues~0 + yvalues, data = xy_data))
  fitlmv_yvalues <- as.numeric(fitcoefs[1])
  lmv_intercept <- as.numeric(coefs[1])
  lmv_yvalues <- as.numeric(coefs[2])
  lm(xy_data$yvalues~xy_data$xvalues)
  
  xlabel <- "%(x_axis_label)s"
  ylabel <- "%(y_axis_label)s"
  plot_title <- "%(plot_title)s"
  rvalue <- cor(xy_data$yvalues, xy_data$xvalues)
  rvalue
  xy_data
  
  #3QDO = SNX27
  #1G9O = NHERF1
  
  valid_xy_data <- xy_data[which(xy_data$xvalues < 6.99),]
  rvalue <- cor(valid_xy_data$yvalues, valid_xy_data$xvalues)
  rvalue
  valid_xy_data
  
  valid_xy_data_NHERF1 <- xy_data[which(xy_data$xvalues < 6.99 & xy_data$PDB == '1G9O'),]
  rvalue_NHERF1 <- cor(valid_xy_data_NHERF1$yvalues, valid_xy_data_NHERF1$xvalues)
  rvalue_NHERF1
  valid_xy_data_NHERF1
  
  coefs_NHERF1 <- coef(lm(xvalues~yvalues, data = valid_xy_data_NHERF1))
  lmv_intercept_NHERF1 <- as.numeric(coefs_NHERF1[1])
  lmv_yvalues_NHERF1 <- as.numeric(coefs_NHERF1[2])
  
  valid_xy_data_SNX27 <- xy_data[which(xy_data$xvalues < 6.99 & xy_data$PDB == '3QDO'),]
  rvalue_SNX27 <- cor(valid_xy_data_SNX27$yvalues, valid_xy_data_SNX27$xvalues)
  rvalue_SNX27
  valid_xy_data_SNX27
  
  coefs_SNX27 <- coef(lm(xvalues~yvalues, data = valid_xy_data_SNX27))
  lmv_intercept_SNX27 <- as.numeric(coefs_SNX27[1])
  lmv_yvalues_SNX27 <- as.numeric(coefs_SNX27[2])
  
  lmv_intercept
  lmv_yvalues
  lmv_intercept_NHERF1
  lmv_yvalues_NHERF1
  lmv_intercept_SNX27
  lmv_yvalues_SNX27
  
  # Set graph limits and the position for the correlation value
  
  minx <- min(0.0, min(xy_data$xvalues) - 0.1)
  miny <- min(0.0, min(xy_data$yvalues) - 0.1)
  maxx <- max(1.0, max(xy_data$xvalues) + 0.1)
  maxy <- max(1.0, max(xy_data$yvalues) + 0.1)
      '''
  
      if min_predicted_ddg != None:
          single_plot_commands += '''
  miny <- min(miny  - 0.2, %(min_predicted_ddg)f  - 0.2)
      '''
      if max_predicted_ddg != None:
          single_plot_commands += '''
  maxy <- max(maxy + 0.5, %(max_predicted_ddg)f  + 0.5)
  
  miny <- -6
  maxy <- 12.5
  
      '''
      if min_experimental_ddg != None:
          single_plot_commands += '''
      minx <- min(minx, %(min_experimental_ddg)f)
          '''
      if max_experimental_ddg != None:
          single_plot_commands += '''
      maxx <- max(maxx, %(max_experimental_ddg)f) + 0.2
          '''
  
      single_plot_commands += '''
  xpos <- minx + 0.2
  ypos <- maxy - 1
  ypos_SNX27 <- ypos - 1
  ypos_NHERF1 <- ypos_SNX27 - 1
  
  lrt <- expression('R'^tst)
  
  p <- qplot(main="", xvalues, yvalues, data=xy_data, xlab=xlabel, ylab=ylabel, shape = PDB, alpha = I(txtalpha)) +
          geom_point(aes(color = PDB), alpha = 0.6) +
          scale_colour_manual(name="", values = c("1G9O"="orange", "3QDO"="blue", "3"="red", "value3"="grey", "value2"="black")) +
          labs(title = "%(plot_title)s") +
          theme(plot.title = element_text(color = "#555555", size=rel(0.75))) +
  
          # Correlation fit lines (global + one per facet
          geom_abline(size = 0.125, color="black", intercept = lmv_intercept, slope = lmv_yvalues, alpha=0.2) +
          geom_abline(size = 0.125, color="orange", intercept = lmv_intercept_NHERF1, slope = lmv_yvalues_NHERF1, alpha=0.4) +
          geom_abline(size = 0.125, color="blue", intercept = lmv_intercept_SNX27, slope = lmv_yvalues_SNX27, alpha=0.4) +
  
          geom_abline(slope=1, intercept=0, linetype=3, size=0.25, alpha=0.4) + # add a diagonal (dotted)
          coord_cartesian(xlim = c(minx, maxx), ylim = c(miny, maxy)) + # set the graph limits
  
          geom_text(hjust = 0, size=1.5, color="#000000", alpha=0.6, data=subset(xy_data, abs(yvalues - xvalues) > 2 & xvalues <= 0), aes(xvalues, yvalues+0.35, label=Origin_of_peptide), check_overlap = TRUE) + # label outliers
          geom_text(hjust = 1, size=1.5, color="#000000", alpha=0.6, data=subset(xy_data, abs(yvalues - xvalues) > 2 & xvalues > 0), aes(xvalues, yvalues+0.35, label=Origin_of_peptide), check_overlap = TRUE) + # label outliers
          geom_text(hjust=0, size=2, colour="black", aes(x = xpos, y = ypos, label = sprintf("R == %%0.2f", round(rvalue, digits = 4))), parse = TRUE) +
          geom_text(hjust=0, size=2, colour="darkorange", aes(x = xpos, y = ypos_NHERF1, label = sprintf("R[NHERF] == %%0.2f", round(rvalue_NHERF1, digits = 4))), parse = TRUE) +
          geom_text(hjust=0, size=2, colour="blue", aes(x = xpos, y = ypos_SNX27, label = sprintf("R[SNX27] == %%0.2f", round(rvalue_SNX27, digits = 4))), parse = TRUE) +
          theme(legend.position = "none")
  #       geom_text(hjust=0, size=2, colour="black", aes(xpos, ypos, fontface="plain", family = "sans", label=paste(sprintf("R = %%0.2f%%s", round(rvalue, digits = 4), lrt), expression('R'[3])  ))) # add correlation text; hjust=0 sets left-alignment
  
  #geom_text(hjust=0, size=3, colour="black", aes(xpos, ypos, fontface="plain", family = "sans", label=sprintf("R = %%0.2f", round(rvalue, digits = 4)))) # add correlation text; hjust=0 sets left-alignment
  #       geom_text(hjust=0, size=3, colour="black", aes(xpos, ypos, fontface="plain", family = "sans", label=sprintf("R = %%0.2f", round(rvalue, digits = 4)))) # add correlation text; hjust=0 sets left-alignment
  
  # Plot graph
  p
  dev.off()
          '''
  
  
      #geom_point(aes(color = C)) +
      #color = "%(series_color)s"
  
      # Create the R script
      plot_type = 'png'
      png_plot_commands = single_plot_commands % locals()
      boxplot_r_script = boxplot_r_script % locals()
      r_script_filename = '{0}.R'.format(file_prefix)
      r_script_filepath = os.path.join(output_directory, r_script_filename)
      write_file(r_script_filepath, boxplot_r_script)
  
      # Run the R script
      run_r_script(r_script_filename, cwd = output_directory)