gain-state.h

/**************************************************************************
 * Copyright (c) 2001-2011 T. M. Murali                                   *
 * Copyright (c) 2011 Phillip Whisenhunt                                  *
 * Copyright (c) 2011 David Badger                                        *
 * Copyright (c) 2010 Jacqueline Addesa                                   *
 *                                                                        *
 * This file is part of Biorithm.                                         *
 *                                                                        *
 * Biorithm is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU General Public License as published by   *
 * the Free Software Foundation, either version 3 of the License, or      *
 * (at your option) any later version.                                    *
 *                                                                        *
 * Biorithm 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 Biorithm.  If not, see <http://www.gnu.org/licenses/>.      *
 *                                                                        *
 **************************************************************************/

// Purpose: classes to represent the state of a node in GAIN.

#ifndef _GAIN_STATE_H
#define _GAIN_STATE_H

#include <fstream>
#include <map>
using namespace std;

#include "old-annotations.h"

// enum MyGainTriState { NOT_ANNOTATED_STATE = -1, HYPOTHETICAL_STATE,
//                              ANNOTATED_STATE};
typedef MyGainAnnotationType MyGainTriStateType;
// allow only two states. i must distinguish MyGainBiStateType from
// MyGainTriStateType because the template specialisations
// MyGainState< MyGainBiStateType >::collapse() and MyGainState<
// MyGainTriStateType >::collapse() must do different things.
enum MyGainBiStateType { BI_HYPOTHETICAL_STATE = 0, BI_ANNOTATED_STATE};
//typedef MyGainTriStateType MyGainBiStateType;
// class MyGainBiStateType : public MyGainTriStateType
// {};

// how much uglier can my code get? i need this conversion in MyOneVersusNoneGainAlgo::computeCrossValidationResults().
inline MyGainTriStateType convertState(MyGainBiStateType in)
{
  if (BI_HYPOTHETICAL_STATE == in)
    return(HYPOTHETICAL_STATE);
  if (BI_ANNOTATED_STATE == in)
    return(ANNOTATED_STATE);
}


// class hierarchy of states.
template < typename StateType >
class MyGainState
{
private:
  // the weight associated with each state.
  map< StateType, MyNT > _stateWeights;
public:
  MyGainState()
      : _stateWeights()
    {}
  MyGainState(StateType state, MyNT w = 1)
      : _stateWeights()
    {
      //_stateWeights[state] = 1;
      _stateWeights[state] = w;
      this->_sanityCheck();
    }
  virtual ~MyGainState()
    {}

  virtual const MyGainState &operator+=(const MyGainState &rhs)
    {
      typename map< StateType, MyNT >::const_iterator itr;
      for (itr = rhs._stateWeights.begin(); itr != rhs._stateWeights.end(); itr++)
        _stateWeights[itr->first] += itr->second;
      this->_sanityCheck();
      return(*this);
    }

  virtual const MyGainState &operator*(MyNT weight)
    {
      typename map< StateType, MyNT >::const_iterator itr;
      for (itr = _stateWeights.begin(); itr != _stateWeights.end(); itr++)
        _stateWeights[itr->first] *= weight;
      this->_sanityCheck();
      return(*this);
    } 

  virtual bool operator==(const MyGainState &rhs) const
    {
      return(isEqual(rhs));
    }

  virtual bool isEqual(const MyGainState &rhs) const
    {
      typename map< StateType, MyNT >::const_iterator itr;
      for (itr = rhs._stateWeights.begin(); itr != rhs._stateWeights.end(); itr++)
        if (_stateWeights.find(itr->first)->second != itr->second)
          return(false);
      return(true);
    }

  virtual bool isEqual(const MyGainState &rhs, MyNT epsilon) const
    {
      typename map< StateType, MyNT >::const_iterator itr;
      for (itr = rhs._stateWeights.begin(); itr != rhs._stateWeights.end(); itr++)
        if (fabs(_stateWeights.find(itr->first)->second - itr->second) > epsilon)
          return(false);
      return(true);
    }

  MyNT getWeight(StateType state) const
    {
      typename map< StateType, MyNT >::const_iterator itr;
      itr = _stateWeights.find(state);
      if (_stateWeights.end() == itr)
        return(0);
      return(itr->second);
    }
  
  StateType collapse(MyNT threshold = 0);
  
  ostream &print(ostream &str) const
    {
      typename map< StateType, MyNT >::const_iterator itr;
      for (itr = _stateWeights.begin(); itr != _stateWeights.end(); itr++)
        str << itr->first << ": " << itr->second << "\t";
//      str << endl;
      return(str);
    }
    
  void _sanityCheck() const
    {
      typename map< StateType, MyNT >::const_iterator itr;
      for (itr = _stateWeights.begin(); itr != _stateWeights.end(); itr++)
        {
          if (itr->second < 0)
            cerr << "MyGainState::_sanityCheck: negative weight for state " << itr->first << endl;
        }
    }
};

template< typename StateType >
class MyGainStateInfo
{

private:

  string id;
  // initial state.
  StateType initialState;
  // state in the current iteration.
  StateType currentState;
  // state in the last iteration.
  StateType previousState;
  // what is the distinction between all these states? realState is
  // the *real* state of the orf. it will never change. initialState
  // is the initial state of the orf, at the beginning of
  // runGainNet(). for hypothetical proteins, realState always
  // equals initialState. the equality holds for all other proteins
  // too, except for those that are being cross-validated. for them
  // initialState will be HYPOTHETICAL_STATE. currentState is the
  // current state in runGainNet() and previousState is the
  // previous value of currentState (the value before currentState
  // changed or the value of initialState).
  //StateType realState;
  
  bool clamped;

//   // input is the weighted sum of neighbour states.
//   MyNT currentInput;
//   MyNT previousInput;
  
//   // output is contribution to neighbours' input from this org. input
//   // and output may not be the same when i am dividing by the degree.
//   MyNT currentOutput;
//   MyNT previousOutput;
  
public:

  string getId() const
  {
    return(id);
  }

  void setId(string _id)
  {
    id = _id;
  }

  StateType getInitialState() const
  {
    return(initialState);
  }
  
  StateType getCurrentState() const
  {
    return(currentState);
  }
  
  StateType getPreviousState() const
  {
    return(previousState);
  }
  
  /*StateType getRealState() const
  {
    return(realState);
  }*/
//   // given an input and a theshold, compute the corr. state.
//   static StateType computeState(MyNT input, MyNT threshold)
//     {
//       if (input < threshold)
//         return(NOT_ANNOTATED_STATE);
//       // this policy is better than not checking if input and
//       // threshold are equal since it makes solitary nodes and nodes
//       // all of whose neighbours are in HYPOTHETICAL_STATE also to be
//       // in HYPOTHETICAL_STATE.
//       if (input == threshold)
//         return(HYPOTHETICAL_STATE);
      
//       // i don't allow a node to go into HYPOTHETICAL_STATE. it can be
//       // in that state only at the beginning. NOT A GOOD POLICY. see
//       // the comment above.
//       return(ANNOTATED_STATE);
//     }
  
  // compute the current state given the current input.
//   StateType computeCurrentState(MyNT threshold) const
//     {
//       return(computeState(currentInput, threshold));
//     }
  
  
//   // flip state from ANNOTATED_STATE to NOT_ANNOTATED_STATE or vice-versa.
//   void flipState()
//     {
//       previousState = currentState;
//       if (ANNOTATED_STATE == currentState)
//         currentState = NOT_ANNOTATED_STATE;
//       else 
//         currentState = ANNOTATED_STATE;
//     }
  
  StateType getState() const
  {
    return(currentState);
  }


  // set the state to be newState.
  void setState(StateType newState)
    {
      previousState = currentState;
      currentState = newState;
    }

  void initialiseState(StateType state)
  {
    //realState = initialState = state;
    initialState = state;
    currentState = previousState = initialState;
#ifdef DEBUG
//    state.print(cout);
#endif// DEBUG    
  }
  
  /*void setRealState(StateType state)
  {
        realState = state;
  }*/

  
  // update the state based on the threshold (if need be).
//   void updateState(MyNT threshold)
//     {
//       previousState = currentState;
//       currentState = computeCurrentState(threshold);
//     }
  
//   // a node has flipped if the sign of (energy - threshold) is
//   // different from the the sign of currentState, i.e., if the product
//   // of these two numbers is negative.
//   bool checkIfFlipped(MyNT threshold) const
//     {
//       // if currentInput equals threshold, then the rhs is always 0 so
//       // the test returns false, which is incorrect, since i want
//       // ANNOTATED_STATE $\equiv$ (currentInput >= threshold).
//       //
//       if (currentInput >= threshold)
//         // new state will be ANNOTATED_STATE so i return true if the
//         // currentState is NOT_ANNOTATED_STATE.
//         return(NOT_ANNOTATED_STATE == currentState);
//       // currentInput < threshold, so i return true if the current
//       // state is ANNOTATED_STATE.
//       return(ANNOTATED_STATE == currentState);
//     }
  

  
  // in the version of the algorithm where the node can have any
  // state, it is better to check if the node has a "changed" state
  // rather than a "flipped" state.
//   bool checkIfChanged(MyNT threshold, StateType *newState = NULL) const
//     {
//       StateType state = computeCurrentState(threshold);
//       if (newState)
//         *newState = state;
//       return(state != currentState);
//     }

  bool hasConverged()
    {
      return(currentState.collapse() == previousState.collapse());
    }
  
  
  bool hasConvergedApproximately(MyNT epsilon = 1e-03) const
    {
      return(currentState.isEqual(previousState, epsilon));
    }
  
  
  bool isClamped() const
  {
    return(clamped);
  }

  void clamp()
  {
    clamped = true;
  }

  void unclamp()
  {
    clamped = false;
  }
  
};

// each value of MyGainBiState is a map whose keys can take the two
// legal values of MyGainBiStateType.
typedef MyGainState< MyGainBiStateType > MyGainBiState;
// notice that the template parameter is the type i just defined. 
typedef MyGainStateInfo< MyGainBiState > MyGainBiStateInfo;

// each value of MyGainTriState is a map whose keys can take the three
// legal values of MyGainTriStateType.
typedef MyGainState< MyGainTriStateType > MyGainTriState;
// notice that the template parameter is the type i just defined. 
typedef MyGainStateInfo< MyGainTriState > MyGainTriStateInfo;

// template specialisation.

#endif // _GAIN-STATE_H