libmems-1.6-dev 1.6.0+4725-4 / usr / include / libMems-1.6 / libMems / HomologyHMM / homology.h

/* Code generated by HMMoC version VERSION, Copyright (C) 2006 Gerton Lunter */
/* Generated from file homology.xml (author: Aaron Darling) on Mon Jul 16 11:09:12 EST 2007 */

/*
This file is a work based on HMMoC VERSION, a hidden Markov model compiler.
Copyright (C) 2006 by Gerton Lunter, Oxford University.

HMMoC and works based on it are 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 2 of the License, or (at your option) any later version.

HMMOC 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 HMMoC; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
*/

#ifndef _homology_h_
#define _homology_h_


#include "dptables.h"
#include "algebras.h"

#include <vector>
#include <iostream>
#include <string>
#include <algorithm>

using std::vector;
using std::cout;

#include <map>

using std::map;

struct Params;

//void run(std::string& sequence, std::string& prediction, double goHomologous = 0.004, double goUnrelated = 0.004, std::vector<double>* emitHomologous = NULL, std::vector<double>* emitUnrelated = NULL);

void run(std::string& sequence, std::string& prediction, const Params& params );

// Here go the state memory clique typedefs:
typedef States<bfloat,2> Statesblock2;
typedef States<bfloat,1> Statesblock1;
typedef States<bfloat,1> Statesblock3;

class HomologyDPTable {
    public:
    // If true, this class' destructor will delete the DP arrays
    bool isInCharge;
    // Pointers to arrays containing ids of states and transitions
    const std::string* const stateId;
    const std::string* const emissionId;
    const std::string* const transitionId;
    const std::string* const transitionFrom;
    const std::string* const transitionTo;
    const std::string* const transitionProb;
    const std::string* const transitionEmit;
    const std::string* const outputId;
    // The actual DP tables, and total sequence lengths (which determine size of DP arrays) follow:
    int iLen;
    DPTable<Statesblock2,1> StateMemoryblock2;
    DPTable<Statesblock1,0> StateMemoryblock1;
    DPTable<Statesblock3,0> StateMemoryblock3;
    // Member functions:
    public:
    // Default copy constructor is used; user has to set isInCharge appropriately afterwards!
    HomologyDPTable(int iLen);
    ~HomologyDPTable();
    // returns probability from DP table, given position and int or std::string state identifier
    bfloat getProb(int iState ,int ) const;
    bfloat getProb(const std::string sState ,int ) const;
    // converts std::string identifier (for state, transition or emission) into integer id
    static int getId(const std::string& sState);
    static const std::string& getTransitionId(int id);
    static const std::string& getEmissionId(int id);
    static const std::string& getStateId(int id);
    static const std::string& getOutputId(int id);
    static void _cleanup() { getId("_cleanup_"); }
};

// give a name to the real type used for this HMM
typedef bfloat HomologyReal;
// define type for a 'short' real -- usually double, but can be logspace for efficiency
typedef double HomologyShortReal;



class HomologyBaumWelch {
    public:
    // Default copy constructor is used.
    // Void constructor:
    HomologyBaumWelch() { resetCounts(); }
    // Not calling resetCounts() across calls allows to aggregate results over multiple datasets
    void resetCounts();
    void scaleCounts(bfloat scale);
    // Translate an identifier (string or integer) to the index into their corresponding Baum-Welch counter array (below)
    // Which array is used for any particular emission/transition depends on its order signature - see documentation for details
    int transitionIndex(int intId) const { return atransitionIdx[intId]; }
    int transitionIndex(string strId) const;
    int emissionIndex(int intId) const { return aemissionIdx[intId]; }
    int emissionIndex(string strId) const;
    // Now follow, in triplets (one for each order signature):
    //  Transition or emission counters;
    //  Array of identifiers; and
    //  Dimension of array (number of counters).
    bfloat transitionBaumWelchCount0[8];
    static int transitionIdentifier0[8];   
    static const int transitionDimension0 = 8;
    bfloat emissionBaumWelchCount0[1];
    static int emissionIdentifier0[1];   
    static const int emissionDimension0 = 1;
    bfloat emissionBaumWelchCount1[8][2];
    static int emissionIdentifier1[2];   
    static const int emissionDimension1 = 2;
    private:
    static int atransitionIdx[8];
    static int aemissionIdx[3];
    static map<const std::string,int> mId;
};



class Path {
    //protected:
    public:
    vector<int> transitions;
    public:
    unsigned int size() const {                             // Number of transitions in path
        return transitions.size(); 
    }
    int transition(int i) const {                           // i-th transition
        return transitions[i]; 
    }
    virtual double prob(int) const = 0;                  // i-th transition*emission probability:
    virtual const vector<int>& emission(int) const = 0;  // i-th emission vector
    virtual int fromState(int) const = 0;                // State at from-end of i-th transition
    virtual int toState(int) const = 0;                  // State at to-end of i-th transition
    virtual int nextFrom(int) const = 0;                 // index of next sibling, -1 if no more (always -1 for simple paths)
    virtual int nextTo(int) const = 0;                   // index of child (always i+1 for simple paths), or -1 if no more
    virtual ~Path() {} 
};

ostream& operator<<(ostream& os, const Path& p);

class SimplePath: public Path {
    public:
    vector<double> probs;
    vector<vector<int> > emissions;
    vector<int> froms;
    vector<int> tos;
    public:
    void addEdge(int transition, double prob, vector<int>& emission, int from, int to);
    double prob(int index) const;
    int nextFrom(int index) const;
    int nextTo(int index) const;
    const vector<int>& emission(int index) const;
    int fromState(int index) const;
    int toState(int index) const;
    void reverse();
};
struct Params {
    double iStartHomologous;
    double iGoHomologous;
    double iGoUnrelated;
    double iGoStopFromUnrelated;
    double iGoStopFromHomologous;
    double aEmitHomologous[8];
    double aEmitUnrelated[8];
};

bfloat Forward(HomologyDPTable** ppOutTable,Params iPar,char *aSeq,int iLen);

bfloat Backward(HomologyBaumWelch& bw,HomologyDPTable* pInTable,HomologyDPTable** ppOutTable,Params iPar,char *aSeq,int iLen);

bfloat Viterbi_recurse(HomologyDPTable** ppOutTable,Params iPar,char *aSeq,int iLen);

Path& Viterbi_trace(HomologyDPTable* pInTable,Params iPar,char *aSeq,int iLen);

#endif // _homology_h_

/* --- end of HMMoC-generated file --- */