////////////////////////////////////////////////////////////
// ScoreBPSEQ.cc
////////////////////////////////////////////////////////////

#include <iostream>
#include <map>
#include <cassert>
#include <cctype>
#include <vector>
#include <string>
#include <sstream>
#include <fstream>
#include <cstdio>
#include <cstdlib>

using namespace std;

bool toggle_complementary_only = false;
vector<string> filenames;

////////////////////////////////////////////////////////////
// Complementary
////////////////////////////////////////////////////////////

bool Complementary (char c, char d){
  if (c == 'A' && d == 'U') return true;
  if (c == 'G' && d == 'U') return true;
  if (c == 'G' && d == 'C') return true;
  if (d == 'A' && c == 'U') return true;
  if (d == 'G' && c == 'U') return true;
  if (d == 'G' && c == 'C') return true;
  return false;
}

////////////////////////////////////////////////////////////
// Trim()
////////////////////////////////////////////////////////////

string Trim (const string &s){
  int left = 0, right = (int) s.length();
  while (left < right && isspace(s[left])) left++;
  while (left < right && isspace(s[right-1])) right--;
  return s.substr(left,right-left);
}

////////////////////////////////////////////////////////////
// ReadBPSEQ()
////////////////////////////////////////////////////////////

vector<int> ReadBPSEQ (const string &filename){
  ifstream infile (filename.c_str());
  if (infile.fail()) return vector<int>();
  
  vector<int> mapping (1, 0);
  vector<char> letters (1, 0);
  int a, b;
  char ch;

  while (infile >> a >> ch >> b){
    assert (a == (int) mapping.size());
    mapping.push_back (b);
    letters.push_back (ch);
  }
  infile.close();
  
  for (int i = 1; i < (int) mapping.size(); i++){

    // eliminate non-complementary base pairs, if desired
    
    if (toggle_complementary_only && mapping[i] != 0 && 
	!Complementary(letters[i], letters[mapping[i]])){
      mapping[mapping[i]] = 0;
      mapping[i] = 0;
    }
  }
  
  return mapping;
}

////////////////////////////////////////////////////////////
// ScoreFile()
////////////////////////////////////////////////////////////

vector<double> ScoreFile (const string &ref_filename, 
                          const string &test_filename){
  
  // read BPSEQ files
  
  vector<int> ref = ReadBPSEQ (ref_filename);
  vector<int> test = ReadBPSEQ (test_filename);

  if (ref.size() == 0) return vector<double>();
  if (ref.size() != test.size()){
    cerr << "ERROR: " << ref_filename << " (" << ref.size() << ") and " 
	 << test_filename << " (" << test.size() << ") have different lengths!" << endl;
    exit(1);
  }
  
  int num_ref_pairs = 0;
  int num_test_pairs = 0;
  int correctly_paired = 0;

  for (int i = 1; i < (int) test.size(); i++){
    if (test[i] != 0) num_test_pairs++;
    if (ref[i] != 0) num_ref_pairs++;
    if (test[i] != 0 && test[i] == ref[i]) correctly_paired++;
  }

  vector<double> ret;
  ret.push_back ((num_ref_pairs == 0) ? 1 : (double) correctly_paired / num_ref_pairs);
  ret.push_back ((num_test_pairs == 0) ? 1 : (double) correctly_paired / num_test_pairs);
  return ret;
}

////////////////////////////////////////////////////////////
// Version()
//
// Display program version.
////////////////////////////////////////////////////////////
 
void Version(){
  cerr << "ScoreBPSEQ version 1.00 - Compare two BPSEQ files" << endl << endl
       << "Written by Chuong B. Do" << endl;
  exit(0);
}

////////////////////////////////////////////////////////////
// Usage()
////////////////////////////////////////////////////////////

void Usage(){
  cerr << endl
       << "Usage: score_bpseq [OPTION]... REF_BPSEQ TEST_BPSEQ" << endl
       << endl
       << "       where [OPTION]...  is a list of zero or more optional arguments" << endl
       << "             REF_BPSEQ    is the name of the reference BPSEQ file" << endl
       << "             TEST_BPSEQ   is the name of the test BPSEQ file" << endl
       << endl
       << "Miscellaneous arguments:" << endl
       << "  --version               display program version information" << endl
       << "  --complementary         include only complementary (A-U, C-G, G-U) base pairs" << endl
       << endl;
  exit (1);
}

////////////////////////////////////////////////////////////
// ParseParameters()
////////////////////////////////////////////////////////////

void ParseParameters (int argc, char **argv){
  
  if (argc < 3) Usage();
  filenames.clear();

  for (int argno = 1; argno < argc; argno++){
    if (!strcmp (argv[argno], "--version")){
      Version();
    } else if (!strcmp (argv[argno], "--complementary")){
      toggle_complementary_only = true;
    } else {
      filenames.push_back (argv[argno]);
    }
  }
  
  if (filenames.size() != 2){
    cerr << "ERROR: Incorrect number of filenames." << endl;
    exit (1);
  }   
}

////////////////////////////////////////////////////////////
// main()
////////////////////////////////////////////////////////////

int main (int argc, char **argv){
  ParseParameters (argc, argv);
  
  vector<double> results = ScoreFile (filenames[0], filenames[1]);
  cout << "sensitivity=" << results[0] << "; ppv=" << results[1] << ";" << endl;
}