/* massXpert - the true massist's program.
   --------------------------------------
   Copyright(C) 2006,2007 Filippo Rusconi

   http://www.massxpert.org/massXpert

   This file is part of the massXpert project.

   The massxpert project is the successor to the "GNU polyxmass"
   project that is an official GNU project package(see
   www.gnu.org). The massXpert project is not endorsed by the GNU
   project, although it is released ---in its entirety--- under the
   GNU General Public License. A huge part of the code in massXpert
   is actually a C++ rewrite of code in GNU polyxmass. As such
   massXpert was started at the Centre National de la Recherche
   Scientifique(FRANCE), that granted me the formal authorization to
   publish it under this Free Software License.

   This software is free software; you can redistribute it and/or
   modify it under the terms of the GNU  General Public
   License version 3, as published by the Free Software Foundation.
   

   This software is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this software; if not, write to the

   Free Software Foundation, Inc.,

   51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
*/


/////////////////////// Local includes
#include "cleaver.hpp"


namespace massXpert
{

  Cleaver::Cleaver(Polymer *polymer,
		    const PolChemDef *polChemDef, 
		    const CleaveOptions &cleaveOptions,
		    const CalcOptions &calcOptions, 
		    const IonizeRule &ionizeRule)
    : mp_polymer(polymer), 
      mp_polChemDef(polChemDef),
      m_cleaveOptions(cleaveOptions),
      m_calcOptions(calcOptions),
      m_ionizeRule(ionizeRule)
  {
    Q_ASSERT(mp_polymer && mp_polChemDef);
    mp_oligomerList = 0;
  }
  
  
  Cleaver::Cleaver(const Cleaver &other)
    : mp_polymer(other.mp_polymer), 
      mp_polChemDef(other.mp_polChemDef), 
      m_cleaveOptions(other.m_cleaveOptions),
      m_calcOptions(other.m_calcOptions),
      m_ionizeRule(other.m_ionizeRule)
  {
    Q_ASSERT(mp_polymer && mp_polChemDef);
    mp_oligomerList = 0;
  }


  Cleaver::~Cleaver()
  {
    // We are not owner of the oligomer list, do not free it!
  }


  void 
  Cleaver::setOligomerList(OligomerList *oligomerList)
  {
    Q_ASSERT(oligomerList);
    
    mp_oligomerList = oligomerList;
  }


  OligomerList *
  Cleaver::oligomerList()
  {
    return mp_oligomerList;
  }


  bool
  Cleaver::cleave(bool reset)
  {
    if (!mp_oligomerList)
      qFatal("%s@%d -- The oligomer list is 0.", __FILE__, __LINE__);
    
    // If the polymer sequence is empty, just return.
    if (!mp_polymer->size())
      return true;
  
    // Ensure that the cleavage pattern was already parsed.
    
    if (!m_cleaveOptions.motifList()->size())
      {
	if (!m_cleaveOptions.parse())
	  {
	    qDebug() << __FILE__ << __LINE__ 
		     << "Failed to parse the cleavage options";
	    
	    return false;
	  }
      }
    
    //   qDebug() << __FILE__ << __LINE__
    // 	    << "number of motifs:" 
    // 	    << mp_cleaveOptions->motifList()->size();
    
    if (!fillIndexLists())
      {
	qDebug() << __FILE__ << __LINE__ 
		 << "Index lists(cleave/nocleave) are empty." 
	  "No oligomer generated.";
	
	// We can return true, as no error condition was found but not
	// oligomers were generated.

	return true;
      }
    
    if (resolveCleavageNoCleavage() == -1)
      {
	qDebug() << __FILE__ << __LINE__ 
		 << "Failed to resolve cleavage/nocleavage";
	
	return false;
      }

    removeDuplicatesCleavage();
        
    qSort(m_cleaveIndexList.begin(), m_cleaveIndexList.end());

//     for (int debugIter = 0; debugIter < m_cleaveIndexList.size(); 
// 	 ++debugIter)
//       {
// 	qDebug() << __FILE__ << __LINE__
// 		 << "Index:" << m_cleaveIndexList.at(debugIter);
//       }
    
    if (reset)
      emptyOligomerList();
    
    for (int iter = 0 ; iter <= m_cleaveOptions.partials(); ++iter)
      {
	if (cleavePartial(iter) == -1)
	  {
	    qDebug() << __FILE__ << __LINE__ 
		 << "Failed to perform partial cleavage at index:"
		     << iter;
	    
	    return false;
	  }
      }
    
    // At this point we have the list of lists of oligomers, one list of
    // oligomers for each partial cleavage.

    while (m_cleaveIndexList.size())
      m_cleaveIndexList.removeFirst();
  
    while (m_noCleaveIndexList.size())
      m_noCleaveIndexList.removeFirst();
  
    return true;
  }


  int
  Cleaver::fillIndexLists()
  {
    QList<CleaveMotif *> *motifList = m_cleaveOptions.motifList();
  

    while (m_cleaveIndexList.size())
      m_cleaveIndexList.removeFirst();
  
    while (m_noCleaveIndexList.size())
      m_noCleaveIndexList.removeFirst();

    // Since version 2.3.0, the cleavage might be performed on a
    // selected portion of a sequence only, not necessarily on the
    // whole polymer sequence.
  
    CoordinateList coordinateList = m_calcOptions.coordinateList();
    Coordinates *coordinates = coordinateList.first();
    
    int startIndex = coordinates->start();
    int endIndex = coordinates->end();

    for (int iter = 0; iter < motifList->size(); ++iter)
      {
	CleaveMotif *cleaveMotif = motifList->at(iter);
	Q_ASSERT(cleaveMotif);
      
	int index = startIndex -1;
      
	while (1)
	  {
	    index = findCleaveMotif(*cleaveMotif, index + 1, endIndex);
	  
	    if (index == -1)
	      break;
	  
	    // Do not forget: The position at which the motif is found
	    // in the polymer sequence is not necessarily the position
	    // at which the cleavage will effectively occur. Indeed,
	    // let's say that we found such motif in the polymer
	    // sequence: "KKRKGP". This motif was extracted from a
	    // cleave spec that had a pattern like this: "KKRK/GP". What
	    // we see here is that the cleavage occurs after the fourth
	    // monomer! And we must realize that the 'index' returned
	    // above corresponds to the index of the first 'K' in
	    // "KKRKGP" motif that was found in the polymer
	    // sequence. Thus we have to take into account the offset
	    //(+4, in our example, WHICH IS A POSITION and not an
	    // index, which is why we need to remove 1 below) of the
	    // cleavage:
	 
	    int cleavageIndex = index + cleaveMotif->offset() - 1;
	  
	    if (cleavageIndex < 0)
	      continue;

	    if (cleavageIndex >= endIndex)
	      break;
	  
// 	    qDebug() << __FILE__ << __LINE__ 
// 		      << "Found new cleavage index:"
// 		      << cleavageIndex;
	    
	    if (cleaveMotif->isForCleave())
	      {
		m_cleaveIndexList.append(cleavageIndex);
		
// 		qDebug() << __FILE__ << __LINE__ 
// 			  << "For cleavage, index:" << cleavageIndex;
	      }
	    else
	      {
		m_noCleaveIndexList.append(cleavageIndex);
		
// 		qDebug() << __FILE__ << __LINE__ 
// 			  << "Not for cleavage";
	      }
	    
	  }
	// End of 
	// while (1)
      }
    // End of
    // for (int iter = 0; iter < motifList->size(); ++iter)

    // Note that returning 0 is not an error condition, because a
    // sequence where no site is found whatsoever will result in 0.
    return m_cleaveIndexList.size() +  m_noCleaveIndexList.size();
  }


  int
  Cleaver::resolveCleavageNoCleavage()
  {
    for (int iter = 0; iter < m_noCleaveIndexList.size(); ++iter)
      {
	int noCleaveIndex = m_noCleaveIndexList.at(iter);
      
	for(int jter = 0; jter < m_cleaveIndexList.size(); ++jter)
	  {
	    int cleaveIndex = m_cleaveIndexList.at(jter);
	  
	    if (noCleaveIndex == cleaveIndex)
	      m_cleaveIndexList.removeAt(jter);
	  }
      }
  
    return m_cleaveIndexList.size();  
  }


  int
  Cleaver::removeDuplicatesCleavage()
  {
    for (int iter = 0; iter < m_cleaveIndexList.size(); ++iter)
      {
	int index = m_cleaveIndexList.at(iter);

	int foundItemIndex = m_cleaveIndexList.indexOf(index, iter + 1);
	
	if (foundItemIndex != -1)
	  {
	    m_cleaveIndexList.removeAt(foundItemIndex);
	    --iter;
	  }
      }

    return m_cleaveIndexList.size();
  }
  
  
  int 
  Cleaver::findCleaveMotif(CleaveMotif &cleaveMotif, 
			   int startIndex, int endIndex)
  {
    bool noGood = false;

    int firstIndex = 0;
  
    QList<const Monomer *> monomerList = mp_polymer->monomerList();
    const QStringList &codeList = cleaveMotif.codeList();
  

    //We have to iterate in the polymer sequence starting at 'index', in
    //search for a sequence element identical to the sequence in the
    //'cleaveMotif'.
  
    // This means that if 

    // cleavemotif->m_motifList [0] = "Lys"

    // cleavemotif->m_motifList [1] = "Pro"
  
    // the, we want to search in the polymer list of monomers the same
    // sequence by iterating in this list from index 'index' onwards,
    // and we stop searching when the list's end is found or if 

    // list [n] = "Lys" and

    // list [n+1] = "Pro".
  

    if (mp_polymer->size() == 0)
      return 0;

    if (codeList.size() == 0)
      return -1;
  
    if (startIndex < 0)
      return -1;
    
    if (endIndex >= mp_polymer->size())
      return -1;
        
    // Seed the routine by setting 'first' to the first motif in the
    // codeList (in our example this is "Lys").

    QString firstCode = codeList.first();
  
    // And now iterate (starting from 'index') in the polymer
    // sequence's list of monomers in search for a monomer having the
    // proper code ("Lys").
  
    int iterIndex = startIndex;
 
    while (iterIndex < endIndex)
      {
	const Monomer *monomer = monomerList.at(iterIndex);

	if (monomer->code() != firstCode)
	  {
	    // The currently iterated code is not the one we
	    // search. So go one code further in the sequence.
	    
	    ++iterIndex;
	    continue;
	  }
	
	// If we are here, then that means that we actually found on
	// monomer code in the sequence that matches the one we are
	// looking for.
	
	firstIndex = iterIndex;
	noGood = false;
	
	// Now that we have anchored our search at firstIndex in the
	// polymer sequence, continue with next monomer and check if
	// it matches the next monomer in the motif we are looking
	// for.
  
	for (int iter = 1; iter < codeList.size(); ++iter)
	  {
	    if (iterIndex + iter >= endIndex)
	      {
		noGood = true;
		break;
	      }
	      
	    QString nextCode = codeList.at(iter);
	      
	    monomer = monomerList.at(iterIndex + iter);

	    if (monomer->code() == nextCode)
	      continue;
	    else
	      {
		noGood = true;
		break;
	      }
	  }
	// End of
	// for (int iter = 1; iter < codeList.size(); ++iter)

	if (noGood)
	  {
	    ++iterIndex;
	    continue;
	  }
	else
	  {
	    return firstIndex;
	  }
      }
    // End of 
    // while (iterIndex < monomerList.size())
		  
    return -1;
  }


  bool 
  Cleaver::accountCleaveRule(CleaveRule *cleaveRule, 
			     CleaveOligomer *oligomer)
  {
    Q_ASSERT(cleaveRule);
    Q_ASSERT(oligomer);
  
    const QList<Atom *> &refList = 
      oligomer->polymer()->polChemDef()->atomList();

    // For each Coordinates element in the oligomer, we have to ensure
    // we apply the formula(s) that is/are required.

    int coordinatesCount = static_cast<CoordinateList *>(oligomer)->size();

//     qDebug() << __FILE__ << __LINE__
// 	      << "Coordinates count:" << coordinatesCount;
    
    for (int iter = 0; iter < coordinatesCount; ++iter)
      {
	Coordinates *coordinates = 
	  static_cast<CoordinateList *>(oligomer)->at(iter);
	
	if (!cleaveRule->leftCode().isEmpty())
	  {
	    // The formula has to be there.
	    
	    Formula ruleFormula = cleaveRule->leftFormula();
	    Q_ASSERT(!ruleFormula.formula().isEmpty());
	    
	    // What is the monomer at the left end of current oligomer ?
	    const Monomer &monomer = oligomer->atLeftEnd();
	    
	    if (monomer.code() == cleaveRule->leftCode())
	      {
// 		qDebug() << __FILE__ << __LINE__
// 			  << "Matched left code:" << cleaveRule->leftCode();
		
		// But, this is not going to be real true, if the
		// monomer is acutally the left-end monomer of the
		// polymer sequence, because then that would mean that
		// there was no cleavage on this monomer, thus no rule
		// to apply.
		
		if (!coordinates->start())
		  {
		    // The monomer is not the left-end monomer, so the
		    // match is real. Account for the formula !
		    
		    if (!ruleFormula.accountMasses(refList, oligomer))
		      return false;
		  }
	      }
	  }
  
	if (!cleaveRule->rightCode().isEmpty())
	  {
	    // The formula has to be there.
      
	    Formula ruleFormula = cleaveRule->rightFormula();
	    Q_ASSERT(!ruleFormula.formula().isEmpty());
      
	    // What is the monomer at the right end of current oligomer ?
	    const Monomer &monomer = oligomer->atRightEnd();
      
	    if (monomer.code() == cleaveRule->rightCode())
	      {
// 		qDebug() << __FILE__ << __LINE__
// 			  << "Matched right code:" << cleaveRule->rightCode();
	  
		// But, this is not going to be real true, if the
		// monomer is acutally the right-end monomer of the
		// polymer sequence, because then that would mean that
		// there was no cleavage on this monomer, thus no rule
		// to apply.

		if (coordinates->end() != mp_polymer->size() - 1)
		  {
		    // The monomer is not the right-end monomer, so
		    // the match is real. Account for the formula !
		    
		    if (!ruleFormula.accountMasses(refList, oligomer))
		      return false;
		  }
	      }
	  }
      }
      
    return true;
  }


  int
  Cleaver::cleavePartial(int partialCleavageValue)
  {
    bool oligomerIsPolymer = false;
  
    int iter = 0;
  
    static int leftIndex = 0;
    static int rightIndex = 0;
  
    Q_ASSERT(partialCleavageValue >= 0);

    OligomerList partialOligomerList;

    // Since version 2.3.0, the cleavage might be performed on a
    // selected portion of a sequence only, not necessarily on the
    // whole polymer sequence. We have to know these values because
    // when we create oligomer we'll have to set the coordinates of
    // the monomers correctly.
  
    CoordinateList coordinateList = m_calcOptions.coordinateList();
    Coordinates *coordinates = coordinateList.first();
    
    int startIndex = coordinates->start();
    int endIndex = coordinates->end();

    leftIndex = startIndex;
    rightIndex = 0;

    // Iterate in the array of indices where the cleavages should occur.
    
    for (iter = 0; iter < m_cleaveIndexList.size(); ++iter)
      {
	// Make sure, if the partial cleavage is very large, for
	// example, that it will not lead us to access the polymer
	// sequence at a position larger than its upper boundary. 

	// Imagine cutting a polymer with only one Met residue with
	// cyanogen bromide: m_cleaveIndexList will contain a single
	// element: the index at which the methionine occurs in the
	// polymer sequence(and there is a single one). Now, Imagine
	// that we are asked to perform a cleavage with
	// 'partialCleavageValue' of 2. The way we do it is that we fetch the
	// index in the list of cleavage indices(m_cleaveIndexList) two
	// positions farther than the position we are iterating:

	// int partCleave = iter + partialCleavageValue;
      
	// Now, if m_cleaveIndexList contains a single element, asking
	// for this m_cleaveIndexList.at(iter + partialCleavageValue) will
	// go out of the boundaries of the list, since it has a single
	// item and partialCleavageValue is 2. This is what we are willing to
	// avoid.
      
	int partCleave = iter + partialCleavageValue;
      
	if (partCleave >= m_cleaveIndexList.size())
	  {
	    if (iter == 0)
	      oligomerIsPolymer = true;
	  
	    break;
	  }
      
	rightIndex = m_cleaveIndexList.at(partCleave);

	QString name = QString("%1#%2")
	  .arg(partialCleavageValue)
	  .arg(iter + 1);
 
      	CleaveOligomer *oligomer = 
	  new CleaveOligomer(mp_polymer, 
			      name, m_cleaveOptions.name(),
			      Ponderable(),
			      leftIndex, rightIndex, partialCleavageValue);
	
	partialOligomerList.append(oligomer);

	leftIndex = m_cleaveIndexList.at(iter) + 1;
      }
    // End of 
    // for (int iter = 0; iter < m_cleaveIndexList.size(); iter=+)
  
    // Do not forget the right-end oligomer ! And be sure to determine
    // what's its real left end index !
  
    if (oligomerIsPolymer)
      leftIndex = startIndex;
    else
      leftIndex = m_cleaveIndexList.at(--iter) + 1;
  
    // 'iter' is used to construct the name of the oligomer, so we have
    // to increment it once because we did not have the opportunity to
    // increment it between the last but one oligomer and this one.
  
    ++iter;
  
    // Remove 1 because the oligomer is described using indices and not
    // positions.
  
    rightIndex = endIndex;
  
    QString name = QString("%1#%2")
      .arg(partialCleavageValue)
      .arg(iter + 1);
  
    CleaveOligomer *oligomer = 
      new CleaveOligomer(mp_polymer, 
			  name, m_cleaveOptions.name(),
			  Ponderable(),
			  leftIndex, rightIndex, partialCleavageValue);
  
    partialOligomerList.append(oligomer);
  
    // At this point all the skeleton oligomers have been computed for
    // the given partialCleavageValue. We still have to perform the
    // cross-link analysis prior to both calculate the masses and
    // perform the ionization of all the generated oligomers. Note that
    // making cross-link analysis is only useful in case the cleavage is
    // full(partialCleavageValue == 0).

    if (!partialCleavageValue)
      {
	if (m_calcOptions.monomerEntities() & MXT_MONOMER_CHEMENT_CROSS_LINK)
	  {
	    if (analyzeCrossLinks(&partialOligomerList) == -1)
	      {
		return false;  
	      }
	  }
      }
    
    // Finally, we can now perform the mass calculations and the
    // ionization. We will use each oligomer in the oligomerList as a
    // template for creating new oligomers(with different z values) and
    // all the new oligomers will be appended to waitOligomerList. Each
    // time a template oligomer will have been used, it will be removed
    // from oligomerList. Once all the oligomers in oligomerList will
    // have been used, and thus removed, all the newly allocated
    // oligomers in waitOligomerList will be moved to oligomerList.

    OligomerList waitOligomerList;

    while (partialOligomerList.size())
      {
	CleaveOligomer *iterOligomer = 
	  static_cast<CleaveOligomer *>(partialOligomerList.takeFirst());
      
	// We do not ask that the oligomer be ionized yet, because we
	// have to first account for potential cleavage rules! Thus we
	// pass an uninitialized ionization rule with IonizeRule().
	// iterOligomer->calculateMasses(*mp_calcOptions,
	// *mp_ionizeRule); This was a bug in the release versions up
	// to 1.6.1.
	iterOligomer->calculateMasses(&m_calcOptions);
      
	// At this point we should test if the oligomer has to be
	// processed using cleavage rules.

	for(int jter = 0; jter < m_cleaveOptions.ruleList()->size(); ++jter)
	  {
	    //  Note that the accounting of the cleavage rule is
	    //  performed as if the oligomer was charged 1. This is why
	    //  we have to ionize the oligomer only after we have
	    //  completed the determination of its atomic composition.

	    CleaveRule *cleaveRule = m_cleaveOptions.ruleList()->at(jter);
	  
// 	    qDebug() << __FILE__ << __LINE__
// 		      << "Accounting for cleaverule:"
// 		      << cleaveRule->name();
	    
// 	    qDebug() << __FILE__ << __LINE__
// 		      << "Oligomer mono mass before:" << iterOligomer->mono();
	    
	    if (!accountCleaveRule(cleaveRule, iterOligomer))
	      return -1;

// 	    qDebug() << __FILE__ << __LINE__
// 		      << "Oligomer mono mass after:" << iterOligomer->mono();
	  }
      
	// At this point we can finally ionize the oligomer ! Remember
	// that we have to ionize the oligomer as expected in the
	// cleavage options. Because the ionization changes the values
	// in the oligomer, and we need a new oligomer each time, we
	// duplicate the oligomer each time we need it.

	int startIonizeLevel = m_cleaveOptions.startIonizeLevel();
	int endIonizeLevel = m_cleaveOptions.endIonizeLevel() + 1;
	IonizeRule ionizeRule(m_ionizeRule);

	for(int kter = startIonizeLevel; kter < endIonizeLevel; ++kter)
	  {
	    ionizeRule.setLevel(kter);

	    CleaveOligomer *newOligomer = new CleaveOligomer(*iterOligomer);
	  
	    if (newOligomer->ionize(ionizeRule) == -1)
	      {
		delete newOligomer;
	      
		return -1;
	      }
	
	    // The name was set already during the creation of the
	    // template oligomer. All we have to add to the name is the
	    // ionization level.

	    QString name = iterOligomer->name() + 
	      QString("#z=%3").arg(newOligomer->charge());
	  
	    newOligomer->setName(name);
	  
	    waitOligomerList.append(newOligomer);
	  }
      
	// We can delete the template oligomer that was already removed
	// from the oligomerList(use of QList::takeFirst()).
	delete iterOligomer;
      }
  
    // At this point we should transfer all the oligomers from the
    // waitOligomerList to the initial oligomerList.

    while (waitOligomerList.size())
      {
	Oligomer *iterOligomer = waitOligomerList.takeFirst();
      
	partialOligomerList.append(iterOligomer);
      }

    // Finally transfer all the oligomers generated for this partial
    // cleavage to the list of ALL the oligomers.

    int oligomerCount = partialOligomerList.size();
    
    while (partialOligomerList.size())
      {
	Oligomer *iterOligomer = partialOligomerList.takeFirst();
      
	mp_oligomerList->append(iterOligomer);
      }
    
    return oligomerCount;
  }
  

  int
  Cleaver::analyzeCrossLinks(OligomerList *oligomerList)
  {
    Q_ASSERT(oligomerList);

    QList<Oligomer *> crossLinkedOligomerList;

    // General overview:

    // Iterate in the polymer's list of cross-links and for each
    // cross-link find the oligomer that contains the first monomer
    // involved in the cross-link. This first found oligomer should
    // serve as a seed to pull-down all the oligomers cross-linked to
    // it.

    const CrossLinkList &crossLinkList = mp_polymer->crossLinkList();
      
    for (int iter = 0; iter < crossLinkList.size(); ++iter)
      {
	CrossLink *crossLink = crossLinkList.at(iter);
      
	// With that crossLink, find an oligomer that encompasses the
	// first monomer of the cross-link.
      
	const Monomer *firstMonomer = crossLink->firstMonomer();
      
	Q_ASSERT(firstMonomer);
      
	// What oligomer does encompass that monomer ?
      
	int foundIndex = 0;
      
	Oligomer *firstOligomer = 
	  oligomerList->findOligomerEncompassing(firstMonomer, &foundIndex);
      
	if (firstOligomer)
	  {
	    // At this point we should turn this oligomer into a
	    // cross-linked oligomer, so that we can continue performing its
	    // cross-link analysis. To do that we allocate a list of
	    // oligomers for this cross-linked oligomer, were we'll store
	    // this first oligomer and then all the "pulled-down" oligomers.
          
	    // Remove the cross-link from the main list of oligomers so
	    // that we do not stumble upon it in the next analysis
	    // steps.
	    oligomerList->removeAt(foundIndex);

	    // Set the cross-linked oligomer apart.
	    crossLinkedOligomerList.append(firstOligomer);

	    // Finally deeply scrutinize the oligomer.
	    analyzeCrossLinkedOligomer(firstOligomer, oligomerList); 
	  }
	else
	  {
	    // 	  qDebug() << __FILE__ << __LINE__ 
	    // 		    << "Cross-link at index" << iter
	    // 		    << "did not find any oligomer for its first monomer "
	    // 	    "partner";
	  }
      }
  
    // At this point we have terminated analyzing all the oligomers for
    // the partial cleavage. All we have to do is move all the
    // crossLinked oligomers from the crossLinkedOligomerList to
    // oligomerList.

    while (crossLinkedOligomerList.size())
      {
	Oligomer *oligomer = crossLinkedOligomerList.takeAt(0);
      
	oligomerList->append(oligomer);
      }
  
    crossLinkedOligomerList.clear();

    // Return the number of cross-linked/non-cross-linked oligomers
    // alltogether.

    return oligomerList->size();
  }


  int
  Cleaver::analyzeCrossLinkedOligomer(Oligomer *oligomer,
				       OligomerList *oligomerList)
  {
    Q_ASSERT(oligomer);
    Q_ASSERT(oligomerList);
  
    const CrossLinkList &crossLinkList =  mp_polymer->crossLinkList();
  
    OligomerList clearanceOligomerList;

    // 'oligomer' is the first oligomer in the cross-link series of
    // oligomers. It is the "seeding" oligomer with which to pull-down
    // all the others. Prepend to its name the "cl-" string to let it
    // know it is cross-linked.

    QString name = oligomer->name();
    name.prepend("cl-");
    oligomer->setName(name);
  
    // Iterate in the 'oligomer' and for each monomer get any
    // cross-linked oligomer out of the list of cross-links.

    for (int iter = oligomer->startIndex();
	 iter < oligomer->endIndex() + 1; ++iter)
      {
	const Monomer *monomer = mp_polymer->at(iter);

	// What crossLinks do involve this monomer ?

	QList<int> crossLinkIndices;
      
	int ret = 
	  crossLinkList.crossLinksInvolvingMonomer(monomer,
						    &crossLinkIndices);
      
	if (ret)
	  {
	    // At least one cross-link involves the monomer currently
	    // iterated in the oligomer being analysed.

	    int index = 0;
	  
	    foreach(index, crossLinkIndices)
	      {
		CrossLink *crossLink = crossLinkList.at(index);

		// 	      qDebug() << __FILE__ << __LINE__ 
		// 			<< crossLink->name();
	      
		// First off, we can add the cross-link to the list of
		// cross-links of the oligomer(we'll need them to be
		// able to perform mass calculations). Note that this is
		// only copying the pointer to the actual cross-link in
		// the polymer's list of cross-links. Note also that a
		// cross-link might not be found more than once(the
		// call below first checks that the cross-link is not
		// already in the list).

		if (!oligomer->addCrossLink(crossLink))
		  {
		    // 		  qDebug() << __FILE__ << __LINE__ 
		    // 			    << "The cross-link:" 
		    // 			    << crossLink->name() 
		    // 			    << "was already in the"
		    // 			    << oligomer 
		    // 			    << "oligomer's list of cross-links: " 
		    // 		    "not duplicated.";
		  }
		else
		  {
		    // 		  qDebug() << __FILE__ << __LINE__ 
		    // 			    << "The cross-link:" 
		    // 			    << crossLink->name() 
		    // 			    << "was added to the"
		    // 			    << oligomer 
		    // 			    << "oligomer's list of cross-links.";
		  }

		const Monomer *iterMonomer = 0;
	      
		foreach(iterMonomer, *(crossLink->monomerList()))
		  {
		    // 		  qDebug() << __FILE__ << __LINE__ 
		    // 			    << iterMonomer->name();
		  
		    int foundIndex = 0;

		    Oligomer *foundOligomer = 
		      oligomerList->findOligomerEncompassing(iterMonomer,
							      &foundIndex);
		  
		    if (foundOligomer)
		      {
			// 		      qDebug() << __FILE__ << __LINE__ 
			// 				<< foundOligomer->name() << foundIndex;

			// One oligomer in the original oligomer list
			// encompasses a monomer that seems to be
			// cross-linked to the 'monomer' being iterated
			// in in the currently analyzed oligomer. Move
			// that oligomer to the clearance list of
			// oligomer that will need to be further
			// analyzed later.
		            
			oligomerList->removeAt(foundIndex);
		      
			clearanceOligomerList.append(foundOligomer);

			// Update the name of the oligomer with the name
			// of the new foundOligomer.
		      
			QString name = QString("%1+%2")
			  .arg(oligomer->name())
			  .arg(foundOligomer->name());
			oligomer->setName(name);
		      }
		  }
	      }
	    // End of
	    // foreach(index, crossLinkIndices)
	  }
      }
  
    // At this point we have one oligomer which we know is cross-linked
    // at least once(with another oligomer or the cross-link is between
    // two or more monomers in the same oligomer, think cyan fluorescent
    // protein). If monomers in that same oligomer were cross-linked to
    // other monomers in other oligomers, then these oligomers should by
    // now have been moved from the original list of oligomers
    //(oligomerList) to the clearance list of oligomers
    //(clearanceOligomerList). We have to iterate in each oligomer of that
    // clearance list and for each of its monomers, check if it has a
    // cross-link to any oligomer still in the original oligomerList
    //(this is what I call "pull-down" stuff). Found oligomers are
    // appended to the clearanceOligomerList.

    while (clearanceOligomerList.size())
      {
	Oligomer *iterOligomer = clearanceOligomerList.first();
      
	for(int iter = iterOligomer->startIndex();
	     iter < iterOligomer->endIndex() + 1; ++iter)
	  {
	    const Monomer *monomer = mp_polymer->at(iter);

	    // 	  qDebug() << __FILE__ << __LINE__ 
	    // 		    << monomer->name();

	    // What crossLinks do involve this monomer ?

	    QList<int> crossLinkIndices;
      
	    int ret = 
	      crossLinkList.crossLinksInvolvingMonomer(monomer,
							&crossLinkIndices);
      
	    if (ret)
	      {
		// At least one cross-link involves the monomer currently
		// iterated in the iterOligomer being analysed.

		int index = 0;
	  
		foreach(index, crossLinkIndices)
		  {
		    CrossLink *crossLink = crossLinkList.at(index);

		    // 		  qDebug() << __FILE__ << __LINE__ 
		    // 			    << crossLink->name();

		    // First off, we can add the cross-link to the list of
		    // cross-links of the oligomer(we'll need them to be
		    // able to perform mass calculations). Note that this is
		    // only copying the pointer to the actual cross-link in
		    // the polymer's list of cross-links. Note also that a
		    // cross-link might not be found more than once(the
		    // call below first checks that the cross-link is not
		    // already in the list).

		    if (!oligomer->addCrossLink(crossLink))
		      {
			// 		      qDebug() << __FILE__ << __LINE__ 
			// 				<< "The cross-link:" 
			// 				<< crossLink->name() 
			// 				<< "was already in the"
			// 				<< oligomer 
			// 				<< "oligomer's list of cross-links: " 
			// 			"not duplicated.";
		      }
		    else
		      {
			// 		      qDebug() << __FILE__ << __LINE__ 
			// 				<< "The cross-link:" 
			// 				<< crossLink->name() 
			// 				<< "was added to the"
			// 				<< oligomer 
			// 				<< "oligomer's list of cross-links.";
		      }
		  
		    const Monomer *iterMonomer = 0;

		    foreach(iterMonomer, *(crossLink->monomerList()))
		      {
			// 		      qDebug() << __FILE__ << __LINE__ 
			// 				<< iterMonomer->name();
		      
			int foundIndex = 0;
		      
			Oligomer *foundOligomer = 
			  oligomerList->findOligomerEncompassing(iterMonomer,
								  &foundIndex);
		  
			if (foundOligomer)
			  {
			    // 			  qDebug() << __FILE__ << __LINE__ 
			    // 				    << foundOligomer->name() << foundIndex;
			  
			    // One oligomer in the original oligomer list
			    // encompasses a monomer that seems to be
			    // cross-linked to the 'monomer' being iterated
			    // in in the currently analyzed oligomer. Move
			    // that oligomer to the clearance list of
			    // oligomer that will need to be further
			    // analyzed later.
		            
			    oligomerList->removeAt(foundIndex);
		      
			    clearanceOligomerList.append(foundOligomer);
			  
			    // Update the name of the oligomer with the name
			    // of the new foundOligomer.
			  
			    QString name = QString("%1+%2")
			      .arg(oligomer->name())
			      .arg(foundOligomer->name());
			    oligomer->setName(name);
			  }
		      }
		  }
		// End of
		// foreach(index, crossLinkIndices)
	      }
	    // End of(ret) ie cross-links involved monomer
	  }
	// End of
	//   for (int iter = iterOligomer->startIndex();
	//   iter < iterOligomer->endIndex() + 1; ++iter)

	// At this point this quarantinized oligomer might be removed
	// from the clearance clearanceOligomerList and its coordinates
	// be appended to the 'oligomer' list of coordinates. Then, the
	// quanrantinized oligomer might be destroyed.

	clearanceOligomerList.removeFirst();
      
	oligomer->appendCoordinates(iterOligomer);

	delete iterOligomer;
      }
  
    // At this point, all the oligomers in the clearance oligomer list
    // have all been dealt with, return the number of cross-linked
    // oligomers in this oligomer.

    //   for (int iter = 0; iter < oligomer->crossLinkList()->size(); ++iter)
    //     qDebug() << __FILE__ << __LINE__ 
    // 	      << "Finished for oligomer:" << 
    //       oligomer->crossLinkList()->at(iter)->name();
  
  
    return static_cast<CoordinateList *>(oligomer)->size();
  }


  void 
  Cleaver::emptyOligomerList()
  {
    while (!mp_oligomerList->isEmpty())
      {
	delete mp_oligomerList->first();
      }
  }

} // namespace massXpert