#include <stdio.h>
#include <string>
#include <vector>
#include <set>
#include <iterator>
#include <fstream>
#include <iostream>
#include <sys/timeb.h>

#include "RecoReadout/RecoElem.hh"
#include "RecoReadout/Reaction.hh"
#include "MExplorer/xduGen.hh"

#include "sbml/SBMLReader.h"
#include "sbml/SBMLTypes.h"


using namespace std;

//************************************************************************


RecoElem* SBMLtoRecoElem(Model_t* m); 
void inreaction(Reaction_t* r, vector<Reaction>& rxns);
void printList (RecoElem& elem);
//************************************************************************




//************************************************************************
int
main(int argc, char* argv[] ) {
//************************************************************************

 
  SBMLDocument_t *d;
  Model_t *m;

  if (argc != 2)
  {
    printf("usage: readSBML <filename>\n");
    return 1;
  }


  d     = readSBML(argv[1]);

  m = d->model;
  if (0==m) {
    cout << "Did not read model " << argv[1] << endl;
    exit( 1 );
  }

  RecoElem* elem = SBMLtoRecoElem (m); 

  cout << "************************" << endl;
  printList(*elem);
  cout << "************************" << endl;

   std::string path="";
  xduGen( *elem, path );

  return 0;
}




//************************************************************************
void printList (RecoElem& elem)
//************************************************************************
{

  // I don't think there's any subelements. Model w/ one subelement, and 1 kid

  vector<RecoElem>& subElem   = elem.getSubElements();
  vector<Reaction>& rxns      = elem.getReactions();

  int nrxns= rxns.size ();
  int nkids= subElem.size ();

  cout << "MY NAME: " << elem.getName() << endl;
  cout << "  No. of reactions " << nrxns << endl;
  cout << "  No. of kids: " << nkids << endl;


  for (int l=0; l< nrxns ; l++) {
    cout << rxns[l].asString( ) << endl; 
  }


  for (int i=0; i<nkids; i++){
    cout << "in deeper  ****  " << endl;
    printList(subElem[i]);
    cout << "back up *** " << endl;
  }

  return;
}





//**********************************************************************
RecoElem*
SBMLtoRecoElem(Model_t* m) {
//************************************************************************

  string name;
  int id, numreactant, numproduct;

  Reaction_t* r;  
  Species_t* s;
  List_t* rlist = m->reaction;
 

  //not sure what to put as an id, 0 for now
  name=m->name;
  id=0;

  RecoElem* elem= new RecoElem( name, id);

  vector<Reaction>& rxns      = elem->getReactions();
 

  //put in reactions
  int reactionsize= List_size(rlist);
  //cout << "number of reactions " << reactionsize << endl;
  
  

  for (int j=0; j< reactionsize ; j++)
    {
      
      //  cout << "Reaction : " << j <<endl;
      r=Model_getReaction(m,j);
      //cout << r->name << endl;
      inreaction(r,rxns);

      // cout << endl;
    }


  return elem;
  

}
//************************************************************************
void
inreaction(Reaction_t* r, vector<Reaction>& rxns) {
//************************************************************************

  int id;
  int dir;

  int numreactant=Reaction_getNumReactants(r);
  int numproduct=Reaction_getNumProducts(r);
  SpeciesReference_t* sr, *pr;
   
  //id not sure what to put here
  //direction is reversible  or not. false, true

  Reaction::Direction reDir;
  id=0;
  dir=r->reversible;
  Reaction rxn( id, Reaction::forward );
  if ( dir ) {
    rxn.setDirection( Reaction::both );
  }

  vector<string>& reactant = rxn.getInputMetaboliteNames();
  vector<string>& product  = rxn.getOutputMetaboliteNames();

  vector<int>& inMetabN = rxn.getInputStoich();
  vector<int>& outMetabN = rxn.getOutputStoich();
  
  //add reactants
  //add products
  
 

  //figuring out the common denominator
  int rdem = 1;
  int pdem = 1;
  int demstoich;

  for (int j=0; j< numreactant ; j++ )
    {
      sr=Reaction_getReactant(r,j);
      rdem= rdem * sr->denominator;
      
    }

  
  for (int m=0; m<numproduct; m++)
    {
      pr=Reaction_getProduct(r,m);
      pdem= pdem * pr->denominator;

    }


  demstoich = rdem * pdem;
  //  cout << "common denominator: " << demstoich;

   for (int k=0; k < numreactant ; k++ )
	 {
	   sr=Reaction_getReactant(r,k);
	   // cout << "test reactant: " << sr->species << endl;

		  
	   reactant.push_back (sr->species);
           inMetabN.push_back( demstoich * sr->stoichiometry/sr->denominator);
	 }

   //cout << rxn.asString()[0];

   for (int l=0; l<numproduct; l++)
	 {

	   pr=Reaction_getProduct(r,l);
	   //  cout << "product: " << pr->species  << endl;
	 
	   product.push_back(pr->species);
	   outMetabN.push_back(demstoich * pr->stoichiometry/ pr->denominator );
	 }
   
   rxns.push_back(rxn);
  
  return;
}