/* Code generated by HMMoC version VERSION, Copyright (C) 2006 Gerton Lunter */
/* Generated from file readalign.xml (author:  Gerton Lunter ) on Mon Dec 20 11:41:40 GMT 2010 */

/*
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
*/


#include "readalign.h"

#include "banding.h"

const extern string _AlignstateId[];
const extern string _AlignemissionId[];
const extern string _AligntransitionId[];
const extern string _AligntransF[];
const extern string _AligntransT[];
const extern string _AligntransP[];
const extern string _AligntransE[];
const extern string _AlignoutputId[];
const extern string _Alignempty;
const extern int _AlignstateNum;
const extern int _AlignemitNum;
const extern int _AligntransNum;
const extern int _AlignoutputNum;

AlignDPTable::AlignDPTable(int iLen1,int iLen2) : isInCharge(true), stateId(_AlignstateId), emissionId(_AlignemissionId), transitionId(_AligntransitionId), transitionFrom(_AligntransF), transitionTo(_AligntransT), transitionProb(_AligntransP), transitionEmit(_AligntransE), outputId(_AlignoutputId) {
    // init code:
    this->iLen1 = iLen1;
    this->iLen2 = iLen2;
    StateMemoryblock4.allocate(1+iLen2);
    StateMemoryblock5.allocate();
    StateMemoryblock3.allocate(0+iLen1, 1+iLen2);
    StateMemoryblock2.allocate(1+iLen2);
    StateMemoryblock1.allocate();
}


AlignDPTable::~AlignDPTable() {
    if (!isInCharge) {
        // make sure data does not get deleted:
        StateMemoryblock4.absolve();
        StateMemoryblock5.absolve();
        StateMemoryblock3.absolve();
        StateMemoryblock2.absolve();
        StateMemoryblock1.absolve();
    } // if(!isInCharge)
} // destructor

const string& AlignDPTable::getTransitionId(int id) { return id>=0 && id<_AligntransNum ? _AligntransitionId[id] : _Alignempty; }
const string& AlignDPTable::getEmissionId(int id) { return id>=0 && id<_AlignemitNum ? _AlignemissionId[id] : _Alignempty; }
const string& AlignDPTable::getStateId(int id) { return id>=0 && id<_AlignstateNum ? _AlignstateId[id] : _Alignempty; }
const string& AlignDPTable::getOutputId(int id) { return id>=0 && id<_AlignoutputNum ? _AlignoutputId[id] : _Alignempty; }
int AlignDPTable::getId(const string& sId)
{
    static bool bInit = false;
    static map<string,int>* pmId;
    if (!bInit) {
        pmId = new map<string,int>();
        for (int i=0;i<_AlignstateNum;i++) {
            (*pmId)[_AlignstateId[i]] = i;         // add state identifiers
        }
        for (int i=0; i<_AlignemitNum; i++) {
            (*pmId)[_AlignemissionId[i]] = i;      // add emission identifiers
        }
        for (int i=0; i<_AligntransNum; i++) {  
            (*pmId)[_AligntransitionId[i]] = i;    // add transition identifiers
        }
        for (int i=0; i<_AlignoutputNum; i++) {
            (*pmId)[_AlignoutputId[i]] = i;        // finally, add output identifiers
        }
        bInit = true;
    }
    map<string,int>::iterator iter = pmId->find(sId);
    if (iter == pmId->end()) {
        if (sId == "_cleanup_") {
            delete pmId;
            } else {
            cout << "AlignDPTable::getId: WARNING: identifier '" << sId << "' not found." << endl;
        }
        return -1;
    }
    return iter->second;
}


logspace AlignDPTable::getProb(const string sState ,int iPos0,int iPos1) const
{
    return getProb(getId(sState) ,iPos0,iPos1);
}


logspace AlignDPTable::getProb(int iState ,int iPos0,int iPos1) const
{
    const logspace *CurStateMemoryblock1Secondary;
    const logspace *CurStateMemoryblock2Secondary;
    const logspace *CurStateMemoryblock3Secondary;
    const logspace *CurStateMemoryblock4Secondary;
    const logspace *CurStateMemoryblock5Secondary;
    static const int blockTable[] = {0, 1, 2, 2, 2, 3, 4};
    static const int stateTable[] = {0, 0, 0, 1, 2, 0, 0};
    switch (blockTable[iState]) {
        default:
        return 0.0;
        break;
        case 0:
        if ((iPos0+0>=0)&&(iPos0+0<=0)&&(iPos1+0>=0)&&(iPos1+0<=0)) {
            CurStateMemoryblock1Secondary = this->StateMemoryblock1.read();
            return CurStateMemoryblock1Secondary[stateTable[iState]];
        } 
        else { 
            return 0.0;
            
        }
        break;
        case 1:
        if ((iPos0+0>=0)&&(iPos0+0<=0)&&(iPos1+0>=0)&&(iPos1+0<=iLen2+0)) {
            CurStateMemoryblock2Secondary = this->StateMemoryblock2.read((iPos1-(0))-(0));
            return CurStateMemoryblock2Secondary[stateTable[iState]];
        } 
        else { 
            return 0.0;
            
        }
        break;
        case 2:
        if ((iPos0+0>=1)&&(iPos0+0<=iLen1+0)&&(iPos1+0>=0)&&(iPos1+0<=iLen2+0)) {
            CurStateMemoryblock3Secondary = this->StateMemoryblock3.read((iPos0-(0))-(1), (iPos1-(0))-(0));
            return CurStateMemoryblock3Secondary[stateTable[iState]];
        } 
        else { 
            return 0.0;
            
        }
        break;
        case 3:
        if ((iPos0+0>=iLen1+0)&&(iPos0+0<=iLen1+0)&&(iPos1+0>=0)&&(iPos1+0<=iLen2+0)) {
            CurStateMemoryblock4Secondary = this->StateMemoryblock4.read((iPos1-(0))-(0));
            return CurStateMemoryblock4Secondary[stateTable[iState]];
        } 
        else { 
            return 0.0;
            
        }
        break;
        case 4:
        if ((iPos0+0>=iLen1+0)&&(iPos0+0<=iLen1+0)&&(iPos1+0>=iLen2+0)&&(iPos1+0<=iLen2+0)) {
            CurStateMemoryblock5Secondary = this->StateMemoryblock5.read();
            return CurStateMemoryblock5Secondary[stateTable[iState]];
        } 
        else { 
            return 0.0;
            
        }
    } // switch
} // DPTable...::getProb(int,...)

const extern string _AlignBandingstateId[];
const extern string _AlignBandingemissionId[];
const extern string _AlignBandingtransitionId[];
const extern string _AlignBandingtransF[];
const extern string _AlignBandingtransT[];
const extern string _AlignBandingtransP[];
const extern string _AlignBandingtransE[];
const extern string _AlignBandingoutputId[];
const extern string _AlignBandingempty;
const extern int _AlignBandingstateNum;
const extern int _AlignBandingemitNum;
const extern int _AlignBandingtransNum;
const extern int _AlignBandingoutputNum;

AlignBandingDPTable::AlignBandingDPTable(int iLen1,int iLen2) : isInCharge(true), stateId(_AlignBandingstateId), emissionId(_AlignBandingemissionId), transitionId(_AlignBandingtransitionId), transitionFrom(_AlignBandingtransF), transitionTo(_AlignBandingtransT), transitionProb(_AlignBandingtransP), transitionEmit(_AlignBandingtransE), outputId(_AlignBandingoutputId) {
    // init code:
    this->iLen1 = iLen1;
    this->iLen2 = iLen2;
    StateMemoryblock4.allocate(1+iLen2);
    StateMemoryblock5.allocate();
    StateMemoryblock3withbanding.allocate(0+iLen1, 1+iLen2);
    StateMemoryblock2.allocate(1+iLen2);
    StateMemoryblock1.allocate();
}


AlignBandingDPTable::~AlignBandingDPTable() {
    if (!isInCharge) {
        // make sure data does not get deleted:
        StateMemoryblock4.absolve();
        StateMemoryblock5.absolve();
        StateMemoryblock3withbanding.absolve();
        StateMemoryblock2.absolve();
        StateMemoryblock1.absolve();
    } // if(!isInCharge)
} // destructor

const string& AlignBandingDPTable::getTransitionId(int id) { return id>=0 && id<_AlignBandingtransNum ? _AlignBandingtransitionId[id] : _AlignBandingempty; }
const string& AlignBandingDPTable::getEmissionId(int id) { return id>=0 && id<_AlignBandingemitNum ? _AlignBandingemissionId[id] : _AlignBandingempty; }
const string& AlignBandingDPTable::getStateId(int id) { return id>=0 && id<_AlignBandingstateNum ? _AlignBandingstateId[id] : _AlignBandingempty; }
const string& AlignBandingDPTable::getOutputId(int id) { return id>=0 && id<_AlignBandingoutputNum ? _AlignBandingoutputId[id] : _AlignBandingempty; }
int AlignBandingDPTable::getId(const string& sId)
{
    static bool bInit = false;
    static map<string,int>* pmId;
    if (!bInit) {
        pmId = new map<string,int>();
        for (int i=0;i<_AlignBandingstateNum;i++) {
            (*pmId)[_AlignBandingstateId[i]] = i;         // add state identifiers
        }
        for (int i=0; i<_AlignBandingemitNum; i++) {
            (*pmId)[_AlignBandingemissionId[i]] = i;      // add emission identifiers
        }
        for (int i=0; i<_AlignBandingtransNum; i++) {  
            (*pmId)[_AlignBandingtransitionId[i]] = i;    // add transition identifiers
        }
        for (int i=0; i<_AlignBandingoutputNum; i++) {
            (*pmId)[_AlignBandingoutputId[i]] = i;        // finally, add output identifiers
        }
        bInit = true;
    }
    map<string,int>::iterator iter = pmId->find(sId);
    if (iter == pmId->end()) {
        if (sId == "_cleanup_") {
            delete pmId;
            } else {
            cout << "AlignBandingDPTable::getId: WARNING: identifier '" << sId << "' not found." << endl;
        }
        return -1;
    }
    return iter->second;
}


logspace AlignBandingDPTable::getProb(const string sState ,int iPos0,int iPos1) const
{
    return getProb(getId(sState) ,iPos0,iPos1);
}


logspace AlignBandingDPTable::getProb(int iState ,int iPos0,int iPos1) const
{
    const logspace *CurStateMemoryblock1Secondary;
    const logspace *CurStateMemoryblock2Secondary;
    const logspace *CurStateMemoryblock3withbandingSecondary;
    const logspace *CurStateMemoryblock4Secondary;
    const logspace *CurStateMemoryblock5Secondary;
    static const int blockTable[] = {0, 1, 2, 2, 2, 3, 4};
    static const int stateTable[] = {0, 0, 0, 1, 2, 0, 0};
    switch (blockTable[iState]) {
        default:
        return 0.0;
        break;
        case 0:
        if ((iPos0+0>=0)&&(iPos0+0<=0)&&(iPos1+0>=0)&&(iPos1+0<=0)) {
            CurStateMemoryblock1Secondary = this->StateMemoryblock1.read();
            return CurStateMemoryblock1Secondary[stateTable[iState]];
        } 
        else { 
            return 0.0;
            
        }
        break;
        case 1:
        if ((iPos0+0>=0)&&(iPos0+0<=0)&&(iPos1+0>=0)&&(iPos1+0<=iLen2+0)) {
            CurStateMemoryblock2Secondary = this->StateMemoryblock2.read((iPos1-(0))-(0));
            return CurStateMemoryblock2Secondary[stateTable[iState]];
        } 
        else { 
            return 0.0;
            
        }
        break;
        case 2:
        if ((iPos0+0>=1)&&(iPos0+0<=iLen1+0)&&(iPos1+0>=0)&&(iPos1+0<=iLen2+0)) {
            CurStateMemoryblock3withbandingSecondary = this->StateMemoryblock3withbanding.read((iPos0-(0))-(1), (iPos1-(0))-(0));
            return CurStateMemoryblock3withbandingSecondary[stateTable[iState]];
        } 
        else { 
            return 0.0;
            
        }
        break;
        case 3:
        if ((iPos0+0>=iLen1+0)&&(iPos0+0<=iLen1+0)&&(iPos1+0>=0)&&(iPos1+0<=iLen2+0)) {
            CurStateMemoryblock4Secondary = this->StateMemoryblock4.read((iPos1-(0))-(0));
            return CurStateMemoryblock4Secondary[stateTable[iState]];
        } 
        else { 
            return 0.0;
            
        }
        break;
        case 4:
        if ((iPos0+0>=iLen1+0)&&(iPos0+0<=iLen1+0)&&(iPos1+0>=iLen2+0)&&(iPos1+0<=iLen2+0)) {
            CurStateMemoryblock5Secondary = this->StateMemoryblock5.read();
            return CurStateMemoryblock5Secondary[stateTable[iState]];
        } 
        else { 
            return 0.0;
            
        }
    } // switch
} // DPTable...::getProb(int,...)

logspace hmmocMax(logspace i, logspace j) { return i>j ? i : j; }
void hmmocMaxInPlace(logspace& i, logspace j) { if (i<j) i=j; }
ostream& operator<<(ostream& os, const Path& p)
{
    for (unsigned int i=0; i<p.size(); i++) {
        os << p.fromState(i) << "--{";
            for (unsigned int j=0; j<p.emission(i).size(); j++) {
                if (j) os<<",";
                os<<p.emission(i)[j];
            }
        os<<"};"<<p.prob(i)<<"-->"<<p.toState(i)<<endl;
    }
    return os;
}

void SimplePath::addEdge(int tr, double p, vector<int>& e, int f, int t) {
    transitions.push_back(tr);
    probs.push_back(p);
    emissions.push_back(e);
    froms.push_back(f);
    tos.push_back(t);
}

void SimplePath::reverse() 
{
    std::reverse(transitions.begin(),transitions.end());
    std::reverse(probs.begin(),probs.end());
    std::reverse(emissions.begin(),emissions.end());
    std::reverse(froms.begin(),froms.end());
    std::reverse(tos.begin(),tos.end());
}

double SimplePath::prob(int i) const {
    return probs[i];
}

int SimplePath::nextFrom(int i) const {
    if (i+1 < (int)transitions.size())
    return i+1;
    else
    return -1;
}

int SimplePath::nextTo(int i) const {
    return -1;
}

const vector<int>& SimplePath::emission(int i) const {
    return emissions[i];
}

int SimplePath::fromState(int i) const {
    return froms[i];
}

int SimplePath::toState(int i) const {
    return tos[i];
}

const string _AlignstateId[] = {"start","pad1","delete","insert","match","pad2","end"};
const string _AlignemissionId[] = {"emit1","emit12","emit2free","empty"};
const string _AligntransitionId[] = {"trSS","trPP1","trPM","trPI","trMM","trMI","trMD","trMP","trIM","trII","trIE","trIP","trDM","trDD","trPP2","trPE"};
const string _AligntransF[] = {"start","pad1","pad1","pad1","match","match","match","match","insert","insert","insert","insert","delete","delete","pad2","pad2"};
const string _AligntransT[] = {"pad1","pad1","match","insert","match","insert","delete","pad2","match","insert","delete","pad2","match","delete","pad2","end"};
const string _AligntransP[] = {"probSkip","probSkip","probMatch_Begin","probGapOpen_Begin","probMatch","probGapOpen","probGapOpen","probSkip","probGapEnd","probGapExtend","probGapOpen","probSkip","probGapEnd","probGapExtend","probSkip","probSkip"};
const string _AligntransE[] = {"empty","emit2free","emit12","emit1","emit12","emit1","emit2free","empty","emit12","emit1","emit2free","empty","emit12","emit2free","emit2free","empty"};
const string _AlignoutputId[] = {"sequence1","sequence2"};
const string _Alignempty = "";
const int _AlignstateNum = 7;
const int _AlignemitNum = 4;
const int _AligntransNum = 16;
const int _AlignoutputNum = 2;




logspace Viterbi_recurse(AlignDPTable** ppOutTable,char* pQuality1,char* pSequence1,char* pSequence2,int iGapExtendPhred,int iGapOpenPhred,int iInsNucPrior,int iLen1,int iLen2,int iStartMean,int iStartSD) {
    logspace iTransition[16];
    logspace *CurStateMemoryblock4To;
    const logspace *CurStateMemoryblock4From;
    const logspace *CurStateMemoryblock5From;
    logspace *CurStateMemoryblock3To;
    const logspace *CurStateMemoryblock3From;
    logspace *CurStateMemoryblock2To;
    const logspace *CurStateMemoryblock2From;
    logspace *CurStateMemoryblock1To;
    const logspace *CurStateMemoryblock1From;
    int iPrevSlowCoord;
    /* initialization of various probabilities */
    
    static const double db = 10.0 / log(0.1);
    static const logspace iOne = 1.0;
    logspace iEquilibrium = exp( iInsNucPrior / db );     // likelihood of any unaligned read base
    logspace aQuality1[ iLen1 ];
    logspace aHalfQuality1[ iLen1 ];
    for (int i=0; i<iLen1; i++) {
        if (pQuality1[i] < (char)33) {
            cout << "Internal error: Negative quality at position "<<i<<endl;
        }
        aQuality1[i] = exp( (pQuality1[i]-33) / db );
        aHalfQuality1[i] = exp( (pQuality1[i]-33) / (2.0 * db) );
        if (pSequence1[i] == 'N') {
            aQuality1[i] = iOne;
            aHalfQuality1[i] = iOne;
        }
    }
    logspace iGapExtend = exp( iGapExtendPhred / db );
    logspace iGapOpen = exp( iGapOpenPhred / db ) * (iOne - iGapExtend);
    double iLocationFactor = -1.0 / (2.0*iStartSD*iStartSD);
    int iSymbol[2];
    if (false && iSymbol[0] == iSymbol[0]) {}   // avoid 'unused variable' warnings
    logspace iEmission[1];
    /* temporary storage for ordinary reals */
    register logspace iTempResult[1];
    /* temporary storage for extended-exponent reals */
    register logspace iTempProb[1];
    AlignDPTable dp(iLen1,iLen2);
    iTransition[0] = iOne;
    
    iTransition[1] = iOne;
    
    iTransition[4] = iOne;
    
    iTransition[5] = iGapOpen;
    
    iTransition[6] = iGapOpen;
    
    iTransition[7] = iOne;
    
    iTransition[8] = iOne;
    
    iTransition[9] = iGapExtend;
    
    iTransition[10] = iGapOpen;
    
    iTransition[11] = iOne;
    
    iTransition[12] = iOne;
    
    iTransition[13] = iGapExtend;
    
    iTransition[14] = iOne;
    
    iTransition[15] = iOne;
    dp.StateMemoryblock5.write()[0] = 1.0;
    dp.StateMemoryblock5.written();
    iPrevSlowCoord = -1;
    for (int iPos1=(iLen2+1)-1; iPos1>=0; --iPos1) {
        for (int iPos0=(iLen1+1)-1; iPos0>=0; --iPos0) {
            if ((iPos0+0>=iLen1+0)&&(iPos1+0>=iLen2+0)) {
            }
            if ((iPos0+0>=iLen1+0)) {
                if ((iPos1+0<=iLen2+-1)) {
                    iSymbol[0] = pSequence2[iPos1+0];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock4To = dp.StateMemoryblock4.write((iPos1-(0))-(0));
                iEmission[0] = iOne;
                if ((iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock4From = dp.StateMemoryblock4.read((iPos1-(-1))-(0));
                    CurStateMemoryblock4To[0] = ((iTransition[14])*(iEmission[0]))*CurStateMemoryblock4From[0];
                }
                iEmission[0] = iOne;
                if ((iPos1+0>=iLen2+0)) {
                    CurStateMemoryblock5From = dp.StateMemoryblock5.read();
                    hmmocMaxInPlace( CurStateMemoryblock4To[0], ((iTransition[15])*(iEmission[0]))*CurStateMemoryblock5From[0] );
                }
                dp.StateMemoryblock4.written();
            }
            if ((iPos0+0>=1)) {
                if ((iPos0+0<=iLen1+-1)) {
                    iSymbol[0] = pSequence1[iPos0+0];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                if ((iPos1+0<=iLen2+-1)) {
                    iSymbol[1] = pSequence2[iPos1+0];
                } 
                else { 
                    iSymbol[1] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock3To = dp.StateMemoryblock3.write((iPos0-(0))-(1), (iPos1-(0))-(0));
                if ((iSymbol[0] == iSymbol[1]) || (iSymbol[1] == 'N')) 
                iTempResult[0] = iOne;
                else 
                iTempResult[0] = aQuality1[ (iPos0)-0 ];
                iEmission[0] = iTempResult[0];
                if ((iPos0+1<=iLen1+0)&&(iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(-1))-(1), (iPos1-(-1))-(0));
                    CurStateMemoryblock3To[1] = ((iTransition[8])*(iEmission[0]))*CurStateMemoryblock3From[2];
                    CurStateMemoryblock3To[0] = ((iTransition[12])*(iEmission[0]))*CurStateMemoryblock3From[2];
                    CurStateMemoryblock3To[2] = ((iTransition[4])*(iEmission[0]))*CurStateMemoryblock3From[2];
                }
                iEmission[0] = iOne;
                if ((iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(0))-(1), (iPos1-(-1))-(0));
                    hmmocMaxInPlace( CurStateMemoryblock3To[1], ((iTransition[10])*(iEmission[0]))*CurStateMemoryblock3From[0] );
                    hmmocMaxInPlace( CurStateMemoryblock3To[0], ((iTransition[13])*(iEmission[0]))*CurStateMemoryblock3From[0] );
                    hmmocMaxInPlace( CurStateMemoryblock3To[2], ((iTransition[6])*(iEmission[0]))*CurStateMemoryblock3From[0] );
                }
                iEmission[0] = ((iEquilibrium > aHalfQuality1[ (iPos0)-0 ]) ? iEquilibrium : aHalfQuality1[ (iPos0)-0 ])  /* Was: iEquilibrium */;
                if ((iPos0+1<=iLen1+0)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(-1))-(1), (iPos1-(0))-(0));
                    hmmocMaxInPlace( CurStateMemoryblock3To[1], ((iTransition[9])*(iEmission[0]))*CurStateMemoryblock3From[1] );
                    hmmocMaxInPlace( CurStateMemoryblock3To[2], ((iTransition[5])*(iEmission[0]))*CurStateMemoryblock3From[1] );
                }
                iEmission[0] = iOne;
                if ((iPos0+0>=iLen1+0)) {
                    CurStateMemoryblock4From = dp.StateMemoryblock4.read((iPos1-(0))-(0));
                    hmmocMaxInPlace( CurStateMemoryblock3To[1], ((iTransition[11])*(iEmission[0]))*CurStateMemoryblock4From[0] );
                    hmmocMaxInPlace( CurStateMemoryblock3To[2], ((iTransition[7])*(iEmission[0]))*CurStateMemoryblock4From[0] );
                }
                dp.StateMemoryblock3.written();
            }
            if ((iPos0+0<=0)) {
                if ((iPos0+0<=iLen1+-1)) {
                    iSymbol[0] = pSequence1[iPos0+0];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                if ((iPos1+0<=iLen2+-1)) {
                    iSymbol[1] = pSequence2[iPos1+0];
                } 
                else { 
                    iSymbol[1] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock2To = dp.StateMemoryblock2.write((iPos1-(0))-(0));
                if ((iSymbol[0] == iSymbol[1]) || (iSymbol[1] == 'N')) 
                iTempResult[0] = iOne;
                else 
                iTempResult[0] = aQuality1[ (iPos0)-0 ];
                iEmission[0] = iTempResult[0];
                if ((iPos0+1<=iLen1+0)&&(iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(-1))-(1), (iPos1-(-1))-(0));
                    CurStateMemoryblock2To[0] = (((logspace::doubleexp( (((iPos1)-1 - iStartMean) * ((iPos1)-1 - iStartMean)) * iLocationFactor )))*(iEmission[0]))*CurStateMemoryblock3From[2];
                }
                iEmission[0] = iOne;
                if ((iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock2From = dp.StateMemoryblock2.read((iPos1-(-1))-(0));
                    hmmocMaxInPlace( CurStateMemoryblock2To[0], ((iTransition[1])*(iEmission[0]))*CurStateMemoryblock2From[0] );
                }
                iEmission[0] = ((iEquilibrium > aHalfQuality1[ (iPos0)-0 ]) ? iEquilibrium : aHalfQuality1[ (iPos0)-0 ])  /* Was: iEquilibrium */;
                if ((iPos0+1<=iLen1+0)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(-1))-(1), (iPos1-(0))-(0));
                    hmmocMaxInPlace( CurStateMemoryblock2To[0], (((iGapOpen * logspace::doubleexp( (((iPos1)-0 - iStartMean) * ((iPos1)-0 - iStartMean)) * iLocationFactor )))*(iEmission[0]))*CurStateMemoryblock3From[1] );
                }
                dp.StateMemoryblock2.written();
            }
            if ((iPos0+0<=0)&&(iPos1+0<=0)) {
                CurStateMemoryblock1To = dp.StateMemoryblock1.write();
                iEmission[0] = iOne;
                if (1) {
                    CurStateMemoryblock2From = dp.StateMemoryblock2.read((iPos1-(0))-(0));
                    CurStateMemoryblock1To[0] = ((iTransition[0])*(iEmission[0]))*CurStateMemoryblock2From[0];
                }
                dp.StateMemoryblock1.written();
            }
            iPrevSlowCoord = iPos1;
        }
    }
    iPrevSlowCoord = -1;
    {
        int iPos1=0;
        if (iPos1==iPos1) {} // avoid 'unused variable' warnings
        {
            int iPos0=0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            CurStateMemoryblock1From = dp.StateMemoryblock1.read();
            iTempProb[0] = CurStateMemoryblock1From[0];
        }
    }
    *ppOutTable = new AlignDPTable(dp);
    // make sure tables don't get deleted
    dp.isInCharge = false;
    return iTempProb[0];
};





//EDITED
Path& Viterbi_trace(AlignDPTable* pInTable,char* pQuality1,char* pSequence1,char* pSequence2,int iGapExtendPhred,int iGapOpenPhred,int iInsNucPrior,int iLen1,int iLen2,int iStartMean,int iStartSD, double ambiguity_bias) {
    logspace iTransition[16];
    const logspace *CurStateMemoryblock1To;
    const logspace *CurStateMemoryblock2To;
    const logspace *CurStateMemoryblock3To;
    const logspace *CurStateMemoryblock4To;
    const logspace *CurStateMemoryblock5To;
    // EDITED
    logspace iAmb = ambiguity_bias;
    int iPrevSlowCoord;
    SimplePath* pPath = new SimplePath();
    vector<int> emit;
    /* initialization of various probabilities */
    
    static const double db = 10.0 / log(0.1);
    static const logspace iOne = 1.0;
    logspace iEquilibrium = exp( iInsNucPrior / db );     // likelihood of any unaligned read base
    logspace aQuality1[ iLen1 ];
    logspace aHalfQuality1[ iLen1 ];
    for (int i=0; i<iLen1; i++) {
        if (pQuality1[i] < (char)33) {
            cout << "Internal error: Negative quality at position "<<i<<endl;
        }
        aQuality1[i] = exp( (pQuality1[i]-33) / db );
        aHalfQuality1[i] = exp( (pQuality1[i]-33) / (2.0 * db) );
        if (pSequence1[i] == 'N') {
            aQuality1[i] = iOne;
            aHalfQuality1[i] = iOne;
        }
    }
    logspace iGapExtend = exp( iGapExtendPhred / db );
    logspace iGapOpen = exp( iGapOpenPhred / db ) * (iOne - iGapExtend);
    double iLocationFactor = -1.0 / (2.0*iStartSD*iStartSD);
    int iSymbol[2];
    if (false && iSymbol[0] == iSymbol[0]) {}   // avoid 'unused variable' warnings
    logspace iEmission[1];
    /* temporary vector storage */
    logspace iTempVector[11];
    /* temporary int vector storage */
    int iTempIntVec[7];
    /* temporary storage for ordinary reals */
    register logspace iTempResult[1];
    iTransition[0] = iOne;
    
    iTransition[1] = iOne;
    
    iTransition[4] = iOne;
    
    iTransition[5] = iGapOpen;
    
    iTransition[6] = iGapOpen;
    
    iTransition[7] = iOne;
    
    iTransition[8] = iOne;
    
    iTransition[9] = iGapExtend;
    
    iTransition[10] = iGapOpen;
    
    iTransition[11] = iOne;
    
    iTransition[12] = iOne;
    
    iTransition[13] = iGapExtend;
    
    iTransition[14] = iOne;
    
    iTransition[15] = iOne;
    static const int stateTable[] = {1, 1, 4, 3, 4, 3, 2, 5, 4, 3, 2, 5, 4, 2, 5, 6};
    static const int stateFromTable[] = {0, 1, 1, 1, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 5, 5};
    static const int iPos0Table[] = {0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0};
    static const int iPos1Table[] = {0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0};
    AlignDPTable dp(*pInTable);
    // make sure tables don't get deleted
    dp.isInCharge = false;
    iPrevSlowCoord = -1;
    {
        int iPos1=0;
        if (iPos1==iPos1) {} // avoid 'unused variable' warnings
        {
            int iPos0=0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            iTempIntVec[0] = 0;
            while (iTempIntVec[0] != 6) {
                iTempIntVec[1] = 2;
                if ((iPos1+0<=iLen2+-1)) {
                    iSymbol[0] = pSequence2[iPos1+0];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                if ((iPos0+0<=iLen1+-1)) {
                    iSymbol[1] = pSequence1[iPos0+0];
                } 
                else { 
                    iSymbol[1] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock1To = dp.StateMemoryblock1.read();
                CurStateMemoryblock2To = dp.StateMemoryblock2.read((iPos1-(0))-(0));
                if ((iPos0+0<=0)&&(iPos1+1<=iLen2+0)) {
                    iEmission[0] = iOne;
                    CurStateMemoryblock2To = dp.StateMemoryblock2.read((iPos1-(-1))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 1:
                        iTempVector[iTempIntVec[1]] = iTransition[1]*iEmission[0]*CurStateMemoryblock2To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[1]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 1;
                        break;
                    }
                }
                if ((iPos0+0<=0)) {
                    iEmission[0] = iOne;
                    CurStateMemoryblock2To = dp.StateMemoryblock2.read((iPos1-(0))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 0:
                        iTempVector[iTempIntVec[1]] = iTransition[0]*iEmission[0]*CurStateMemoryblock2To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[0]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 0;
                        break;
                    }
                }
                CurStateMemoryblock3To = dp.StateMemoryblock3.read((iPos0-(0))-(1), (iPos1-(0))-(0));
                if ((iPos0+1<=iLen1+0)&&(iPos1+1<=iLen2+0)) {
                    if ((iSymbol[1] == iSymbol[0]) || (iSymbol[0] == 'N')) 
                    iTempResult[0] = iOne;
                    else 
                    iTempResult[0] = aQuality1[ (iPos0)-0 ];
                    iEmission[0] = iTempResult[0];
                    CurStateMemoryblock3To = dp.StateMemoryblock3.read((iPos0-(-1))-(1), (iPos1-(-1))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 3:
                        iTempVector[iTempIntVec[1]] = iTransition[8]*iEmission[0]*CurStateMemoryblock3To[2];
                        iTempVector[iTempIntVec[1]+4] = iTransition[8]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 8;
                        break;
                        case 1:
                        iTempVector[iTempIntVec[1]] = (logspace::doubleexp( (((iPos1)-1 - iStartMean) * ((iPos1)-1 - iStartMean)) * iLocationFactor ))*iEmission[0]*CurStateMemoryblock3To[2];
                        iTempVector[iTempIntVec[1]+4] = (logspace::doubleexp( (((iPos1)-1 - iStartMean) * ((iPos1)-1 - iStartMean)) * iLocationFactor ))*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 2;
                        break;
                        case 4:
                        iTempVector[iTempIntVec[1]] = iTransition[4]*iEmission[0]*CurStateMemoryblock3To[2];
                        iTempVector[iTempIntVec[1]+4] = iTransition[4]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 4;
                        break;
                        case 2:
                        iTempVector[iTempIntVec[1]] = iTransition[12]*iEmission[0]*CurStateMemoryblock3To[2];
                        iTempVector[iTempIntVec[1]+4] = iTransition[12]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 12;
                        break;
                    }
                }
                if ((iPos0+0>=1)&&(iPos1+1<=iLen2+0)) {
                    iEmission[0] = iOne;
                    CurStateMemoryblock3To = dp.StateMemoryblock3.read((iPos0-(0))-(1), (iPos1-(-1))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 4:
                        iTempVector[iTempIntVec[1]] = iTransition[6]*iEmission[0]*CurStateMemoryblock3To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[6]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 6;
                        break;
                        case 2:
                        iTempVector[iTempIntVec[1]] = iTransition[13]*iEmission[0]*CurStateMemoryblock3To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[13]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 13;
                        break;
                        case 3:
                        iTempVector[iTempIntVec[1]] = iTransition[10]*iEmission[0]*CurStateMemoryblock3To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[10]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 10;
                        break;
                    }
                }
                if ((iPos0+1<=iLen1+0)) {
                    iEmission[0] = ((iEquilibrium > aHalfQuality1[ (iPos0)-0 ]) ? iEquilibrium : aHalfQuality1[ (iPos0)-0 ])  /* Was: iEquilibrium */;
                    CurStateMemoryblock3To = dp.StateMemoryblock3.read((iPos0-(-1))-(1), (iPos1-(0))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 1:
                        iTempVector[iTempIntVec[1]] = (iGapOpen * logspace::doubleexp( (((iPos1)-0 - iStartMean) * ((iPos1)-0 - iStartMean)) * iLocationFactor ))*iEmission[0]*CurStateMemoryblock3To[1];
                        iTempVector[iTempIntVec[1]+4] = (iGapOpen * logspace::doubleexp( (((iPos1)-0 - iStartMean) * ((iPos1)-0 - iStartMean)) * iLocationFactor ))*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 3;
                        break;
                        case 4:
                        iTempVector[iTempIntVec[1]] = iTransition[5]*iEmission[0]*CurStateMemoryblock3To[1];
                        iTempVector[iTempIntVec[1]+4] = iTransition[5]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 5;
                        break;
                        case 3:
                        iTempVector[iTempIntVec[1]] = iTransition[9]*iEmission[0]*CurStateMemoryblock3To[1];
                        iTempVector[iTempIntVec[1]+4] = iTransition[9]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 9;
                        break;
                    }
                }
                CurStateMemoryblock4To = dp.StateMemoryblock4.read((iPos1-(0))-(0));
                if ((iPos0+0>=iLen1+0)&&(iPos1+1<=iLen2+0)) {
                    iEmission[0] = iOne;
                    CurStateMemoryblock4To = dp.StateMemoryblock4.read((iPos1-(-1))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 5:
                        iTempVector[iTempIntVec[1]] = iTransition[14]*iEmission[0]*CurStateMemoryblock4To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[14]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 14;
                        break;
                    }
                }
                if ((iPos0+0>=iLen1+0)) {
                    iEmission[0] = iOne;
                    CurStateMemoryblock4To = dp.StateMemoryblock4.read((iPos1-(0))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 4:
                        iTempVector[iTempIntVec[1]] = iTransition[7]*iEmission[0]*CurStateMemoryblock4To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[7]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 7;
                        break;
                        case 3:
                        iTempVector[iTempIntVec[1]] = iTransition[11]*iEmission[0]*CurStateMemoryblock4To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[11]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 11;
                        break;
                    }
                }
                CurStateMemoryblock5To = dp.StateMemoryblock5.read();
                if ((iPos0+0>=iLen1+0)&&(iPos1+0>=iLen2+0)) {
                    iEmission[0] = iOne;
                    CurStateMemoryblock5To = dp.StateMemoryblock5.read();
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 5:
                        iTempVector[iTempIntVec[1]] = iTransition[15]*iEmission[0]*CurStateMemoryblock5To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[15]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 15;
                        break;
                    }
                }
                iTempVector[0] = 0.0;
                for (int i=2; i<iTempIntVec[1]; i++) {
		    //EDITED
                    if (iTempVector[i]*iAmb>iTempVector[0]) {
                        iTempVector[0]=iTempVector[i];
                        iTempIntVec[0] = i;
                    }
                }
                emit.resize(2);
                emit[0] = iPos0Table[iTempIntVec[iTempIntVec[0]]];
                emit[1] = iPos1Table[iTempIntVec[iTempIntVec[0]]];
                pPath->addEdge(iTempIntVec[iTempIntVec[0]],iTempVector[iTempIntVec[0]+4],emit,stateFromTable[iTempIntVec[iTempIntVec[0]]],stateTable[iTempIntVec[iTempIntVec[0]]]);
                iPos0 += iPos0Table[iTempIntVec[iTempIntVec[0]]];
                iPos1 += iPos1Table[iTempIntVec[iTempIntVec[0]]];
                iTempIntVec[0] = stateTable[iTempIntVec[iTempIntVec[0]]];
            }
        }
    }
    return *pPath;
};



const string _AlignBandingstateId[] = {"start","pad1","delete","insert","match","pad2","end"};
const string _AlignBandingemissionId[] = {"emit1","emit12","emit2free","empty"};
const string _AlignBandingtransitionId[] = {"trSS","trPP1","trPM","trPI","trMM","trMI","trMD","trMP","trIM","trII","trIE","trIP","trDM","trDD","trPP2","trPE"};
const string _AlignBandingtransF[] = {"start","pad1","pad1","pad1","match","match","match","match","insert","insert","insert","insert","delete","delete","pad2","pad2"};
const string _AlignBandingtransT[] = {"pad1","pad1","match","insert","match","insert","delete","pad2","match","insert","delete","pad2","match","delete","pad2","end"};
const string _AlignBandingtransP[] = {"probSkip","probSkip","probMatch_Begin","probGapOpen_Begin","probMatch","probGapOpen","probGapOpen","probSkip","probGapEnd","probGapExtend","probGapOpen","probSkip","probGapEnd","probGapExtend","probSkip","probSkip"};
const string _AlignBandingtransE[] = {"empty","emit2free","emit12","emit1","emit12","emit1","emit2free","empty","emit12","emit1","emit2free","empty","emit12","emit2free","emit2free","empty"};
const string _AlignBandingoutputId[] = {"sequence1","sequence2"};
const string _AlignBandingempty = "";
const int _AlignBandingstateNum = 7;
const int _AlignBandingemitNum = 4;
const int _AlignBandingtransNum = 16;
const int _AlignBandingoutputNum = 2;




logspace ViterbiBanding_recurse(AlignBandingDPTable** ppOutTable,char* pQuality1,char* pSequence1,char* pSequence2,int iGapExtendPhred,int iGapOpenPhred,int iInsNucPrior,int iLen1,int iLen2,int iStartMean,int iStartSD,int iWidth) {
    logspace iTransition[16];
    logspace *CurStateMemoryblock4To;
    const logspace *CurStateMemoryblock4From;
    const logspace *CurStateMemoryblock5From;
    logspace *CurStateMemoryblock3withbandingTo;
    const logspace *CurStateMemoryblock3withbandingFrom;
    logspace *CurStateMemoryblock2To;
    const logspace *CurStateMemoryblock2From;
    logspace *CurStateMemoryblock1To;
    const logspace *CurStateMemoryblock1From;
    int iPrevSlowCoord;
    /* initialization of various probabilities */
    
    static const double db = 10.0 / log(0.1);
    static const logspace iOne = 1.0;
    logspace iEquilibrium = exp( iInsNucPrior / db );     // likelihood of any unaligned read base
    logspace aQuality1[ iLen1 ];
    logspace aHalfQuality1[ iLen1 ];
    for (int i=0; i<iLen1; i++) {
        if (pQuality1[i] < (char)33) {
            cout << "Internal error: Negative quality at position "<<i<<endl;
        }
        aQuality1[i] = exp( (pQuality1[i]-33) / db );
        aHalfQuality1[i] = exp( (pQuality1[i]-33) / (2.0 * db) );
        if (pSequence1[i] == 'N') {
            aQuality1[i] = iOne;
            aHalfQuality1[i] = iOne;
        }
    }
    logspace iGapExtend = exp( iGapExtendPhred / db );
    logspace iGapOpen = exp( iGapOpenPhred / db ) * (iOne - iGapExtend);
    double iLocationFactor = -1.0 / (2.0*iStartSD*iStartSD);
    DiagonalBanding bandingInstance( iLen1, iLen2, iWidth );
    int iSymbol[2];
    if (false && iSymbol[0] == iSymbol[0]) {}   // avoid 'unused variable' warnings
    logspace iEmission[1];
    /* temporary storage for ordinary reals */
    register logspace iTempResult[1];
    /* temporary storage for extended-exponent reals */
    register logspace iTempProb[1];
    AlignBandingDPTable dp(iLen1,iLen2);
    iTransition[0] = iOne;
    
    iTransition[1] = iOne;
    
    iTransition[4] = iOne;
    
    iTransition[5] = iGapOpen;
    
    iTransition[6] = iGapOpen;
    
    iTransition[7] = iOne;
    
    iTransition[8] = iOne;
    
    iTransition[9] = iGapExtend;
    
    iTransition[10] = iGapOpen;
    
    iTransition[11] = iOne;
    
    iTransition[12] = iOne;
    
    iTransition[13] = iGapExtend;
    
    iTransition[14] = iOne;
    
    iTransition[15] = iOne;
    dp.StateMemoryblock5.write()[0] = 1.0;
    dp.StateMemoryblock5.written();
    iPrevSlowCoord = -1;
    for (int iPos1=(iLen2+1)-1; iPos1>=0; --iPos1) {
        {
            int iPos0=iLen1+0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            if ((iPos1+0>=iLen2+0)) {
            }
            if (1) {
                if ((iPos1+0<=iLen2+-1)) {
                    iSymbol[0] = pSequence2[iPos1+0];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock4To = dp.StateMemoryblock4.write((iPos1-(0))-(0));
                iEmission[0] = iOne;
                if ((iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock4From = dp.StateMemoryblock4.read((iPos1-(-1))-(0));
                    CurStateMemoryblock4To[0] = ((iTransition[14])*(iEmission[0]))*CurStateMemoryblock4From[0];
                }
                iEmission[0] = iOne;
                if ((iPos1+0>=iLen2+0)) {
                    CurStateMemoryblock5From = dp.StateMemoryblock5.read();
                    hmmocMaxInPlace( CurStateMemoryblock4To[0], ((iTransition[15])*(iEmission[0]))*CurStateMemoryblock5From[0] );
                }
                dp.StateMemoryblock4.written();
            }
            iPrevSlowCoord = iPos1;
        }
    }
    iPrevSlowCoord = -1;
    {
        Banding<2>::Position& position = bandingInstance.backwardIterator();
        int iCheckSlowCoordTraversal = -1;  
        do {
            if (iCheckSlowCoordTraversal != -1 && iCheckSlowCoordTraversal < position[1]) {
                cout << "WARNING: Banding (backward): Slowest coordinate be nonincreasing.  Perhaps forgot to specify speed of output coordinates?" << endl;
            }
            iCheckSlowCoordTraversal = position[1];
            if ((position[0]+0>=1)&&(position[0]+0<=iLen1+0)&&(position[1]+0>=0)&&(position[1]+0<=iLen2+0)) {
                if (1) {
                    if ((position[0]+0<=iLen1+-1)) {
                        iSymbol[0] = pSequence1[position[0]+0];
                    } 
                    else { 
                        iSymbol[0] = 'A' /* dummy value */;
                        
                    }
                    if ((position[1]+0<=iLen2+-1)) {
                        iSymbol[1] = pSequence2[position[1]+0];
                    } 
                    else { 
                        iSymbol[1] = 'A' /* dummy value */;
                        
                    }
                    CurStateMemoryblock3withbandingTo = dp.StateMemoryblock3withbanding.write((position[0]-(0))-(1), (position[1]-(0))-(0));
                    if ((iSymbol[0] == iSymbol[1]) || (iSymbol[1] == 'N')) 
                    iTempResult[0] = iOne;
                    else 
                    iTempResult[0] = aQuality1[ (position[0])-0 ];
                    iEmission[0] = iTempResult[0];
                    if ((position[0]+1<=iLen1+0)&&(position[1]+1<=iLen2+0)) {
                        CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((position[0]-(-1))-(1), (position[1]-(-1))-(0));
                        CurStateMemoryblock3withbandingTo[2] = ((iTransition[4])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[2];
                        CurStateMemoryblock3withbandingTo[0] = ((iTransition[12])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[2];
                        CurStateMemoryblock3withbandingTo[1] = ((iTransition[8])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[2];
                    }
                    iEmission[0] = iOne;
                    if ((position[1]+1<=iLen2+0)) {
                        CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((position[0]-(0))-(1), (position[1]-(-1))-(0));
                        hmmocMaxInPlace( CurStateMemoryblock3withbandingTo[2], ((iTransition[6])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[0] );
                        hmmocMaxInPlace( CurStateMemoryblock3withbandingTo[0], ((iTransition[13])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[0] );
                        hmmocMaxInPlace( CurStateMemoryblock3withbandingTo[1], ((iTransition[10])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[0] );
                    }
                    iEmission[0] = ((iEquilibrium > aHalfQuality1[ (position[0])-0 ]) ? iEquilibrium : aHalfQuality1[ (position[0])-0 ])  /* Was: iEquilibrium */;
                    if ((position[0]+1<=iLen1+0)) {
                        CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((position[0]-(-1))-(1), (position[1]-(0))-(0));
                        hmmocMaxInPlace( CurStateMemoryblock3withbandingTo[2], ((iTransition[5])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[1] );
                        hmmocMaxInPlace( CurStateMemoryblock3withbandingTo[1], ((iTransition[9])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[1] );
                    }
                    iEmission[0] = iOne;
                    if ((position[0]+0>=iLen1+0)) {
                        CurStateMemoryblock4From = dp.StateMemoryblock4.read((position[1]-(0))-(0));
                        hmmocMaxInPlace( CurStateMemoryblock3withbandingTo[2], ((iTransition[7])*(iEmission[0]))*CurStateMemoryblock4From[0] );
                        hmmocMaxInPlace( CurStateMemoryblock3withbandingTo[1], ((iTransition[11])*(iEmission[0]))*CurStateMemoryblock4From[0] );
                    }
                    dp.StateMemoryblock3withbanding.written();
                }
                iPrevSlowCoord = position[1];
            } 
            else { 
                bandingInstance.warning();
                
            }
        } while (bandingInstance.hasNextBackward());
    }
    iPrevSlowCoord = -1;
    for (int iPos1=(iLen2+1)-1; iPos1>=0; --iPos1) {
        {
            int iPos0=0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            if (1) {
                if ((iPos0+0<=iLen1+-1)) {
                    iSymbol[0] = pSequence1[iPos0+0];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                if ((iPos1+0<=iLen2+-1)) {
                    iSymbol[1] = pSequence2[iPos1+0];
                } 
                else { 
                    iSymbol[1] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock2To = dp.StateMemoryblock2.write((iPos1-(0))-(0));
                if ((iSymbol[0] == iSymbol[1]) || (iSymbol[1] == 'N')) 
                iTempResult[0] = iOne;
                else 
                iTempResult[0] = aQuality1[ (iPos0)-0 ];
                iEmission[0] = iTempResult[0];
                if ((iPos0+1<=iLen1+0)&&(iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((iPos0-(-1))-(1), (iPos1-(-1))-(0));
                    CurStateMemoryblock2To[0] = (((logspace::doubleexp( (((iPos1)-1 - iStartMean) * ((iPos1)-1 - iStartMean)) * iLocationFactor )))*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[2];
                }
                iEmission[0] = iOne;
                if ((iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock2From = dp.StateMemoryblock2.read((iPos1-(-1))-(0));
                    hmmocMaxInPlace( CurStateMemoryblock2To[0], ((iTransition[1])*(iEmission[0]))*CurStateMemoryblock2From[0] );
                }
                iEmission[0] = ((iEquilibrium > aHalfQuality1[ (iPos0)-0 ]) ? iEquilibrium : aHalfQuality1[ (iPos0)-0 ])  /* Was: iEquilibrium */;
                if ((iPos0+1<=iLen1+0)) {
                    CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((iPos0-(-1))-(1), (iPos1-(0))-(0));
                    hmmocMaxInPlace( CurStateMemoryblock2To[0], (((iGapOpen * logspace::doubleexp( (((iPos1)-0 - iStartMean) * ((iPos1)-0 - iStartMean)) * iLocationFactor )))*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[1] );
                }
                dp.StateMemoryblock2.written();
            }
            if ((iPos1+0<=0)) {
                CurStateMemoryblock1To = dp.StateMemoryblock1.write();
                iEmission[0] = iOne;
                if (1) {
                    CurStateMemoryblock2From = dp.StateMemoryblock2.read((iPos1-(0))-(0));
                    CurStateMemoryblock1To[0] = ((iTransition[0])*(iEmission[0]))*CurStateMemoryblock2From[0];
                }
                dp.StateMemoryblock1.written();
            }
            iPrevSlowCoord = iPos1;
        }
    }
    iPrevSlowCoord = -1;
    {
        int iPos1=0;
        if (iPos1==iPos1) {} // avoid 'unused variable' warnings
        {
            int iPos0=0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            CurStateMemoryblock1From = dp.StateMemoryblock1.read();
            iTempProb[0] = CurStateMemoryblock1From[0];
        }
    }
    *ppOutTable = new AlignBandingDPTable(dp);
    // make sure tables don't get deleted
    dp.isInCharge = false;
    return iTempProb[0];
};





//EDITED
Path& ViterbiBanding_trace(AlignBandingDPTable* pInTable,char* pQuality1,char* pSequence1,char* pSequence2,int iGapExtendPhred,int iGapOpenPhred,int iInsNucPrior,int iLen1,int iLen2,int iStartMean,int iStartSD,double ambiguity_bias) {
    logspace iTransition[16];
    const logspace *CurStateMemoryblock1To;
    const logspace *CurStateMemoryblock2To;
    const logspace *CurStateMemoryblock3withbandingTo;
    const logspace *CurStateMemoryblock4To;
    const logspace *CurStateMemoryblock5To;
    // EDITED
    logspace iAmb = ambiguity_bias;
    int iPrevSlowCoord;
    SimplePath* pPath = new SimplePath();
    vector<int> emit;
    /* initialization of various probabilities */
    
    static const double db = 10.0 / log(0.1);
    static const logspace iOne = 1.0;
    logspace iEquilibrium = exp( iInsNucPrior / db );     // likelihood of any unaligned read base
    logspace aQuality1[ iLen1 ];
    logspace aHalfQuality1[ iLen1 ];
    for (int i=0; i<iLen1; i++) {
        if (pQuality1[i] < (char)33) {
            cout << "Internal error: Negative quality at position "<<i<<endl;
        }
        aQuality1[i] = exp( (pQuality1[i]-33) / db );
        aHalfQuality1[i] = exp( (pQuality1[i]-33) / (2.0 * db) );
        if (pSequence1[i] == 'N') {
            aQuality1[i] = iOne;
            aHalfQuality1[i] = iOne;
        }
    }
    logspace iGapExtend = exp( iGapExtendPhred / db );
    logspace iGapOpen = exp( iGapOpenPhred / db ) * (iOne - iGapExtend);
    double iLocationFactor = -1.0 / (2.0*iStartSD*iStartSD);
    int iSymbol[2];
    if (false && iSymbol[0] == iSymbol[0]) {}   // avoid 'unused variable' warnings
    logspace iEmission[1];
    /* temporary vector storage */
    logspace iTempVector[11];
    /* temporary int vector storage */
    int iTempIntVec[7];
    /* temporary storage for ordinary reals */
    register logspace iTempResult[1];
    iTransition[0] = iOne;
    
    iTransition[1] = iOne;
    
    iTransition[4] = iOne;
    
    iTransition[5] = iGapOpen;
    
    iTransition[6] = iGapOpen;
    
    iTransition[7] = iOne;
    
    iTransition[8] = iOne;
    
    iTransition[9] = iGapExtend;
    
    iTransition[10] = iGapOpen;
    
    iTransition[11] = iOne;
    
    iTransition[12] = iOne;
    
    iTransition[13] = iGapExtend;
    
    iTransition[14] = iOne;
    
    iTransition[15] = iOne;
    static const int stateTable[] = {1, 1, 4, 3, 4, 3, 2, 5, 4, 3, 2, 5, 4, 2, 5, 6};
    static const int stateFromTable[] = {0, 1, 1, 1, 4, 4, 4, 4, 3, 3, 3, 3, 2, 2, 5, 5};
    static const int iPos0Table[] = {0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0};
    static const int iPos1Table[] = {0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0};
    AlignBandingDPTable dp(*pInTable);
    // make sure tables don't get deleted
    dp.isInCharge = false;
    iPrevSlowCoord = -1;
    {
        int iPos1=0;
        if (iPos1==iPos1) {} // avoid 'unused variable' warnings
        {
            int iPos0=0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            iTempIntVec[0] = 0;
            while (iTempIntVec[0] != 6) {
                iTempIntVec[1] = 2;
                if ((iPos1+0<=iLen2+-1)) {
                    iSymbol[0] = pSequence2[iPos1+0];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                if ((iPos0+0<=iLen1+-1)) {
                    iSymbol[1] = pSequence1[iPos0+0];
                } 
                else { 
                    iSymbol[1] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock1To = dp.StateMemoryblock1.read();
                CurStateMemoryblock2To = dp.StateMemoryblock2.read((iPos1-(0))-(0));
                if ((iPos0+0<=0)&&(iPos1+1<=iLen2+0)) {
                    iEmission[0] = iOne;
                    CurStateMemoryblock2To = dp.StateMemoryblock2.read((iPos1-(-1))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 1:
                        iTempVector[iTempIntVec[1]] = iTransition[1]*iEmission[0]*CurStateMemoryblock2To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[1]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 1;
                        break;
                    }
                }
                if ((iPos0+0<=0)) {
                    iEmission[0] = iOne;
                    CurStateMemoryblock2To = dp.StateMemoryblock2.read((iPos1-(0))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 0:
                        iTempVector[iTempIntVec[1]] = iTransition[0]*iEmission[0]*CurStateMemoryblock2To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[0]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 0;
                        break;
                    }
                }
                CurStateMemoryblock3withbandingTo = dp.StateMemoryblock3withbanding.read((iPos0-(0))-(1), (iPos1-(0))-(0));
                if ((iPos0+1<=iLen1+0)&&(iPos1+1<=iLen2+0)) {
                    if ((iSymbol[1] == iSymbol[0]) || (iSymbol[0] == 'N')) 
                    iTempResult[0] = iOne;
                    else 
                    iTempResult[0] = aQuality1[ (iPos0)-0 ];
                    iEmission[0] = iTempResult[0];
                    CurStateMemoryblock3withbandingTo = dp.StateMemoryblock3withbanding.read((iPos0-(-1))-(1), (iPos1-(-1))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 1:
                        iTempVector[iTempIntVec[1]] = (logspace::doubleexp( (((iPos1)-1 - iStartMean) * ((iPos1)-1 - iStartMean)) * iLocationFactor ))*iEmission[0]*CurStateMemoryblock3withbandingTo[2];
                        iTempVector[iTempIntVec[1]+4] = (logspace::doubleexp( (((iPos1)-1 - iStartMean) * ((iPos1)-1 - iStartMean)) * iLocationFactor ))*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 2;
                        break;
                        case 2:
                        iTempVector[iTempIntVec[1]] = iTransition[12]*iEmission[0]*CurStateMemoryblock3withbandingTo[2];
                        iTempVector[iTempIntVec[1]+4] = iTransition[12]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 12;
                        break;
                        case 4:
                        iTempVector[iTempIntVec[1]] = iTransition[4]*iEmission[0]*CurStateMemoryblock3withbandingTo[2];
                        iTempVector[iTempIntVec[1]+4] = iTransition[4]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 4;
                        break;
                        case 3:
                        iTempVector[iTempIntVec[1]] = iTransition[8]*iEmission[0]*CurStateMemoryblock3withbandingTo[2];
                        iTempVector[iTempIntVec[1]+4] = iTransition[8]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 8;
                        break;
                    }
                }
                if ((iPos0+0>=1)&&(iPos1+1<=iLen2+0)) {
                    iEmission[0] = iOne;
                    CurStateMemoryblock3withbandingTo = dp.StateMemoryblock3withbanding.read((iPos0-(0))-(1), (iPos1-(-1))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 2:
                        iTempVector[iTempIntVec[1]] = iTransition[13]*iEmission[0]*CurStateMemoryblock3withbandingTo[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[13]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 13;
                        break;
                        case 4:
                        iTempVector[iTempIntVec[1]] = iTransition[6]*iEmission[0]*CurStateMemoryblock3withbandingTo[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[6]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 6;
                        break;
                        case 3:
                        iTempVector[iTempIntVec[1]] = iTransition[10]*iEmission[0]*CurStateMemoryblock3withbandingTo[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[10]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 10;
                        break;
                    }
                }
                if ((iPos0+1<=iLen1+0)) {
                    iEmission[0] = ((iEquilibrium > aHalfQuality1[ (iPos0)-0 ]) ? iEquilibrium : aHalfQuality1[ (iPos0)-0 ])  /* Was: iEquilibrium */;
                    CurStateMemoryblock3withbandingTo = dp.StateMemoryblock3withbanding.read((iPos0-(-1))-(1), (iPos1-(0))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 1:
                        iTempVector[iTempIntVec[1]] = (iGapOpen * logspace::doubleexp( (((iPos1)-0 - iStartMean) * ((iPos1)-0 - iStartMean)) * iLocationFactor ))*iEmission[0]*CurStateMemoryblock3withbandingTo[1];
                        iTempVector[iTempIntVec[1]+4] = (iGapOpen * logspace::doubleexp( (((iPos1)-0 - iStartMean) * ((iPos1)-0 - iStartMean)) * iLocationFactor ))*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 3;
                        break;
                        case 4:
                        iTempVector[iTempIntVec[1]] = iTransition[5]*iEmission[0]*CurStateMemoryblock3withbandingTo[1];
                        iTempVector[iTempIntVec[1]+4] = iTransition[5]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 5;
                        break;
                        case 3:
                        iTempVector[iTempIntVec[1]] = iTransition[9]*iEmission[0]*CurStateMemoryblock3withbandingTo[1];
                        iTempVector[iTempIntVec[1]+4] = iTransition[9]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 9;
                        break;
                    }
                }
                CurStateMemoryblock4To = dp.StateMemoryblock4.read((iPos1-(0))-(0));
                if ((iPos0+0>=iLen1+0)&&(iPos1+1<=iLen2+0)) {
                    iEmission[0] = iOne;
                    CurStateMemoryblock4To = dp.StateMemoryblock4.read((iPos1-(-1))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 5:
                        iTempVector[iTempIntVec[1]] = iTransition[14]*iEmission[0]*CurStateMemoryblock4To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[14]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 14;
                        break;
                    }
                }
                if ((iPos0+0>=iLen1+0)) {
                    iEmission[0] = iOne;
                    CurStateMemoryblock4To = dp.StateMemoryblock4.read((iPos1-(0))-(0));
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 4:
                        iTempVector[iTempIntVec[1]] = iTransition[7]*iEmission[0]*CurStateMemoryblock4To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[7]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 7;
                        break;
                        case 3:
                        iTempVector[iTempIntVec[1]] = iTransition[11]*iEmission[0]*CurStateMemoryblock4To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[11]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 11;
                        break;
                    }
                }
                CurStateMemoryblock5To = dp.StateMemoryblock5.read();
                if ((iPos0+0>=iLen1+0)&&(iPos1+0>=iLen2+0)) {
                    iEmission[0] = iOne;
                    CurStateMemoryblock5To = dp.StateMemoryblock5.read();
                    switch (iTempIntVec[0]) {
                        default:
                        break;
                        case 5:
                        iTempVector[iTempIntVec[1]] = iTransition[15]*iEmission[0]*CurStateMemoryblock5To[0];
                        iTempVector[iTempIntVec[1]+4] = iTransition[15]*iEmission[0];
                        iTempIntVec[iTempIntVec[1]++] = 15;
                        break;
                    }
                }
                iTempVector[0] = 0.0;
                for (int i=2; i<iTempIntVec[1]; i++) {
		    //EDITED
                    if (iTempVector[i]*iAmb>iTempVector[0]) {
                        iTempVector[0]=iTempVector[i];
                        iTempIntVec[0] = i;
                    }
                }
                emit.resize(2);
                emit[0] = iPos0Table[iTempIntVec[iTempIntVec[0]]];
                emit[1] = iPos1Table[iTempIntVec[iTempIntVec[0]]];
                pPath->addEdge(iTempIntVec[iTempIntVec[0]],iTempVector[iTempIntVec[0]+4],emit,stateFromTable[iTempIntVec[iTempIntVec[0]]],stateTable[iTempIntVec[iTempIntVec[0]]]);
                iPos0 += iPos0Table[iTempIntVec[iTempIntVec[0]]];
                iPos1 += iPos1Table[iTempIntVec[iTempIntVec[0]]];
                iTempIntVec[0] = stateTable[iTempIntVec[iTempIntVec[0]]];
            }
        }
    }
    return *pPath;
};





logspace Forward(AlignDPTable** ppOutTable,char* pQuality1,char* pSequence1,char* pSequence2,int iGapExtendPhred,int iGapOpenPhred,int iInsNucPrior,int iLen1,int iLen2,int iStartMean,int iStartSD) {
    logspace iTransition[16];
    logspace *CurStateMemoryblock2To;
    const logspace *CurStateMemoryblock2From;
    const logspace *CurStateMemoryblock1From;
    logspace *CurStateMemoryblock3To;
    const logspace *CurStateMemoryblock3From;
    logspace *CurStateMemoryblock4To;
    const logspace *CurStateMemoryblock4From;
    logspace *CurStateMemoryblock5To;
    const logspace *CurStateMemoryblock5From;
    int iPrevSlowCoord;
    /* initialization of various probabilities */
    
    static const double db = 10.0 / log(0.1);
    static const logspace iOne = 1.0;
    logspace iEquilibrium = exp( iInsNucPrior / db );     // likelihood of any unaligned read base
    logspace aQuality1[ iLen1 ];
    logspace aHalfQuality1[ iLen1 ];
    for (int i=0; i<iLen1; i++) {
        if (pQuality1[i] < (char)33) {
            cout << "Internal error: Negative quality at position "<<i<<endl;
        }
        aQuality1[i] = exp( (pQuality1[i]-33) / db );
        aHalfQuality1[i] = exp( (pQuality1[i]-33) / (2.0 * db) );
        if (pSequence1[i] == 'N') {
            aQuality1[i] = iOne;
            aHalfQuality1[i] = iOne;
        }
    }
    logspace iGapExtend = exp( iGapExtendPhred / db );
    logspace iGapOpen = exp( iGapOpenPhred / db ) * (iOne - iGapExtend);
    double iLocationFactor = -1.0 / (2.0*iStartSD*iStartSD);
    int iSymbol[2];
    if (false && iSymbol[0] == iSymbol[0]) {}   // avoid 'unused variable' warnings
    logspace iEmission[1];
    /* temporary storage for ordinary reals */
    register logspace iTempResult[1];
    /* temporary storage for extended-exponent reals */
    register logspace iTempProb[1];
    AlignDPTable dp(iLen1,iLen2);
    iTransition[0] = iOne;
    
    iTransition[1] = iOne;
    
    iTransition[4] = iOne;
    
    iTransition[5] = iGapOpen;
    
    iTransition[6] = iGapOpen;
    
    iTransition[7] = iOne;
    
    iTransition[8] = iOne;
    
    iTransition[9] = iGapExtend;
    
    iTransition[10] = iGapOpen;
    
    iTransition[11] = iOne;
    
    iTransition[12] = iOne;
    
    iTransition[13] = iGapExtend;
    
    iTransition[14] = iOne;
    
    iTransition[15] = iOne;
    dp.StateMemoryblock1.write()[0] = 1.0;
    dp.StateMemoryblock1.written();
    iPrevSlowCoord = -1;
    for (int iPos1=0; iPos1<iLen2+1; ++iPos1) {
        for (int iPos0=0; iPos0<iLen1+1; ++iPos0) {
            if ((iPos0+0<=0)&&(iPos1+0<=0)) {
            }
            if ((iPos0+0<=0)) {
                if ((iPos1+-1>=0)) {
                    iSymbol[0] = pSequence2[iPos1+-1];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock2To = dp.StateMemoryblock2.write((iPos1-(0))-(0));
                iEmission[0] = iOne;
                if ((iPos1+-1>=0)) {
                    CurStateMemoryblock2From = dp.StateMemoryblock2.read((iPos1-(1))-(0));
                    CurStateMemoryblock2To[0] = ((iTransition[1])*(iEmission[0]))*CurStateMemoryblock2From[0];
                }
                iEmission[0] = iOne;
                if ((iPos1+0<=0)) {
                    CurStateMemoryblock1From = dp.StateMemoryblock1.read();
                    CurStateMemoryblock2To[0] += ((iTransition[0])*(iEmission[0]))*CurStateMemoryblock1From[0];
                }
                dp.StateMemoryblock2.written();
            }
            if ((iPos0+0>=1)) {
                if ((iPos0+-1>=0)) {
                    iSymbol[0] = pSequence1[iPos0+-1];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                if ((iPos1+-1>=0)) {
                    iSymbol[1] = pSequence2[iPos1+-1];
                } 
                else { 
                    iSymbol[1] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock3To = dp.StateMemoryblock3.write((iPos0-(0))-(1), (iPos1-(0))-(0));
                if ((iSymbol[0] == iSymbol[1]) || (iSymbol[1] == 'N')) 
                iTempResult[0] = iOne;
                else 
                iTempResult[0] = aQuality1[ (iPos0)-1 ];
                iEmission[0] = iTempResult[0];
                if ((iPos0+-1<=0)&&(iPos1+-1>=0)) {
                    CurStateMemoryblock2From = dp.StateMemoryblock2.read((iPos1-(1))-(0));
                    CurStateMemoryblock3To[2] = (((logspace::doubleexp( (((iPos1)-1 - iStartMean) * ((iPos1)-1 - iStartMean)) * iLocationFactor )))*(iEmission[0]))*CurStateMemoryblock2From[0];
                }
                if ((iPos0+-1>=1)&&(iPos1+-1>=0)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(1))-(1), (iPos1-(1))-(0));
                    CurStateMemoryblock3To[2] += ((iTransition[8])*(iEmission[0]))*CurStateMemoryblock3From[1];
                    CurStateMemoryblock3To[2] += ((iTransition[12])*(iEmission[0]))*CurStateMemoryblock3From[0];
                    CurStateMemoryblock3To[2] += ((iTransition[4])*(iEmission[0]))*CurStateMemoryblock3From[2];
                }
                iEmission[0] = iOne;
                if ((iPos1+-1>=0)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(0))-(1), (iPos1-(1))-(0));
                    CurStateMemoryblock3To[0] = ((iTransition[10])*(iEmission[0]))*CurStateMemoryblock3From[1];
                    CurStateMemoryblock3To[0] += ((iTransition[13])*(iEmission[0]))*CurStateMemoryblock3From[0];
                    CurStateMemoryblock3To[0] += ((iTransition[6])*(iEmission[0]))*CurStateMemoryblock3From[2];
                }
                iEmission[0] = ((iEquilibrium > aHalfQuality1[ (iPos0)-1 ]) ? iEquilibrium : aHalfQuality1[ (iPos0)-1 ])  /* Was: iEquilibrium */;
                if ((iPos0+-1<=0)) {
                    CurStateMemoryblock2From = dp.StateMemoryblock2.read((iPos1-(0))-(0));
                    CurStateMemoryblock3To[1] = (((iGapOpen * logspace::doubleexp( (((iPos1)-0 - iStartMean) * ((iPos1)-0 - iStartMean)) * iLocationFactor )))*(iEmission[0]))*CurStateMemoryblock2From[0];
                }
                if ((iPos0+-1>=1)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(1))-(1), (iPos1-(0))-(0));
                    CurStateMemoryblock3To[1] += ((iTransition[9])*(iEmission[0]))*CurStateMemoryblock3From[1];
                    CurStateMemoryblock3To[1] += ((iTransition[5])*(iEmission[0]))*CurStateMemoryblock3From[2];
                }
                dp.StateMemoryblock3.written();
            }
            if ((iPos0+0>=iLen1+0)) {
                if ((iPos1+-1>=0)) {
                    iSymbol[0] = pSequence2[iPos1+-1];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock4To = dp.StateMemoryblock4.write((iPos1-(0))-(0));
                iEmission[0] = iOne;
                if ((iPos1+-1>=0)) {
                    CurStateMemoryblock4From = dp.StateMemoryblock4.read((iPos1-(1))-(0));
                    CurStateMemoryblock4To[0] = ((iTransition[14])*(iEmission[0]))*CurStateMemoryblock4From[0];
                }
                iEmission[0] = iOne;
                if ((iPos0+0>=1)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(0))-(1), (iPos1-(0))-(0));
                    CurStateMemoryblock4To[0] += ((iTransition[11])*(iEmission[0]))*CurStateMemoryblock3From[1];
                    CurStateMemoryblock4To[0] += ((iTransition[7])*(iEmission[0]))*CurStateMemoryblock3From[2];
                }
                dp.StateMemoryblock4.written();
            }
            if ((iPos0+0>=iLen1+0)&&(iPos1+0>=iLen2+0)) {
                CurStateMemoryblock5To = dp.StateMemoryblock5.write();
                iEmission[0] = iOne;
                if (1) {
                    CurStateMemoryblock4From = dp.StateMemoryblock4.read((iPos1-(0))-(0));
                    CurStateMemoryblock5To[0] = ((iTransition[15])*(iEmission[0]))*CurStateMemoryblock4From[0];
                }
                dp.StateMemoryblock5.written();
            }
            iPrevSlowCoord = iPos1;
        }
    }
    iPrevSlowCoord = -1;
    {
        int iPos1=iLen2+0;
        if (iPos1==iPos1) {} // avoid 'unused variable' warnings
        {
            int iPos0=iLen1+0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            CurStateMemoryblock5From = dp.StateMemoryblock5.read();
            iTempProb[0] = CurStateMemoryblock5From[0];
        }
    }
    *ppOutTable = new AlignDPTable(dp);
    // make sure tables don't get deleted
    dp.isInCharge = false;
    return iTempProb[0];
};





logspace Backward(AlignDPTable** ppOutTable,char* pQuality1,char* pSequence1,char* pSequence2,int iGapExtendPhred,int iGapOpenPhred,int iInsNucPrior,int iLen1,int iLen2,int iStartMean,int iStartSD) {
    logspace iTransition[16];
    logspace *CurStateMemoryblock4To;
    const logspace *CurStateMemoryblock4From;
    const logspace *CurStateMemoryblock5From;
    logspace *CurStateMemoryblock3To;
    const logspace *CurStateMemoryblock3From;
    logspace *CurStateMemoryblock2To;
    const logspace *CurStateMemoryblock2From;
    logspace *CurStateMemoryblock1To;
    const logspace *CurStateMemoryblock1From;
    int iPrevSlowCoord;
    /* initialization of various probabilities */
    
    static const double db = 10.0 / log(0.1);
    static const logspace iOne = 1.0;
    logspace iEquilibrium = exp( iInsNucPrior / db );     // likelihood of any unaligned read base
    logspace aQuality1[ iLen1 ];
    logspace aHalfQuality1[ iLen1 ];
    for (int i=0; i<iLen1; i++) {
        if (pQuality1[i] < (char)33) {
            cout << "Internal error: Negative quality at position "<<i<<endl;
        }
        aQuality1[i] = exp( (pQuality1[i]-33) / db );
        aHalfQuality1[i] = exp( (pQuality1[i]-33) / (2.0 * db) );
        if (pSequence1[i] == 'N') {
            aQuality1[i] = iOne;
            aHalfQuality1[i] = iOne;
        }
    }
    logspace iGapExtend = exp( iGapExtendPhred / db );
    logspace iGapOpen = exp( iGapOpenPhred / db ) * (iOne - iGapExtend);
    double iLocationFactor = -1.0 / (2.0*iStartSD*iStartSD);
    int iSymbol[2];
    if (false && iSymbol[0] == iSymbol[0]) {}   // avoid 'unused variable' warnings
    logspace iEmission[1];
    /* temporary storage for ordinary reals */
    register logspace iTempResult[1];
    /* temporary storage for extended-exponent reals */
    register logspace iTempProb[1];
    AlignDPTable dp(iLen1,iLen2);
    iTransition[0] = iOne;
    
    iTransition[1] = iOne;
    
    iTransition[4] = iOne;
    
    iTransition[5] = iGapOpen;
    
    iTransition[6] = iGapOpen;
    
    iTransition[7] = iOne;
    
    iTransition[8] = iOne;
    
    iTransition[9] = iGapExtend;
    
    iTransition[10] = iGapOpen;
    
    iTransition[11] = iOne;
    
    iTransition[12] = iOne;
    
    iTransition[13] = iGapExtend;
    
    iTransition[14] = iOne;
    
    iTransition[15] = iOne;
    dp.StateMemoryblock5.write()[0] = 1.0;
    dp.StateMemoryblock5.written();
    iPrevSlowCoord = -1;
    for (int iPos1=(iLen2+1)-1; iPos1>=0; --iPos1) {
        for (int iPos0=(iLen1+1)-1; iPos0>=0; --iPos0) {
            if ((iPos0+0>=iLen1+0)&&(iPos1+0>=iLen2+0)) {
            }
            if ((iPos0+0>=iLen1+0)) {
                if ((iPos1+0<=iLen2+-1)) {
                    iSymbol[0] = pSequence2[iPos1+0];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock4To = dp.StateMemoryblock4.write((iPos1-(0))-(0));
                iEmission[0] = iOne;
                if ((iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock4From = dp.StateMemoryblock4.read((iPos1-(-1))-(0));
                    CurStateMemoryblock4To[0] = ((iTransition[14])*(iEmission[0]))*CurStateMemoryblock4From[0];
                }
                iEmission[0] = iOne;
                if ((iPos1+0>=iLen2+0)) {
                    CurStateMemoryblock5From = dp.StateMemoryblock5.read();
                    CurStateMemoryblock4To[0] += ((iTransition[15])*(iEmission[0]))*CurStateMemoryblock5From[0];
                }
                dp.StateMemoryblock4.written();
            }
            if ((iPos0+0>=1)) {
                if ((iPos0+0<=iLen1+-1)) {
                    iSymbol[0] = pSequence1[iPos0+0];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                if ((iPos1+0<=iLen2+-1)) {
                    iSymbol[1] = pSequence2[iPos1+0];
                } 
                else { 
                    iSymbol[1] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock3To = dp.StateMemoryblock3.write((iPos0-(0))-(1), (iPos1-(0))-(0));
                if ((iSymbol[0] == iSymbol[1]) || (iSymbol[1] == 'N')) 
                iTempResult[0] = iOne;
                else 
                iTempResult[0] = aQuality1[ (iPos0)-0 ];
                iEmission[0] = iTempResult[0];
                if ((iPos0+1<=iLen1+0)&&(iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(-1))-(1), (iPos1-(-1))-(0));
                    CurStateMemoryblock3To[1] = ((iTransition[8])*(iEmission[0]))*CurStateMemoryblock3From[2];
                    CurStateMemoryblock3To[0] = ((iTransition[12])*(iEmission[0]))*CurStateMemoryblock3From[2];
                    CurStateMemoryblock3To[2] = ((iTransition[4])*(iEmission[0]))*CurStateMemoryblock3From[2];
                }
                iEmission[0] = iOne;
                if ((iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(0))-(1), (iPos1-(-1))-(0));
                    CurStateMemoryblock3To[1] += ((iTransition[10])*(iEmission[0]))*CurStateMemoryblock3From[0];
                    CurStateMemoryblock3To[0] += ((iTransition[13])*(iEmission[0]))*CurStateMemoryblock3From[0];
                    CurStateMemoryblock3To[2] += ((iTransition[6])*(iEmission[0]))*CurStateMemoryblock3From[0];
                }
                iEmission[0] = ((iEquilibrium > aHalfQuality1[ (iPos0)-0 ]) ? iEquilibrium : aHalfQuality1[ (iPos0)-0 ])  /* Was: iEquilibrium */;
                if ((iPos0+1<=iLen1+0)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(-1))-(1), (iPos1-(0))-(0));
                    CurStateMemoryblock3To[1] += ((iTransition[9])*(iEmission[0]))*CurStateMemoryblock3From[1];
                    CurStateMemoryblock3To[2] += ((iTransition[5])*(iEmission[0]))*CurStateMemoryblock3From[1];
                }
                iEmission[0] = iOne;
                if ((iPos0+0>=iLen1+0)) {
                    CurStateMemoryblock4From = dp.StateMemoryblock4.read((iPos1-(0))-(0));
                    CurStateMemoryblock3To[1] += ((iTransition[11])*(iEmission[0]))*CurStateMemoryblock4From[0];
                    CurStateMemoryblock3To[2] += ((iTransition[7])*(iEmission[0]))*CurStateMemoryblock4From[0];
                }
                dp.StateMemoryblock3.written();
            }
            if ((iPos0+0<=0)) {
                if ((iPos0+0<=iLen1+-1)) {
                    iSymbol[0] = pSequence1[iPos0+0];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                if ((iPos1+0<=iLen2+-1)) {
                    iSymbol[1] = pSequence2[iPos1+0];
                } 
                else { 
                    iSymbol[1] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock2To = dp.StateMemoryblock2.write((iPos1-(0))-(0));
                if ((iSymbol[0] == iSymbol[1]) || (iSymbol[1] == 'N')) 
                iTempResult[0] = iOne;
                else 
                iTempResult[0] = aQuality1[ (iPos0)-0 ];
                iEmission[0] = iTempResult[0];
                if ((iPos0+1<=iLen1+0)&&(iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(-1))-(1), (iPos1-(-1))-(0));
                    CurStateMemoryblock2To[0] = (((logspace::doubleexp( (((iPos1)-1 - iStartMean) * ((iPos1)-1 - iStartMean)) * iLocationFactor )))*(iEmission[0]))*CurStateMemoryblock3From[2];
                }
                iEmission[0] = iOne;
                if ((iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock2From = dp.StateMemoryblock2.read((iPos1-(-1))-(0));
                    CurStateMemoryblock2To[0] += ((iTransition[1])*(iEmission[0]))*CurStateMemoryblock2From[0];
                }
                iEmission[0] = ((iEquilibrium > aHalfQuality1[ (iPos0)-0 ]) ? iEquilibrium : aHalfQuality1[ (iPos0)-0 ])  /* Was: iEquilibrium */;
                if ((iPos0+1<=iLen1+0)) {
                    CurStateMemoryblock3From = dp.StateMemoryblock3.read((iPos0-(-1))-(1), (iPos1-(0))-(0));
                    CurStateMemoryblock2To[0] += (((iGapOpen * logspace::doubleexp( (((iPos1)-0 - iStartMean) * ((iPos1)-0 - iStartMean)) * iLocationFactor )))*(iEmission[0]))*CurStateMemoryblock3From[1];
                }
                dp.StateMemoryblock2.written();
            }
            if ((iPos0+0<=0)&&(iPos1+0<=0)) {
                CurStateMemoryblock1To = dp.StateMemoryblock1.write();
                iEmission[0] = iOne;
                if (1) {
                    CurStateMemoryblock2From = dp.StateMemoryblock2.read((iPos1-(0))-(0));
                    CurStateMemoryblock1To[0] = ((iTransition[0])*(iEmission[0]))*CurStateMemoryblock2From[0];
                }
                dp.StateMemoryblock1.written();
            }
            iPrevSlowCoord = iPos1;
        }
    }
    iPrevSlowCoord = -1;
    {
        int iPos1=0;
        if (iPos1==iPos1) {} // avoid 'unused variable' warnings
        {
            int iPos0=0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            CurStateMemoryblock1From = dp.StateMemoryblock1.read();
            iTempProb[0] = CurStateMemoryblock1From[0];
        }
    }
    *ppOutTable = new AlignDPTable(dp);
    // make sure tables don't get deleted
    dp.isInCharge = false;
    return iTempProb[0];
};





logspace ForwardBanding(AlignBandingDPTable** ppOutTable,char* pQuality1,char* pSequence1,char* pSequence2,int iGapExtendPhred,int iGapOpenPhred,int iInsNucPrior,int iLen1,int iLen2,int iStartMean,int iStartSD,int iWidth) {
    logspace iTransition[16];
    logspace *CurStateMemoryblock2To;
    const logspace *CurStateMemoryblock2From;
    const logspace *CurStateMemoryblock1From;
    logspace *CurStateMemoryblock3withbandingTo;
    const logspace *CurStateMemoryblock3withbandingFrom;
    logspace *CurStateMemoryblock4To;
    const logspace *CurStateMemoryblock4From;
    logspace *CurStateMemoryblock5To;
    const logspace *CurStateMemoryblock5From;
    int iPrevSlowCoord;
    /* initialization of various probabilities */
    
    static const double db = 10.0 / log(0.1);
    static const logspace iOne = 1.0;
    logspace iEquilibrium = exp( iInsNucPrior / db );     // likelihood of any unaligned read base
    logspace aQuality1[ iLen1 ];
    logspace aHalfQuality1[ iLen1 ];
    for (int i=0; i<iLen1; i++) {
        if (pQuality1[i] < (char)33) {
            cout << "Internal error: Negative quality at position "<<i<<endl;
        }
        aQuality1[i] = exp( (pQuality1[i]-33) / db );
        aHalfQuality1[i] = exp( (pQuality1[i]-33) / (2.0 * db) );
        if (pSequence1[i] == 'N') {
            aQuality1[i] = iOne;
            aHalfQuality1[i] = iOne;
        }
    }
    logspace iGapExtend = exp( iGapExtendPhred / db );
    logspace iGapOpen = exp( iGapOpenPhred / db ) * (iOne - iGapExtend);
    double iLocationFactor = -1.0 / (2.0*iStartSD*iStartSD);
    DiagonalBanding bandingInstance( iLen1, iLen2, iWidth );
    int iSymbol[2];
    if (false && iSymbol[0] == iSymbol[0]) {}   // avoid 'unused variable' warnings
    logspace iEmission[1];
    /* temporary storage for ordinary reals */
    register logspace iTempResult[1];
    /* temporary storage for extended-exponent reals */
    register logspace iTempProb[1];
    AlignBandingDPTable dp(iLen1,iLen2);
    iTransition[0] = iOne;
    
    iTransition[1] = iOne;
    
    iTransition[4] = iOne;
    
    iTransition[5] = iGapOpen;
    
    iTransition[6] = iGapOpen;
    
    iTransition[7] = iOne;
    
    iTransition[8] = iOne;
    
    iTransition[9] = iGapExtend;
    
    iTransition[10] = iGapOpen;
    
    iTransition[11] = iOne;
    
    iTransition[12] = iOne;
    
    iTransition[13] = iGapExtend;
    
    iTransition[14] = iOne;
    
    iTransition[15] = iOne;
    dp.StateMemoryblock1.write()[0] = 1.0;
    dp.StateMemoryblock1.written();
    iPrevSlowCoord = -1;
    for (int iPos1=0; iPos1<iLen2+1; ++iPos1) {
        {
            int iPos0=0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            if ((iPos1+0<=0)) {
            }
            if (1) {
                if ((iPos1+-1>=0)) {
                    iSymbol[0] = pSequence2[iPos1+-1];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock2To = dp.StateMemoryblock2.write((iPos1-(0))-(0));
                iEmission[0] = iOne;
                if ((iPos1+-1>=0)) {
                    CurStateMemoryblock2From = dp.StateMemoryblock2.read((iPos1-(1))-(0));
                    CurStateMemoryblock2To[0] = ((iTransition[1])*(iEmission[0]))*CurStateMemoryblock2From[0];
                }
                iEmission[0] = iOne;
                if ((iPos1+0<=0)) {
                    CurStateMemoryblock1From = dp.StateMemoryblock1.read();
                    CurStateMemoryblock2To[0] += ((iTransition[0])*(iEmission[0]))*CurStateMemoryblock1From[0];
                }
                dp.StateMemoryblock2.written();
            }
            iPrevSlowCoord = iPos1;
        }
    }
    iPrevSlowCoord = -1;
    {
        Banding<2>::Position& position = bandingInstance.forwardIterator();
        bool bLastSlowCoordInited = false;
        int iLastSlowCoord = -1;  
        do {
            if (bLastSlowCoordInited) {
                if (iLastSlowCoord > position[1]) {
                    cout << "WARNING: Banding (forward): Slowest coordinate should be nondecreasing.  Perhaps forgot to specify speed of output coordinates?" << endl;
                }
                } else {
                bLastSlowCoordInited = true;
            }
            iLastSlowCoord = position[1];
            if ((position[0]+0>=1)&&(position[0]+0<=iLen1+0)&&(position[1]+0>=0)&&(position[1]+0<=iLen2+0)) {
                if (1) {
                    if (1) {
                        iSymbol[0] = pSequence1[position[0]+-1];
                    } 
                    else { 
                        iSymbol[0] = 'A' /* dummy value */;
                        
                    }
                    if ((position[1]+-1>=0)) {
                        iSymbol[1] = pSequence2[position[1]+-1];
                    } 
                    else { 
                        iSymbol[1] = 'A' /* dummy value */;
                        
                    }
                    CurStateMemoryblock3withbandingTo = dp.StateMemoryblock3withbanding.write((position[0]-(0))-(1), (position[1]-(0))-(0));
                    if ((iSymbol[0] == iSymbol[1]) || (iSymbol[1] == 'N')) 
                    iTempResult[0] = iOne;
                    else 
                    iTempResult[0] = aQuality1[ (position[0])-1 ];
                    iEmission[0] = iTempResult[0];
                    if ((position[0]+-1<=0)&&(position[1]+-1>=0)) {
                        CurStateMemoryblock2From = dp.StateMemoryblock2.read((position[1]-(1))-(0));
                        CurStateMemoryblock3withbandingTo[2] = (((logspace::doubleexp( (((position[1])-1 - iStartMean) * ((position[1])-1 - iStartMean)) * iLocationFactor )))*(iEmission[0]))*CurStateMemoryblock2From[0];
                    }
                    if ((position[0]+-1>=1)&&(position[1]+-1>=0)) {
                        CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((position[0]-(1))-(1), (position[1]-(1))-(0));
                        CurStateMemoryblock3withbandingTo[2] += ((iTransition[12])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[0];
                        CurStateMemoryblock3withbandingTo[2] += ((iTransition[4])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[2];
                        CurStateMemoryblock3withbandingTo[2] += ((iTransition[8])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[1];
                    }
                    iEmission[0] = iOne;
                    if ((position[1]+-1>=0)) {
                        CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((position[0]-(0))-(1), (position[1]-(1))-(0));
                        CurStateMemoryblock3withbandingTo[0] = ((iTransition[13])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[0];
                        CurStateMemoryblock3withbandingTo[0] += ((iTransition[6])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[2];
                        CurStateMemoryblock3withbandingTo[0] += ((iTransition[10])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[1];
                    }
                    iEmission[0] = ((iEquilibrium > aHalfQuality1[ (position[0])-1 ]) ? iEquilibrium : aHalfQuality1[ (position[0])-1 ])  /* Was: iEquilibrium */;
                    if ((position[0]+-1<=0)) {
                        CurStateMemoryblock2From = dp.StateMemoryblock2.read((position[1]-(0))-(0));
                        CurStateMemoryblock3withbandingTo[1] = (((iGapOpen * logspace::doubleexp( (((position[1])-0 - iStartMean) * ((position[1])-0 - iStartMean)) * iLocationFactor )))*(iEmission[0]))*CurStateMemoryblock2From[0];
                    }
                    if ((position[0]+-1>=1)) {
                        CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((position[0]-(1))-(1), (position[1]-(0))-(0));
                        CurStateMemoryblock3withbandingTo[1] += ((iTransition[5])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[2];
                        CurStateMemoryblock3withbandingTo[1] += ((iTransition[9])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[1];
                    }
                    dp.StateMemoryblock3withbanding.written();
                }
                iPrevSlowCoord = position[1];
            } 
            else { 
                bandingInstance.warning();
                
            }
        } while (bandingInstance.hasNextForward());
    }
    iPrevSlowCoord = -1;
    for (int iPos1=0; iPos1<iLen2+1; ++iPos1) {
        {
            int iPos0=iLen1+0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            if (1) {
                if ((iPos1+-1>=0)) {
                    iSymbol[0] = pSequence2[iPos1+-1];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock4To = dp.StateMemoryblock4.write((iPos1-(0))-(0));
                iEmission[0] = iOne;
                if ((iPos1+-1>=0)) {
                    CurStateMemoryblock4From = dp.StateMemoryblock4.read((iPos1-(1))-(0));
                    CurStateMemoryblock4To[0] = ((iTransition[14])*(iEmission[0]))*CurStateMemoryblock4From[0];
                }
                iEmission[0] = iOne;
                if ((iPos0+0>=1)) {
                    CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((iPos0-(0))-(1), (iPos1-(0))-(0));
                    CurStateMemoryblock4To[0] += ((iTransition[7])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[2];
                    CurStateMemoryblock4To[0] += ((iTransition[11])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[1];
                }
                dp.StateMemoryblock4.written();
            }
            if ((iPos1+0>=iLen2+0)) {
                CurStateMemoryblock5To = dp.StateMemoryblock5.write();
                iEmission[0] = iOne;
                if (1) {
                    CurStateMemoryblock4From = dp.StateMemoryblock4.read((iPos1-(0))-(0));
                    CurStateMemoryblock5To[0] = ((iTransition[15])*(iEmission[0]))*CurStateMemoryblock4From[0];
                }
                dp.StateMemoryblock5.written();
            }
            iPrevSlowCoord = iPos1;
        }
    }
    iPrevSlowCoord = -1;
    {
        int iPos1=iLen2+0;
        if (iPos1==iPos1) {} // avoid 'unused variable' warnings
        {
            int iPos0=iLen1+0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            CurStateMemoryblock5From = dp.StateMemoryblock5.read();
            iTempProb[0] = CurStateMemoryblock5From[0];
        }
    }
    *ppOutTable = new AlignBandingDPTable(dp);
    // make sure tables don't get deleted
    dp.isInCharge = false;
    return iTempProb[0];
};





logspace BackwardBanding(AlignBandingDPTable** ppOutTable,char* pQuality1,char* pSequence1,char* pSequence2,int iGapExtendPhred,int iGapOpenPhred,int iInsNucPrior,int iLen1,int iLen2,int iStartMean,int iStartSD,int iWidth) {
    logspace iTransition[16];
    logspace *CurStateMemoryblock4To;
    const logspace *CurStateMemoryblock4From;
    const logspace *CurStateMemoryblock5From;
    logspace *CurStateMemoryblock3withbandingTo;
    const logspace *CurStateMemoryblock3withbandingFrom;
    logspace *CurStateMemoryblock2To;
    const logspace *CurStateMemoryblock2From;
    logspace *CurStateMemoryblock1To;
    const logspace *CurStateMemoryblock1From;
    int iPrevSlowCoord;
    /* initialization of various probabilities */
    
    static const double db = 10.0 / log(0.1);
    static const logspace iOne = 1.0;
    logspace iEquilibrium = exp( iInsNucPrior / db );     // likelihood of any unaligned read base
    logspace aQuality1[ iLen1 ];
    logspace aHalfQuality1[ iLen1 ];
    for (int i=0; i<iLen1; i++) {
        if (pQuality1[i] < (char)33) {
            cout << "Internal error: Negative quality at position "<<i<<endl;
        }
        aQuality1[i] = exp( (pQuality1[i]-33) / db );
        aHalfQuality1[i] = exp( (pQuality1[i]-33) / (2.0 * db) );
        if (pSequence1[i] == 'N') {
            aQuality1[i] = iOne;
            aHalfQuality1[i] = iOne;
        }
    }
    logspace iGapExtend = exp( iGapExtendPhred / db );
    logspace iGapOpen = exp( iGapOpenPhred / db ) * (iOne - iGapExtend);
    double iLocationFactor = -1.0 / (2.0*iStartSD*iStartSD);
    DiagonalBanding bandingInstance( iLen1, iLen2, iWidth );
    int iSymbol[2];
    if (false && iSymbol[0] == iSymbol[0]) {}   // avoid 'unused variable' warnings
    logspace iEmission[1];
    /* temporary storage for ordinary reals */
    register logspace iTempResult[1];
    /* temporary storage for extended-exponent reals */
    register logspace iTempProb[1];
    AlignBandingDPTable dp(iLen1,iLen2);
    iTransition[0] = iOne;
    
    iTransition[1] = iOne;
    
    iTransition[4] = iOne;
    
    iTransition[5] = iGapOpen;
    
    iTransition[6] = iGapOpen;
    
    iTransition[7] = iOne;
    
    iTransition[8] = iOne;
    
    iTransition[9] = iGapExtend;
    
    iTransition[10] = iGapOpen;
    
    iTransition[11] = iOne;
    
    iTransition[12] = iOne;
    
    iTransition[13] = iGapExtend;
    
    iTransition[14] = iOne;
    
    iTransition[15] = iOne;
    dp.StateMemoryblock5.write()[0] = 1.0;
    dp.StateMemoryblock5.written();
    iPrevSlowCoord = -1;
    for (int iPos1=(iLen2+1)-1; iPos1>=0; --iPos1) {
        {
            int iPos0=iLen1+0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            if ((iPos1+0>=iLen2+0)) {
            }
            if (1) {
                if ((iPos1+0<=iLen2+-1)) {
                    iSymbol[0] = pSequence2[iPos1+0];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock4To = dp.StateMemoryblock4.write((iPos1-(0))-(0));
                iEmission[0] = iOne;
                if ((iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock4From = dp.StateMemoryblock4.read((iPos1-(-1))-(0));
                    CurStateMemoryblock4To[0] = ((iTransition[14])*(iEmission[0]))*CurStateMemoryblock4From[0];
                }
                iEmission[0] = iOne;
                if ((iPos1+0>=iLen2+0)) {
                    CurStateMemoryblock5From = dp.StateMemoryblock5.read();
                    CurStateMemoryblock4To[0] += ((iTransition[15])*(iEmission[0]))*CurStateMemoryblock5From[0];
                }
                dp.StateMemoryblock4.written();
            }
            iPrevSlowCoord = iPos1;
        }
    }
    iPrevSlowCoord = -1;
    {
        Banding<2>::Position& position = bandingInstance.backwardIterator();
        int iCheckSlowCoordTraversal = -1;  
        do {
            if (iCheckSlowCoordTraversal != -1 && iCheckSlowCoordTraversal < position[1]) {
                cout << "WARNING: Banding (backward): Slowest coordinate be nonincreasing.  Perhaps forgot to specify speed of output coordinates?" << endl;
            }
            iCheckSlowCoordTraversal = position[1];
            if ((position[0]+0>=1)&&(position[0]+0<=iLen1+0)&&(position[1]+0>=0)&&(position[1]+0<=iLen2+0)) {
                if (1) {
                    if ((position[0]+0<=iLen1+-1)) {
                        iSymbol[0] = pSequence1[position[0]+0];
                    } 
                    else { 
                        iSymbol[0] = 'A' /* dummy value */;
                        
                    }
                    if ((position[1]+0<=iLen2+-1)) {
                        iSymbol[1] = pSequence2[position[1]+0];
                    } 
                    else { 
                        iSymbol[1] = 'A' /* dummy value */;
                        
                    }
                    CurStateMemoryblock3withbandingTo = dp.StateMemoryblock3withbanding.write((position[0]-(0))-(1), (position[1]-(0))-(0));
                    if ((iSymbol[0] == iSymbol[1]) || (iSymbol[1] == 'N')) 
                    iTempResult[0] = iOne;
                    else 
                    iTempResult[0] = aQuality1[ (position[0])-0 ];
                    iEmission[0] = iTempResult[0];
                    if ((position[0]+1<=iLen1+0)&&(position[1]+1<=iLen2+0)) {
                        CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((position[0]-(-1))-(1), (position[1]-(-1))-(0));
                        CurStateMemoryblock3withbandingTo[0] = ((iTransition[12])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[2];
                        CurStateMemoryblock3withbandingTo[2] = ((iTransition[4])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[2];
                        CurStateMemoryblock3withbandingTo[1] = ((iTransition[8])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[2];
                    }
                    iEmission[0] = iOne;
                    if ((position[1]+1<=iLen2+0)) {
                        CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((position[0]-(0))-(1), (position[1]-(-1))-(0));
                        CurStateMemoryblock3withbandingTo[0] += ((iTransition[13])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[0];
                        CurStateMemoryblock3withbandingTo[2] += ((iTransition[6])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[0];
                        CurStateMemoryblock3withbandingTo[1] += ((iTransition[10])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[0];
                    }
                    iEmission[0] = ((iEquilibrium > aHalfQuality1[ (position[0])-0 ]) ? iEquilibrium : aHalfQuality1[ (position[0])-0 ])  /* Was: iEquilibrium */;
                    if ((position[0]+1<=iLen1+0)) {
                        CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((position[0]-(-1))-(1), (position[1]-(0))-(0));
                        CurStateMemoryblock3withbandingTo[2] += ((iTransition[5])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[1];
                        CurStateMemoryblock3withbandingTo[1] += ((iTransition[9])*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[1];
                    }
                    iEmission[0] = iOne;
                    if ((position[0]+0>=iLen1+0)) {
                        CurStateMemoryblock4From = dp.StateMemoryblock4.read((position[1]-(0))-(0));
                        CurStateMemoryblock3withbandingTo[2] += ((iTransition[7])*(iEmission[0]))*CurStateMemoryblock4From[0];
                        CurStateMemoryblock3withbandingTo[1] += ((iTransition[11])*(iEmission[0]))*CurStateMemoryblock4From[0];
                    }
                    dp.StateMemoryblock3withbanding.written();
                }
                iPrevSlowCoord = position[1];
            } 
            else { 
                bandingInstance.warning();
                
            }
        } while (bandingInstance.hasNextBackward());
    }
    iPrevSlowCoord = -1;
    for (int iPos1=(iLen2+1)-1; iPos1>=0; --iPos1) {
        {
            int iPos0=0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            if (1) {
                if ((iPos0+0<=iLen1+-1)) {
                    iSymbol[0] = pSequence1[iPos0+0];
                } 
                else { 
                    iSymbol[0] = 'A' /* dummy value */;
                    
                }
                if ((iPos1+0<=iLen2+-1)) {
                    iSymbol[1] = pSequence2[iPos1+0];
                } 
                else { 
                    iSymbol[1] = 'A' /* dummy value */;
                    
                }
                CurStateMemoryblock2To = dp.StateMemoryblock2.write((iPos1-(0))-(0));
                if ((iSymbol[0] == iSymbol[1]) || (iSymbol[1] == 'N')) 
                iTempResult[0] = iOne;
                else 
                iTempResult[0] = aQuality1[ (iPos0)-0 ];
                iEmission[0] = iTempResult[0];
                if ((iPos0+1<=iLen1+0)&&(iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((iPos0-(-1))-(1), (iPos1-(-1))-(0));
                    CurStateMemoryblock2To[0] = (((logspace::doubleexp( (((iPos1)-1 - iStartMean) * ((iPos1)-1 - iStartMean)) * iLocationFactor )))*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[2];
                }
                iEmission[0] = iOne;
                if ((iPos1+1<=iLen2+0)) {
                    CurStateMemoryblock2From = dp.StateMemoryblock2.read((iPos1-(-1))-(0));
                    CurStateMemoryblock2To[0] += ((iTransition[1])*(iEmission[0]))*CurStateMemoryblock2From[0];
                }
                iEmission[0] = ((iEquilibrium > aHalfQuality1[ (iPos0)-0 ]) ? iEquilibrium : aHalfQuality1[ (iPos0)-0 ])  /* Was: iEquilibrium */;
                if ((iPos0+1<=iLen1+0)) {
                    CurStateMemoryblock3withbandingFrom = dp.StateMemoryblock3withbanding.read((iPos0-(-1))-(1), (iPos1-(0))-(0));
                    CurStateMemoryblock2To[0] += (((iGapOpen * logspace::doubleexp( (((iPos1)-0 - iStartMean) * ((iPos1)-0 - iStartMean)) * iLocationFactor )))*(iEmission[0]))*CurStateMemoryblock3withbandingFrom[1];
                }
                dp.StateMemoryblock2.written();
            }
            if ((iPos1+0<=0)) {
                CurStateMemoryblock1To = dp.StateMemoryblock1.write();
                iEmission[0] = iOne;
                if (1) {
                    CurStateMemoryblock2From = dp.StateMemoryblock2.read((iPos1-(0))-(0));
                    CurStateMemoryblock1To[0] = ((iTransition[0])*(iEmission[0]))*CurStateMemoryblock2From[0];
                }
                dp.StateMemoryblock1.written();
            }
            iPrevSlowCoord = iPos1;
        }
    }
    iPrevSlowCoord = -1;
    {
        int iPos1=0;
        if (iPos1==iPos1) {} // avoid 'unused variable' warnings
        {
            int iPos0=0;
            if (iPos0==iPos0) {} // avoid 'unused variable' warnings
            CurStateMemoryblock1From = dp.StateMemoryblock1.read();
            iTempProb[0] = CurStateMemoryblock1From[0];
        }
    }
    *ppOutTable = new AlignBandingDPTable(dp);
    // make sure tables don't get deleted
    dp.isInCharge = false;
    return iTempProb[0];
};



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