// latticearchive.cpp -- managing lattice archives
//
// F. Seide, V-hansu

#ifndef _CRT_SECURE_NO_WARNINGS
#define _CRT_SECURE_NO_WARNINGS // "secure" CRT not available on all platforms  --add this at the top of all CPP files that give "function or variable may be unsafe" warnings
#endif

#include "Basics.h"
#include "simple_checked_arrays.h"
#include "latticearchive.h"
#include "simplesenonehmm.h" // the model
#include "ssematrix.h"       // the matrices
#include "latticestorage.h"
#include <unordered_map>
#include <list>
#include <stdexcept>

using namespace std;

#define VIRGINLOGZERO (10 * LOGZERO) // used for printing statistics on unseen states
#undef CPU_VERIFICATION

#ifdef _WIN32
int msra::numa::node_override = -1; // for numahelpers.h
#endif

namespace msra { namespace lattices {

// ---------------------------------------------------------------------------
// helper class for allocation lots of small matrices, no free
// ---------------------------------------------------------------------------

class littlematrixheap
{
    static const size_t CHUNKSIZE;
    typedef msra::math::ssematrixfrombuffer matrixfrombuffer;
    std::list<std::vector<float>> heap;
    size_t allocatedinlast; // in last heap element
    size_t totalallocated;
    std::vector<matrixfrombuffer> matrices;

public:
    littlematrixheap(size_t estimatednumentries)
        : totalallocated(0), allocatedinlast(0)
    {
        matrices.reserve(estimatednumentries + 1);
    }
    msra::math::ssematrixbase &newmatrix(size_t rows, size_t cols)
    {
        const size_t elementsneeded = matrixfrombuffer::elementsneeded(rows, cols);
        if (heap.empty() || (heap.back().size() - allocatedinlast) < elementsneeded)
        {
            const size_t nelem = max(CHUNKSIZE, elementsneeded + 3 /*+3 for SSE alignment*/);
            heap.push_back(std::vector<float>(nelem));
            allocatedinlast = 0;
            // make sure starting element is SSE-aligned (the constructor demands that)
            const size_t offelem = (((size_t) &heap.back()[allocatedinlast]) / sizeof(float)) % 4;
            if (offelem != 0)
                allocatedinlast += 4 - offelem;
        }
        auto &buffer = heap.back();
        if (elementsneeded > heap.back().size() - allocatedinlast)
            LogicError("newmatrix: allocation logic screwed up");
        // get our buffer into a handy vector-like thingy
        array_ref<float> vecbuffer(&buffer[allocatedinlast], elementsneeded);
        // allocate in the current heap location
        matrices.resize(matrices.size() + 1);
        if (matrices.size() + 1 > matrices.capacity())
            LogicError("newmatrix: littlematrixheap cannot grow but was constructed with too small number of eements");
        auto &matrix = matrices.back();
        matrix = matrixfrombuffer(vecbuffer, rows, cols);
        allocatedinlast += elementsneeded;
        totalallocated += elementsneeded;
        return matrix;
    }
};

const size_t littlematrixheap::CHUNKSIZE = 256 * 1024; // 1 MB

// ---------------------------------------------------------------------------
// helpers for log-domain addition
// ---------------------------------------------------------------------------

#ifndef LOGZERO
#define LOGZERO -1e30f
#endif

// logadd (loga, logb) -> a += b, or loga = log [ exp(loga) + exp(logb) ]
static void logaddratio(float &loga, float diff)
{
    if (diff < -17.0f)
        return; // log (2^-24), 23-bit mantissa -> cut of after 24th bit
    loga += logf(1.0f + expf(diff));
}
static void logaddratio(double &loga, double diff)
{
    if (diff < -37.0f)
        return; // log (2^-53), 52-bit mantissa -> cut of after 53th bit
    loga += log(1.0 + exp(diff));
}
// loga <- log (exp (loga) + exp (logb)) = log (exp (loga) * (1.0 + exp (logb - loga)) = loga + log (1.0 + exp (logb - loga))
template <typename FLOAT>
static void logadd(FLOAT &loga, FLOAT logb)
{
    if (logb > loga) // we add smaller to bigger
        ::swap(loga, logb);
    if (loga <= LOGZERO) // both are 0
        return;
    logaddratio(loga, logb - loga);
}
template <typename FLOAT>
static void logmax(FLOAT &loga, FLOAT logb) // for testing (max approx)
{
    if (logb > loga)
        loga = logb;
}

template <typename FLOAT>
static FLOAT expdiff(FLOAT a, FLOAT b) // for testing
{
    if (b > a)
        return exp(b) * (exp(a - b) - 1);
    else
        return exp(a) * (1 - exp(b - a));
}

template <typename FLOAT>
static bool islogzero(FLOAT v)
{
    return v < LOGZERO / 2;
} // is this number to be considered 0

// ---------------------------------------------------------------------------
// other helpers go here
// ---------------------------------------------------------------------------

// helper to reconstruct the phonetic transcript
/*static*/ std::string lattice::gettranscript(const_array_ref<aligninfo> units, const msra::asr::simplesenonehmm &hset)
{
    std::string trans;
    foreach_index (k, units) // we exploit that units have fixed boundaries
    {
        if (k > 0)
            trans.push_back(' ');
        trans.append(hset.gethmm(units[k].unit).getname());
    }
    return trans;
}

// ---------------------------------------------------------------------------
// forwardbackwardedge() -- perform state-level forward-backward on a single lattice edge
//
// Results:
//  - gammas(j,t) for valid time ranges (remaining areas are not initialized)
//  - return value is edge acoustic score
// Gammas matrix must have two extra columns as buffer.
// ---------------------------------------------------------------------------

/*static*/ float lattice::forwardbackwardedge(const_array_ref<aligninfo> units, const msra::asr::simplesenonehmm &hset, const msra::math::ssematrixbase &logLLs,
                                              msra::math::ssematrixbase &loggammas, size_t edgeindex)
{
    // alphas and betas are stored in-place inside the loggammas matrix shifted by one?two columns
    assert(loggammas.cols() == logLLs.cols() + 2);
    msra::math::ssematrixstriperef<msra::math::ssematrixbase> logalphas(loggammas, 1, logLLs.cols()); // shifted views into gammas(,) for alphas and betas
    msra::math::ssematrixstriperef<msra::math::ssematrixbase> logbetas(loggammas, 2, logLLs.cols());

    // alphas(j,t) store the sum of all paths up to including state j at time t, including logLL(j,t)
    // betas(j,t) store the sum of all paths exiting from state j at time t, not including logLL(j,t)
    // gammas(j,t) = alphas(j,t) * betas(j,t) / totalLL

    // backward pass   --token passing
    size_t te = logbetas.cols();
    size_t je = logbetas.rows();
    float bwscore = 0.0f; // backward score
    for (size_t k = units.size() - 1; k + 1 > 0; k--)
    {
        const auto &hmm = hset.gethmm(units[k].unit);
        const size_t n = hmm.getnumstates();
        const auto &transP = hmm.gettransP();
        const size_t ts = te - units[k].frames; // end time of current unit
        const size_t js = je - n;               // range of state indices

        // pass in the transition score
        // t = ts: exit transition (last frame only or tee transition)
        float exitscore = 1e30f; // (something impossible)
        if (te == ts)            // tee transition
        {
            exitscore = bwscore + transP(-1, n);
        }
        else // not tee: expand all last states
        {
            for (size_t from = 0 /*no tee possible here*/; from < n; from++)
            {
                const size_t i = js + from; // origin trellis node
                logbetas(i, te - 1) = bwscore + transP(from, n);
            }
        }

        // expand from states j at time t (not yet including LL) to time t-1
        for (size_t t = te - 1; t + 1 > ts /*note: cannot test t >= ts because t < 0 possible*/; t--)
        {
            for (size_t to = 0; to < n; to++)
            {
                const size_t j = js + to;             // source trellis node
                const size_t s = hmm.getsenoneid(to); // senone id for state at position 'to' in the HMM
                const float acLL = logLLs(s, t);
                if (islogzero(acLL))
                    fprintf(stderr, "forwardbackwardedge: WARNING: edge J=%d unit %d (%s) frames [%d,%d) ac score(%d,%d) is zero (%d st, %d fr: %s)\n",
                            (int) edgeindex, (int) k, hmm.getname(), (int) ts, (int) te,
                            (int) s, (int) t,
                            (int) logbetas.rows(), (int) logbetas.cols(), gettranscript(units, hset).c_str());
                const float betajt = logbetas(j, t);   // sum over all all path exiting from (j,t) to end
                const float betajtpll = betajt + acLL; // incorporate acoustic score
                if (t > ts)
                    for (size_t from = 0 /*no transition from entry state*/; from < n; from++)
                    {
                        const size_t i = js + from; // target trellis node
                        const float pathscore = betajtpll + transP(from, to);
                        if (to == 0)
                            logbetas(i, t - 1 /*propagate into preceding frame*/) = pathscore;
                        else
                            logadd(logbetas(i, t - 1 /*propagate into preceding frame*/), pathscore);
                    }
                else // transition to entry state
                {
                    const float pathscore = betajtpll + transP(-1, to);
                    if (to == 0)
                        exitscore = pathscore;
                    else
                        logadd(exitscore, pathscore); // propagate into preceding unit
                }
            }
        }

        bwscore = exitscore;
        if (islogzero(bwscore))
            fprintf(stderr, "forwardbackwardedge: WARNING: edge J=%d unit %d (%s) frames [%d,%d) bw score is zero (%d st, %d fr: %s)\n",
                    (int) edgeindex, (int) k, hmm.getname(), (int) ts, (int) te, (int) logbetas.rows(), (int) logbetas.cols(), gettranscript(units, hset).c_str());

        te = ts;
        je = js;
    }
    assert(te == 0 && je == 0);
    const float totalbwscore = bwscore;

    // forward pass   --regular Viterbi
    // This also computes the gammas right away.
    size_t ts = 0;           // start frame for unit 'k'
    size_t js = 0;           // first row index of unit ' k'
    float fwscore = 0.0f;    // score passed across phone boundaries
    foreach_index (k, units) // we exploit that units have fixed boundaries
    {
        const auto &hmm = hset.gethmm(units[k].unit);
        const size_t n = hmm.getnumstates();
        const auto &transP = hmm.gettransP();
        const size_t te2 = ts + units[k].frames; // end time of current unit
        const size_t je2 = js + n;               // range of state indices

        // expand from states j at time t (including LL) to time t+1
        for (size_t t = ts; t < te2; t++) // note: loop not entered for 0-frame units (tees)
        {
            for (size_t to = 0; to < n; to++)
            {
                const size_t j = js + to; // target trellis node
                const size_t s = hmm.getsenoneid(to);
                const float acLL = logLLs(s, t);
                float alphajtnoll = LOGZERO;
                if (t == ts) // entering score
                {
                    const float pathscore = fwscore + transP(-1, to);
                    alphajtnoll = pathscore;
                }
                else
                    for (size_t from = 0 /*no entering possible*/; from < n; from++)
                    {
                        const size_t i = js + from; // origin trellis node
                        const float alphaitm1 = logalphas(i, t - 1 /*previous frame*/);
                        const float pathscore = alphaitm1 + transP(from, to);
                        logadd(alphajtnoll, pathscore);
                    }
                logalphas(j, t) = alphajtnoll + acLL;
            }
            // update the gammas  --do it here because in next frame, betas get overwritten by alphas (they share memory)
            for (size_t j = js; j < je2; j++)
            {
                if (!islogzero(totalbwscore))
                    loggammas(j, t) = logalphas(j, t) + logbetas(j, t) - totalbwscore;
                else // 0/0 problem, can occur if an ac score is so bad that it is 0 after going through softmax
                    loggammas(j, t) = LOGZERO;
            }
        }
        // t = te2: exit transition (last frame only or tee transition)
        float exitscore;
        if (te2 == ts) // tee transition
        {
            exitscore = fwscore + transP(-1, n);
        }
        else // not tee: expand all last states
        {
            exitscore = LOGZERO;
            for (size_t from = 0 /*no tee possible here*/; from < n; from++)
            {
                const size_t i = js + from;                   // origin trellis node
                const float alphaitm1 = logalphas(i, te2 - 1); // newly computed path score, transiting to t=te2
                const float pathscore = alphaitm1 + transP(from, n);
                logadd(exitscore, pathscore);
            }
        }
        fwscore = exitscore; // score passed on to next unit
        js = je2;
        ts = te2;
    }
    assert(js == logalphas.rows() && ts == logalphas.cols());
    const float totalfwscore = fwscore;

    // in extreme cases, we may have 0 ac probs, which lead to 0 path scores and division by 0 (subtracting LOGZERO)
    // These cases must be handled separately. If the whole path is 0 (0 prob is on the only path at some point) then skip the lattice.
    if (islogzero(totalbwscore) ^ islogzero(totalfwscore))
        fprintf(stderr, "forwardbackwardedge: WARNING: edge J=%d fw and bw 0 score %.10f vs. %.10f (%d st, %d fr: %s)\n",
                (int) edgeindex, (float) totalfwscore, (float) totalbwscore, (int) js, (int) ts, gettranscript(units, hset).c_str());
    if (islogzero(totalbwscore))
    {
        fprintf(stderr, "forwardbackwardedge: WARNING: edge J=%d has zero ac. score (%d st, %d fr: %s)\n",
                (int) edgeindex, (int) js, (int) ts, gettranscript(units, hset).c_str());
        return LOGZERO;
    }

    if (fabsf(totalfwscore - totalbwscore) / ts > 1e-4f)
        fprintf(stderr, "forwardbackwardedge: WARNING: edge J=%d fw and bw score %.10f vs. %.10f (%d st, %d fr: %s)\n",
                (int) edgeindex, (float) totalfwscore, (float) totalbwscore, (int) js, (int) ts, gettranscript(units, hset).c_str());

    // we return the full path score
    return totalfwscore;
}

// ---------------------------------------------------------------------------
// alignedge() -- perform Viterbi alignment on a single edge
//
// This is an alternative to forwardbackwardedge() that just uses the best path.
// Results:
//  - if not returnsenoneids -> 'binary gammas(j,t)' for valid time ranges (remaining areas are not initialized); MMI-compatible
//  - if returnsenoneids ->  loggammas(0,t) will contain the senone ids directly instead (for sMBR mode)
//  - return value is edge acoustic score
// Gammas matrix must have two extra columns as buffer.
// ---------------------------------------------------------------------------

/*static*/ float lattice::alignedge(const_array_ref<aligninfo> units, const msra::asr::simplesenonehmm &hset, const msra::math::ssematrixbase &logLLs,
                                    msra::math::ssematrixbase &loggammas, size_t edgeindex /*for diagnostic messages*/, const bool returnsenoneids,
                                    array_ref<unsigned short> thisedgealignmentsj)
{
    // alphas and betas are stored in-place inside the loggammas matrix shifted by one?two columns
    assert(loggammas.cols() == logLLs.cols() + 2);
    msra::math::ssematrixstriperef<msra::math::ssematrixbase> backpointers(loggammas, 0, logLLs.cols());
    msra::math::ssematrixstriperef<msra::math::ssematrixbase> pathscores(loggammas, 2, logLLs.cols());

    // pathscores(j,t) store the sum of all paths up to including state j at time t, including logLL(j,t)
    // backpointers(j,t) are the relative states that it came from
    // gammas(j,t) <- 1 if on best path, 0 otherwise

    const int invalidbp = -2;

    // Viterbi alignment
    size_t ts = 0;           // start frame for unit 'k'
    size_t js = 0;           // first row index of unit 'k'
    float fwscore = 0.0f;    // score passed across phone boundaries
    int fwbackpointer = -1;  // bp passed across phone boundaries, -1 means start of utterance
    foreach_index (k, units) // we exploit that units have fixed boundaries
    {
        const auto &hmm = hset.gethmm(units[k].unit);
        const size_t n = hmm.getnumstates();
        const auto &transP = hmm.gettransP();
        const size_t te = ts + units[k].frames;    // end time of current unit
        const size_t je = js + hmm.getnumstates(); // range of state indices

        // expand from states j at time t (including LL) to time t+1
        for (size_t t = ts; t < te; t++) // note: loop not entered for 0-frame units (tees)
        {
            for (size_t j = js; j < je; j++)
            {
                const size_t to = j - js; // relative state
                const size_t s = hmm.getsenoneid(to);
                pathscores(j, t) = LOGZERO;
                backpointers(j, t) = invalidbp;
                if (t == ts) // entering score
                {
                    const float pathscore = fwscore + transP(-1, to);
                    pathscores(j, t) = pathscore;
                    backpointers(j, t) = (float) fwbackpointer;
                }
                else
                    for (size_t i = js; i < je; i++)
                    {
                        const size_t from = i - js;
                        const float alphaitm1 = pathscores(i, t - 1 /*previous frame*/);
                        const float pathscore = alphaitm1 + transP(from, to);
                        if (pathscore > pathscores(j, t))
                        {
                            pathscores(j, t) = pathscore;
                            backpointers(j, t) = (float) i;
                        }
                    }
                const float acLL = logLLs(s, t);
                pathscores(j, t) += acLL;
            }
        }
        // t = te: exit transition (last frame only or tee transition)
        float exitscore = LOGZERO;
        int exitbackpointer = invalidbp;
        if (te == ts) // tee transition
        {
            exitscore = fwscore + transP(-1, n);
            exitbackpointer = fwbackpointer;
        }
        else // not tee: expand all last states
        {
            for (size_t i = js; i < je; i++)
            {
                const size_t from = i - js;
                const float alphaitm1 = pathscores(i, te - 1); // newly computed path score, transiting to t=te
                const float pathscore = alphaitm1 + transP(from, n);
                if (pathscore > exitscore)
                {
                    exitscore = pathscore;
                    exitbackpointer = (int) i;
                }
            }
        }
        if (exitbackpointer == invalidbp)
            LogicError("exitbackpointer came up empty");
        fwscore = exitscore;             // score passed on to next unit
        fwbackpointer = exitbackpointer; // and accompanying backpointer
        js = je;
        ts = te;
    }
    assert(js == pathscores.rows() && ts == pathscores.cols());

    // in extreme cases, we may have 0 ac probs, which lead to 0 path scores and division by 0 (subtracting LOGZERO)
    // These cases must be handled separately. If the whole path is 0 (0 prob is on the only path at some point) then skip the lattice.
    if (islogzero(fwscore))
    {
        fprintf(stderr, "alignedge: WARNING: edge J=%d has zero ac. score (%d st, %d fr: %s)\n",
                (int) edgeindex, (int) js, (int) ts, gettranscript(units, hset).c_str());
        return LOGZERO;
    }

    // traceback & gamma update
    size_t te = backpointers.cols();
    size_t je = backpointers.rows();
    int j = fwbackpointer;
    for (size_t k = units.size() - 1; k + 1 > 0; k--) // go in units because we also need to clear out the column
    {
        const auto &hmm = hset.gethmm(units[k].unit);
        const size_t ts2 = te - units[k].frames;    // end time of current unit
        const size_t js2 = je - hmm.getnumstates(); // range of state indices
        for (size_t t = te - 1; t + 1 > ts2; t--)
        {
            if (j < (int)js2 || j >= (int) je)
                LogicError("invalid backpointer resulting in state index out of range");

            int bp = (int) backpointers(j, t); // save the backpointer before overwriting it (gammas and backpointers are aliases of each other)
                                               // thisedgealignmentsj[t] = (unsigned short)hmm.getsenoneid(j - js2);
            if (!returnsenoneids)              // return binary gammas (for MMI; this mode is compatible with softalignmode)
                for (size_t i = js2; i < je; i++)
                    loggammas(i, t) = ((int) i == j) ? 0.0f : LOGZERO;
            else // return senone id (for sMBR; note: NOT compatible with softalignmode; calling code must know this)
                thisedgealignmentsj[t] = (unsigned short) hmm.getsenoneid(j - js2);

            if (bp == invalidbp)
                LogicError("deltabackpointer not initialized");
            j = bp; // trace back one step
        }

        te = ts2;
        je = js2;
    }
    if (j != -1)
        LogicError("invalid backpointer resulting in not reaching start of utterance when tracing back");
    assert(je == 0 && te == 0);

    // we return the full path score
    return fwscore;
}

// ---------------------------------------------------------------------------
// forwardbackwardlattice() -- lattice-level forward/backward
//
// This computes word posteriors, and also returns the per-node alphas and betas.
// Per-edge acoustic scores are passed in via a lambda, as this function is
// intended for use at multiple places with different scores.
// (Specifically, we also use it to determine a pruning threshold, based on
// the original lattice's ac. scores, before even bothering to compute the
// new ac. scores.)
// ---------------------------------------------------------------------------

double lattice::forwardbackwardlattice(const std::vector<float> &edgeacscores, parallelstate &parallelstate, std::vector<double> &logpps,
                                       std::vector<double> &logalphas, std::vector<double> &logbetas,
                                       const float lmf, const float wp, const float amf, const float boostingfactor, const bool sMBRmode,
                                       const_array_ref<size_t> &uids, const edgealignments &thisedgealignments,
                                       std::vector<double> &logEframescorrect, std::vector<double> &Eframescorrectbuf, double &logEframescorrecttotal) const
{ // ^^ TODO: remove this
    // --- hand off to parallelized (CUDA) implementation if available
    if (parallelstate.enabled())
    {
        double totalfwscore = parallelforwardbackwardlattice(parallelstate, edgeacscores, thisedgealignments, lmf, wp, amf, boostingfactor, logpps, logalphas, logbetas, sMBRmode, uids, logEframescorrect, Eframescorrectbuf, logEframescorrecttotal);

        return totalfwscore;
    }
    // if we get here, we have no CUDA, and do it the good ol' way

    // allocate return values
    logpps.resize(edges.size()); // this is our primary return value

    // TODO: these are return values as well, but really shouldn't anymore; only used in some older baseline code we some day may want to compare against
    logalphas.assign(nodes.size(), LOGZERO);
    logalphas.front() = 0.0f;
    logbetas.assign(nodes.size(), LOGZERO);
    logbetas.back() = 0.0f;

    // --- sMBR version

    if (sMBRmode)
    {
        logEframescorrect.resize(edges.size());
        Eframescorrectbuf.resize(edges.size());

        std::vector<double> logaccalphas(nodes.size(), LOGZERO); // [i] expected frames-correct count over all paths from start to node i
        std::vector<double> logaccbetas(nodes.size(), LOGZERO);  // [i] likewise
        std::vector<double> logframescorrectedge(edges.size());  // raw counts of correct frames in each edge

        // forward pass
        foreach_index (j, edges)
        {
            if (islogzero(edgeacscores[j])) // indicates that this edge is pruned
                continue;
            const auto &e = edges[j];
            const double inscore = logalphas[e.S];
            const double edgescore = (e.l * lmf + wp + edgeacscores[j]) / amf;
            const double pathscore = inscore + edgescore;
            logadd(logalphas[e.E], pathscore);

            size_t ts = nodes[e.S].t;
            size_t te = nodes[e.E].t;
            size_t framescorrect = 0; // count raw number of correct frames
            for (size_t t = ts; t < te; t++)
                framescorrect += (thisedgealignments[j][t - ts] == uids[t]);
            logframescorrectedge[j] = (framescorrect > 0) ? log((double) framescorrect) : LOGZERO; // remember for backward pass
            double loginaccs = logaccalphas[e.S] - logalphas[e.S];
            logadd(loginaccs, logframescorrectedge[j]);
            double logpathacc = loginaccs + logalphas[e.S] + edgescore;
            logadd(logaccalphas[e.E], logpathacc);
        }
        foreach_index (j, logaccalphas)
            logaccalphas[j] -= logalphas[j];

        const double totalfwscore = logalphas.back();
        const double totalfwacc = logaccalphas.back();
        if (islogzero(totalfwscore))
        {
            fprintf(stderr, "forwardbackward: WARNING: no path found in lattice (%d nodes/%d edges)\n", (int) nodes.size(), (int) edges.size());
            return LOGZERO; // failed, do not use resulting matrix
        }

        // backward pass and computation of state-conditioned frames-correct count
        for (size_t j = edges.size() - 1; j + 1 > 0; j--)
        {
            if (islogzero(edgeacscores[j])) // indicates that this edge is pruned
                continue;
            const auto &e = edges[j];
            const double inscore = logbetas[e.E];
            const double edgescore = (e.l * lmf + wp + edgeacscores[j]) / amf;
            const double pathscore = inscore + edgescore;
            logadd(logbetas[e.S], pathscore);

            double loginaccs = logaccbetas[e.E] - logbetas[e.E];
            logadd(loginaccs, logframescorrectedge[j]);
            double logpathacc = loginaccs + logbetas[e.E] + edgescore;
            logadd(logaccbetas[e.S], logpathacc);

            // sum up to get final expected frames-correct count per state == per edge (since we assume hard state alignment)
            double logpp = logalphas[e.S] + edgescore + logbetas[e.E] - totalfwscore;
            if (logpp > 1e-2)
                fprintf(stderr, "forwardbackward: WARNING: edge J=%d log posterior %.10f > 0\n", (int) j, (float) logpp);
            if (logpp > 0.0)
                logpp = 0.0;
            logpps[j] = logpp;
            double tmplogeframecorrect = logframescorrectedge[j];
            logadd(tmplogeframecorrect, logaccalphas[e.S]);
            logadd(tmplogeframecorrect, logaccbetas[e.E] - logbetas[e.E]);
            Eframescorrectbuf[j] = exp(tmplogeframecorrect);
        }
        foreach_index (j, logaccbetas)
            logaccbetas[j] -= logbetas[j];
        const double totalbwscore = logbetas.front();
        const double totalbwacc = logaccbetas.front();
        if (fabs(totalfwscore - totalbwscore) / info.numframes > 1e-4)
            fprintf(stderr, "forwardbackward: WARNING: lattice fw and bw scores %.10f vs. %.10f (%d nodes/%d edges)\n", (float) totalfwscore, (float) totalbwscore, (int) nodes.size(), (int) edges.size());

        if (fabs(totalfwacc - totalbwacc) / info.numframes > 1e-4)
            fprintf(stderr, "forwardbackwardlatticesMBR: WARNING: lattice fw and bw accs %.10f vs. %.10f (%d nodes/%d edges)\n", (float) totalfwacc, (float) totalbwacc, (int) nodes.size(), (int) edges.size());

        logEframescorrecttotal = totalbwacc;
        return totalbwscore;
    }

    // --- MMI version

    // forward pass
    foreach_index (j, edges)
    {
        const auto &e = edges[j];
        const double inscore = logalphas[e.S];
        const double edgescore = (e.l * lmf + wp + edgeacscores[j]) / amf; // note: edgeacscores[j] == LOGZERO if edge was pruned
        const double pathscore = inscore + edgescore;
        logadd(logalphas[e.E], pathscore);
    }
    const double totalfwscore = logalphas.back();
    if (islogzero(totalfwscore))
    {
        fprintf(stderr, "forwardbackward: WARNING: no path found in lattice (%d nodes/%d edges)\n", (int) nodes.size(), (int) edges.size());
        return LOGZERO; // failed, do not use resulting matrix
    }

    // backward pass
    // this also computes the word posteriors on the fly, since we are at it
    for (size_t j = edges.size() - 1; j + 1 > 0; j--)
    {
        const auto &e = edges[j];
        const double inscore = logbetas[e.E];
        const double edgescore = (e.l * lmf + wp + edgeacscores[j]) / amf;
        const double pathscore = inscore + edgescore;
        logadd(logbetas[e.S], pathscore);

        // compute lattice posteriors on the fly since we are at it
        double logpp = logalphas[e.S] + edgescore + logbetas[e.E] - totalfwscore;
        if (logpp > 1e-2)
            fprintf(stderr, "forwardbackward: WARNING: edge J=%d log posterior %.10f > 0\n", (int) j, (float) logpp);
        if (logpp > 0.0)
            logpp = 0.0;
        logpps[j] = logpp;
    }

    const double totalbwscore = logbetas.front();
    if (fabs(totalfwscore - totalbwscore) / info.numframes > 1e-4)
        fprintf(stderr, "forwardbackward: WARNING: lattice fw and bw scores %.10f vs. %.10f (%d nodes/%d edges)\n", (float) totalfwscore, (float) totalbwscore, (int) nodes.size(), (int) edges.size());

    return totalfwscore;
}

// ---------------------------------------------------------------------------
// forwardbackwardlatticesMBR() -- compute expected frame-accuracy counts,
// both the conditioned one (corresponding to c(q) in Dan Povey's thesis)
// and the global one (which is the sMBR criterion to optimize).
//
// Outputs:
//  - Eframescorrect[j] == expected frames-correct count conditioned on a state of edge[j].
//    We currently assume a hard state alignment. With that, the value turns out
//    to be identical for all states of an edge, so we only store it once per edge.
//  - return value: expected frames-correct count for entire lattice
//
// Call forwardbackwardlattices() first to compute logalphas/betas.
// ---------------------------------------------------------------------------

double lattice::forwardbackwardlatticesMBR(const std::vector<float> &edgeacscores, const msra::asr::simplesenonehmm &hset,
                                           const std::vector<double> &logalphas, const std::vector<double> &logbetas,
                                           const float lmf, const float wp, const float amf, const_array_ref<size_t> &uids,
                                           const edgealignments &thisedgealignments, std::vector<double> &Eframescorrect) const
{
    std::vector<double> accalphas(nodes.size(), 0);  // [i] expected frames-correct count over all paths from start to node i
    std::vector<double> accbetas(nodes.size(), 0);   // [i] likewise
    std::vector<size_t> maxcorrect(nodes.size(), 0); // [i] max correct frames up to this node (oracle)

    std::vector<double> framescorrectedge(edges.size()); // raw counts of correct frames in each edge

    std::vector<int> backpointersformaxcorr(nodes.size(), -2); // keep track of backpointer for the max corr
    backpointersformaxcorr.front() = -1;

    // forward pass
    foreach_index (j, edges)
    {
        if (islogzero(edgeacscores[j])) // indicates that this edge is pruned
            continue;
        const auto &e = edges[j];
        const double inaccs = accalphas[e.S];
        size_t ts = nodes[e.S].t;
        size_t te = nodes[e.E].t;

        size_t framescorrect = 0; // count raw number of correct frames
        for (size_t t = ts; t < te; t++)
            framescorrect += (thisedgealignments[j][t - ts] == uids[t]);
        framescorrectedge[j] = (double) framescorrect; // remember for backward pass

        const double edgescore = (e.l * lmf + wp + edgeacscores[j]) / amf;
        // contribution to end node's path acc = start node's plus edge's correct count, weighted by LL, and divided by sum over LLs
        double pathacc = (inaccs + framescorrectedge[j]) * exp(logalphas[e.S] + edgescore - logalphas[e.E]);
        accalphas[e.E] += pathacc;
        // also keep track of max accuracy, so we can find out whether the lattice contains the correct path
        size_t oracleframescorrect = maxcorrect[e.S] + framescorrect; // keep track of most correct path up to end of this edge
        if (oracleframescorrect > maxcorrect[e.E])
        {
            maxcorrect[e.E] = oracleframescorrect;
            backpointersformaxcorr[size_t(e.E)] = j;
        }
    }
    const double totalfwacc = accalphas.back();

    hset; // just for reference

    // report on ground-truth path
    // TODO: we will later have code that adds this path if needed
    size_t oracleframeacc = maxcorrect.back();
    if (oracleframeacc != info.numframes)
        fprintf(stderr, "forwardbackwardlatticesMBR: ground-truth path missing from lattice (most correct path: %d out of %d frames correct)\n", (unsigned int) oracleframeacc, (int) info.numframes);

    // backward pass and computation of state-conditioned frames-correct count
    for (size_t j = edges.size() - 1; j + 1 > 0; j--)
    {
        if (islogzero(edgeacscores[j])) // indicates that this edge is pruned
            continue;
        const auto &e = edges[j];
        const double inaccs = accbetas[e.E];
        const double edgescore = (e.l * lmf + wp + edgeacscores[j]) / amf;
        double pathacc = (inaccs + framescorrectedge[j]) * exp(logbetas[e.E] + edgescore - logbetas[e.S]);
        accbetas[e.S] += pathacc;

        // sum up to get final expected frames-correct count per state == per edge (since we assume hard state alignment)
        Eframescorrect[j] = (float) (accalphas[e.S] + accbetas[e.E] + framescorrectedge[j]);
    }

    const double totalbwacc = accbetas.front();

    if (fabs(totalfwacc - totalbwacc) / info.numframes > 1e-4)
        fprintf(stderr, "forwardbackwardlatticesMBR: WARNING: lattice fw and bw accs %.10f vs. %.10f (%d nodes/%d edges)\n", (float) totalfwacc, (float) totalbwacc, (int) nodes.size(), (int) edges.size());

    return totalbwacc;
}

// ---------------------------------------------------------------------------
// bestpathlattice() -- lattice-level "forward/backward" that only returns the
// best path, but in the form of word posteriors, which are 1 or 0, just like
// a real lattice-level forward/backward would do.
// We don't really use this; this was only for a contrast experiment.
// ---------------------------------------------------------------------------

double lattice::bestpathlattice(const std::vector<float> &edgeacscores, std::vector<double> &logpps,
                                const float lmf, const float wp, const float amf) const
{
    // forward pass --sortnedness => regular Viterbi
    std::vector<double> logalphas(nodes.size(), LOGZERO);
    std::vector<int> backpointers(nodes.size(), -2);
    logalphas.front() = 0.0f;
    backpointers.front() = -1;
    foreach_index (j, edges)
    {
        const auto &e = edges[j];
        const double inscore = logalphas[e.S];
        const double edgescore = (e.l * lmf + wp + edgeacscores[j]) / amf; // note: edgeacscores[j] == LOGZERO if edge was pruned
        const double pathscore = inscore + edgescore;
        if (pathscore > logalphas[e.E])
        {
            logalphas[e.E] = pathscore;
            backpointers[e.E] = j;
        }
    }

    const double totalfwscore = logalphas.back();
    if (islogzero(totalfwscore))
    {
        fprintf(stderr, "bestpathlattice: WARNING: no path found in lattice (%d nodes/%d edges)\n", (int) nodes.size(), (int) edges.size());
        return LOGZERO; // failed, do not use resulting matrix
    }

    // traceback
    // We encode the result by storing log 1 in edges on the best path, and log 0 else;
    // this makes it naturally compatible with softalign mode
    logpps.resize(edges.size());
    foreach_index (j, edges)
        logpps[j] = LOGZERO;

    int backpos = backpointers[nodes.size() - 1];
    while (backpos >= 0)
    {
        logpps[backpos] = 0.0f; // edge is on best path -> PP = 1.0
        backpos = backpointers[edges[backpos].S];
    }
    assert(backpos == -1);

    return totalfwscore;
}

// ---------------------------------------------------------------------------
// forwardbackwardalign() -- compute the statelevel gammas or viterbi alignments
// the first phase of lattice::forwardbackward
//
// Outputs:
// ---------------------------------------------------------------------------
void lattice::forwardbackwardalign(parallelstate &parallelstate,
                                   const msra::asr::simplesenonehmm &hset, const bool softalignstates,
                                   const double minlogpp, const std::vector<double> &origlogpps,
                                   std::vector<msra::math::ssematrixbase *> &abcs, littlematrixheap &matrixheap,
                                   const bool returnsenoneids,
                                   std::vector<float> &edgeacscores, const msra::math::ssematrixbase &logLLs,
                                   edgealignments &thisedgealignments, backpointers &thisbackpointers, array_ref<size_t> &uids, const_array_ref<size_t> bounds) const
{ // NOTE: this will be removed and replaced by a proper representation of alignments someday
    // do forward-backward or alignment on a per-edge basis. This gives us:
    //  - per-edge gamma[j,t] = P(s(t)==s_j|edge) if forwardbackward, per-edge alignment thisedgealignments[j] if alignment
    //  - per-edge acoustic scores
    const size_t silunitid = hset.gethmmid("sil"); // shall be the same as parallelstate.getsilunitid()
    bool parallelsil = true;
    bool cpuverification = false;

#ifndef PARALLEL_SIL // we use a define to make this marked
    parallelsil = false;
#endif
#ifdef CPU_VERIFICATION
    cpuverification = true;
#endif

    // Phase 1: abcs allocate
    if (!parallelstate.enabled() || !parallelsil || cpuverification) // allocate abcs when 1.parallelstate not enabled (cpu mode); 2. enabled but not PARALLEL_SIL (silence need to be allocate); 3. cpuverfication
    {
        abcs.resize(edges.size(), NULL); // [edge index] -> alpha/beta/gamma matrices for each edge
        size_t countskip = 0;            // if pruning: count how many edges are pruned

        foreach_index (j, edges)
        {
            // determine number of frames
            // TODO: this is not efficient--we only use a block-diagonal-like structure, rest is empty (exploiting the fixed boundaries)
            const size_t edgeframes = nodes[edges[j].E].t - nodes[edges[j].S].t;
            if (edgeframes == 0) // dummy !NULL edge at end of lattice
            {
                if ((size_t) j != edges.size() - 1)
                    RuntimeError("forwardbackwardalign: unxpected 0-frame edge (only allowed at very end)");
                // note: abcs[j] is already initialized to be NULL in this case, which protects us from accidentally using it
            }
            else
            {
                // determine the number of states in an edge
                const auto &aligntokens = getaligninfo(j); // get alignment tokens
                size_t edgestates = 0;

                bool edgehassil = false;
                foreach_index (i, aligntokens)
                    if (aligntokens[i].unit == silunitid)
                        edgehassil = true;

                if (!cpuverification && !edgehassil && parallelstate.enabled()) // !cpuverification, parallel & is non sil, we do not allocate
                {
                    abcs[j] = NULL;
                    continue;
                }

                foreach_index (k, aligntokens)
                    edgestates += hset.gethmm(aligntokens[k].unit).getnumstates();

                // allocate the matrix
                if (minlogpp > LOGZERO && origlogpps[j] < minlogpp)
                    countskip++;
                else
                    abcs[j] = &matrixheap.newmatrix(edgestates, edgeframes + 2); // +2 to have one extra column for betas and one for gammas
            }
        }
        if (minlogpp > LOGZERO)
            fprintf(stderr, "forwardbackwardalign: %d of %d edges pruned\n", (int) countskip, (int) edges.size());
    }

    // Phase 2: alignment on CPU
    if (parallelstate.enabled() && !parallelsil) // silence edge shall be process separately if not cuda and not PARALLEL_SIL
    {
        if (softalignstates)
            LogicError("forwardbackwardalign: parallelized version currently only handles hard alignments");
        if (minlogpp > LOGZERO)
            fprintf(stderr, "forwardbackwardalign: pruning not supported (we won't need it!) :)\n");
        edgeacscores.resize(edges.size());
        for (size_t j = 0; j < edges.size(); j++)
        {
            const auto &aligntokens = getaligninfo(j); // get alignment tokens
            if (aligntokens.size() == 0)
                continue;
            bool edgehassil = false;
            foreach_index (i, aligntokens)
            {
                if (aligntokens[i].unit == silunitid)
                    edgehassil = true;
            }
            if (!edgehassil) // only process sil
                continue;
            const edgeinfowithscores &e = edges[j];
            const size_t ts = nodes[e.S].t;
            const size_t te = nodes[e.E].t;
            const auto edgeLLs = msra::math::ssematrixstriperef<msra::math::ssematrixbase>(const_cast<msra::math::ssematrixbase &>(logLLs), ts, te - ts);
            edgeacscores[j] = alignedge(aligntokens, hset, edgeLLs, *abcs[j], j, true, thisedgealignments[j]);
        }
    }

    // Phase 3: alignment on GPU
    if (parallelstate.enabled())
        parallelforwardbackwardalign(parallelstate, hset, logLLs, edgeacscores, thisedgealignments, thisbackpointers);

    // zhaorui align to reference mlf
    if (bounds.size() > 0)
    {
        size_t framenum = bounds.size();

        msra::math::ssematrixbase *refabcs;
        size_t ts, te, t;
        ts = te = 0;

        vector<aligninfo> refinfo(1);
        vector<unsigned short> refalign(framenum);

        array_ref<aligninfo> refunits(refinfo.data(), 1);
        array_ref<unsigned short> refedgealignmentsj(refalign.data(), framenum);

        while (te < framenum)
        {
            // found one phone's boundary (ts, te)
            t = ts + 1;
            while (t < framenum && bounds[t] == 0)
                t++;
            te = t;

            // make one phone unit
            size_t phoneid = bounds[ts] - 1;
            refunits[0].unit = phoneid;
            refunits[0].frames = te - ts;

            size_t edgestates = hset.gethmm(phoneid).getnumstates();
            littlematrixheap refmatrixheap(1); // for abcs
            refabcs = &refmatrixheap.newmatrix(edgestates, te - ts + 2);
            const auto edgeLLs = msra::math::ssematrixstriperef<msra::math::ssematrixbase>(const_cast<msra::math::ssematrixbase &>(logLLs), ts, te - ts);
            // do alignment
            alignedge((const_array_ref<aligninfo>) refunits, hset, edgeLLs, *refabcs, 0, true, refedgealignmentsj);

            for (t = ts; t < te; t++)
            {
                uids[t] = (size_t) refedgealignmentsj[t - ts];
            }
            ts = te;
        }
    }

    // Phase 4: alignment or forwardbackward on CPU for non parallel mode or verification

    if (!parallelstate.enabled() || cpuverification) // non parallel mode or verification
    {
        edgeacscores.resize(edges.size());
        std::vector<float> edgeacscoresgpu;
        edgealignments thisedgealignmentsgpu(thisedgealignments);
        if (cpuverification)
        {
            parallelstate.getedgeacscores(edgeacscoresgpu);
            parallelstate.copyalignments(thisedgealignmentsgpu);
        }
        foreach_index (j, edges)
        {
            const edgeinfowithscores &e = edges[j];
            const size_t ts = nodes[e.S].t;
            const size_t te = nodes[e.E].t;
            if (ts == te) // dummy !NULL edge at end
                edgeacscores[j] = 0.0f;
            else
            {
                const auto &aligntokens = getaligninfo(j); // get alignment tokens
                const auto edgeLLs = msra::math::ssematrixstriperef<msra::math::ssematrixbase>(const_cast<msra::math::ssematrixbase &>(logLLs), ts, te - ts);
                if (minlogpp > LOGZERO && origlogpps[j] < minlogpp)
                    edgeacscores[j] = LOGZERO; // will kill word level forwardbackward hypothesis
                else if (softalignstates)
                    edgeacscores[j] = forwardbackwardedge(aligntokens, hset, edgeLLs, *abcs[j], j);
                else
                    edgeacscores[j] = alignedge(aligntokens, hset, edgeLLs, *abcs[j], j, returnsenoneids, thisedgealignments[j]);
            }
            if (cpuverification)
            {
                const auto &aligntokens = getaligninfo(j); // get alignment tokens
                bool edgehassil = false;
                foreach_index (i, aligntokens)
                {
                    if (aligntokens[i].unit == silunitid)
                        edgehassil = true;
                }
                if (fabs(edgeacscores[j] - edgeacscoresgpu[j]) > 1e-3)
                {
                    fprintf(stderr, "edge %d, sil ? %d, edgeacscores / edgeacscoresgpu MISMATCH %f v.s. %f, diff %e\n",
                            j, edgehassil ? 1 : 0, (float) edgeacscores[j], (float) edgeacscoresgpu[j],
                            (float) (edgeacscores[j] - edgeacscoresgpu[j]));
                    fprintf(stderr, "aligntokens: ");
                    foreach_index (i, aligntokens)
                        fprintf(stderr, "%d %d; ", i, aligntokens[i].unit);
                    fprintf(stderr, "\n");
                }
                for (size_t t = ts; t < te; t++)
                {
                    if (thisedgealignments[j][t - ts] != thisedgealignmentsgpu[j][t - ts])
                        fprintf(stderr, "edge %d, sil ? %d, time %d, alignment / alignmentgpu MISMATCH %d v.s. %d\n", j, edgehassil ? 1 : 0, (int) (t - ts), thisedgealignments[j][t - ts], thisedgealignmentsgpu[j][t - ts]);
                }
            }
        }
    }
}

// compute the error signal for sMBR mode
void lattice::sMBRerrorsignal(parallelstate &parallelstate,
                              msra::math::ssematrixbase &errorsignal, msra::math::ssematrixbase &errorsignalneg, // output
                              const std::vector<double> &logpps, const float amf,
                              double minlogpp, const std::vector<double> &origlogpps, const std::vector<double> &logEframescorrect,
                              const double logEframescorrecttotal, const edgealignments &thisedgealignments) const
{
    if (parallelstate.enabled()) // parallel version
    {
        /*  time measurement for parallel sMBRerrorsignal
            errorsignalcompute: 19.871935 ms (cuda) v.s. 448.711444 ms (emu) */
        if (minlogpp > LOGZERO)
            fprintf(stderr, "sMBRerrorsignal: pruning not supported (we won't need it!) :)\n");
        parallelsMBRerrorsignal(parallelstate, thisedgealignments, logpps, amf, logEframescorrect, logEframescorrecttotal, errorsignal, errorsignalneg);
        return;
    }

    //  linear mode
    foreach_coord (i, j, errorsignal)
        errorsignal(i, j) = 0.0f; // Note: we don't actually put anything into the numgammas
    foreach_index (j, edges)
    {
        const auto &e = edges[j];
        if (nodes[e.S].t == nodes[e.E].t) // this happens for dummy !NULL edge at end of file
            continue;
        if (minlogpp > LOGZERO && origlogpps[j] < minlogpp) // this is pruned
            continue;

        size_t ts = nodes[e.S].t;
        size_t te = nodes[e.E].t;

        const double diff = logEframescorrect[j] - logEframescorrecttotal;
        // Note: the contribution of the states of an edge to their senones is the same for all states
        // so we compute it once and add it to all; this will not be the case without hard alignments.
        const double pp = exp(logpps[j]); // edge posterior
        const float edgecorrect = (float) (pp * diff) / amf;
        for (size_t t = ts; t < te; t++)
        {
            const size_t s = thisedgealignments[j][t - ts];
            errorsignal(s, t) += edgecorrect;
        }
    }
}

// compute the error signal for MMI mode
void lattice::mmierrorsignal(parallelstate &parallelstate, double minlogpp, const std::vector<double> &origlogpps,
                             std::vector<msra::math::ssematrixbase *> &abcs, const bool softalignstates,
                             const std::vector<double> &logpps, const msra::asr::simplesenonehmm &hset,
                             const edgealignments &thisedgealignments, msra::math::ssematrixbase &errorsignal) const
{
    if (parallelstate.enabled())
    {
        if (minlogpp > LOGZERO)
            fprintf(stderr, "mmierrorsignal: pruning not supported (we won't need it!) :)\n");
        if (softalignstates)
            LogicError("mmierrorsignal: parallel version for softalignstates mode is not supported yet");
        parallelmmierrorsignal(parallelstate, thisedgealignments, logpps, errorsignal);
        return;
    }

    for (size_t j = 0; j < (errorsignal).cols(); j++)
        for (size_t i = 0; i < (errorsignal).rows(); i++)
            errorsignal(i, j) = VIRGINLOGZERO; // set to zero  --note: may be in-place with logLLs, which now get overwritten

    // size_t warnings = 0;   // [v-hansu] check code for mmi; search this comment to see all related codes
    foreach_index (j, edges)
    {
        const auto &e = edges[j];
        if (nodes[e.S].t == nodes[e.E].t) // this happens for dummy !NULL edge at end of file
            continue;
        if (minlogpp > LOGZERO && origlogpps[j] < minlogpp) // this is pruned
            continue;

        const auto &aligntokens = getaligninfo(j); // get alignment tokens
        auto &loggammas = *abcs[j];

        const float edgelogP = (float) logpps[j];
        // if (islogzero (edgelogP))               // we had a 0 prob
        //    continue;

        // accumulate this edge's gamma matrix into target posteriors
        const size_t tedge = nodes[e.S].t;
        size_t ts = 0;                 // time index into gamma matrix
        size_t js = 0;                 // state index into gamma matrix
        foreach_index (k, aligntokens) // we exploit that units have fixed boundaries
        {
            const auto &unit = aligntokens[k];
            const size_t te = ts + unit.frames;
            const auto &hmm = hset.gethmm(unit.unit); // TODO: inline these expressions
            const size_t n = hmm.getnumstates();
            const size_t je = js + n;
            // P(s) = P(s|e) * P(e)
            for (size_t t = ts; t < te; t++)
            {
                const size_t tutt = t + tedge; // time index w.r.t. utterance
                // double logsum = LOGZERO;         // [v-hansu] check code for mmi; search this comment to see all related codes
                for (size_t i = 0; i < n; i++)
                {
                    const size_t j2 = js + i;             // state index for this unit in matrix
                    const size_t s = hmm.getsenoneid(i); // state class index
                    const float gammajt = loggammas(j2, t);
                    const float statelogP = edgelogP + gammajt;
                    logadd(errorsignal(s, tutt), statelogP);
                }
            }
            ts = te;
            js = je;
        }
        assert(ts + 2 == loggammas.cols() && js == loggammas.rows());
    }

    // check normalizedness (is that an actual English word?)
    // also count non-zero probs
    size_t nonzerostates = 0;
    foreach_column (t, errorsignal)
    {
        double logsum = LOGZERO;
        foreach_row (s, errorsignal)
        {
            if (islogzero(errorsignal(s, t)))
                nonzerostates++;
            else
                logadd(logsum, (double) errorsignal(s, t));
            // TODO: count VIRGINLOGZERO, print per frame
        }
        if (fabs(logsum) / errorsignal.rows() > 1e-6)
            fprintf(stderr, "forwardbackward: WARNING: overall posterior column(%d) sum = exp (%.10f) != 1\n", (int) t, logsum);
    }
    fprintf(stderr, "forwardbackward: %.3f%% non-zero state posteriors\n", 100.0f - nonzerostates * 100.0f / errorsignal.rows() / errorsignal.cols());

    // convert to non-log posterior  --that's what we return
    foreach_coord (i, j, errorsignal)
        errorsignal(i, j) = expf(errorsignal(i, j));
}

// compute ground truth's score
// It is critical to get all details consistent with the lattice, to avoid skewing the weights.
// ... TODO: we don't need this to be a class member, actually; try to just make it a 'static' function.
/*static*/ double lattice::scoregroundtruth(const_array_ref<size_t> uids, const_array_ref<htkmlfwordsequence::word> transcript, const std::vector<float> &transcriptunigrams,
                                            const msra::math::ssematrixbase &logLLs, const msra::asr::simplesenonehmm &hset, const float lmf, const float wp, const float amf)
{
    if (transcript[0].firstframe != 0) // TODO: should we store the #frames instead? Then we can validate the total duration
        LogicError("scoregroundtruth: first transcript[] token does not start at frame 0");

    // get the silence models, since they are treated specially
    const size_t numframes = logLLs.cols();
    const auto &sil = hset.gethmm(hset.gethmmid("sil"));
    const auto &sp = hset.gethmm(hset.gethmmid("sp"));
    if (sp.numstates != 1 || sil.numstates != 3)
        RuntimeError("scoregroundtruth: only supports 1-state /sp/ and 3-state /sil/ tied to /sp/");
    const size_t silst = sp.senoneids[0];

    // loop over words
    double pathscore = 0.0;
    foreach_index (i, transcript)
    {
        size_t ts = transcript[i].firstframe;
        size_t te = ((size_t) i + 1 < transcript.size()) ? transcript[i + 1].firstframe : numframes;
        if (ts >= te)
            LogicError("scoregroundtruth: transcript[] tokens out of order");
        // acoustic score: loop over frames
        const msra::asr::simplesenonehmm::transP *prevtransP = NULL; // previous transP
        int prevs = -1;                                              // previous state index
        for (size_t t = ts; t < te; t++)
        {
            size_t senoneid = uids[t];
            // recover the transP and state index
            int s;
            const msra::asr::simplesenonehmm::transP *transP;
            if (senoneid == silst)
            {
                if (prevtransP == &sil.gettransP()) // "silst" may be tied to /sp/, and thus the /sil/ center state may be ambiguous
                {
                    transP = prevtransP; // remain in /sil/
                    s = 1;               // note that this will fail if /sp/ can follow /sil/ and if /sil/ may end with "silst" (both not allowed currently)
                }
                else // "silst" -> we are in /sp/
                {
                    transP = &sp.gettransP();
                    s = 0;
                }
            }
            else // all others must be non-ambiguous
            {
                int transPindex = hset.senonetransP(senoneid);
                int sindex = hset.senonestate(senoneid);
                if (transPindex == -1 || sindex == -1)
                    RuntimeError("scoregroundtruth: failed to resolve ambiguous senone %s", hset.getsenonename(senoneid));
                transP = &hset.transPs[transPindex];
                s = sindex;
            }
            // if changing phoneme then add necessary enter/exit score
            if (transP != prevtransP)
            {
                if (prevtransP) // previous exit transition
                    pathscore += (*prevtransP)(prevs, prevtransP->getnumstates());
                prevs = -1; // enter transition
            }
            // add inter-state transP
            pathscore += (*transP)(prevs, s);
            // acoustic score
            pathscore += logLLs(senoneid, t);
            // remember state for next frame
            prevtransP = transP;
            prevs = s;
        }
        // add the last exit transition
        if (prevtransP) // previous exit transition
            pathscore += (*prevtransP)(prevs, prevtransP->getnumstates());
        // need to add a /sp/ tee transition if we don't end in sp, since our transcript dictionary includes a /sp/ at the end of each entry
        if (uids[te - 1] != sp.senoneids[0])
            pathscore += sp.gettransP()(-1, sp.numstates);
        // lm
        pathscore += transcriptunigrams[i] * lmf + wp;
    }
    if (islogzero(pathscore))
        fprintf(stderr, "scoregroundtruth: ground-truth path has zero probability; some model inconsistency, maybe?\n");
    // account for amf
    pathscore /= amf;
    fprintf(stderr, "scoregroundtruth: ground-truth score %.6f (%d frames)\n", pathscore, (int) numframes);
    return pathscore;
}

// ---------------------------------------------------------------------------
// sMBRdiagnostics() -- helper to print some diagnostics for analyzing sMBR results
// ---------------------------------------------------------------------------

// static  // with 'static', compiler will complain if the function is not used (we only compile it in sometimes for diagnostics)
void sMBRdiagnostics(const msra::math::ssematrixbase &errorsignal, const_array_ref<size_t> uids,
                     const_array_ref<size_t> bestpath, const vector<bool> &refseen, const msra::asr::simplesenonehmm &hset)
{
    // TODO:
    //  - print best positive runner-up
    //  - WARN tag if neg > pos
    //  - check the sum and warn if not 0
    //  - indicate whether best path state is correct or not
    size_t numcor = 0;
    size_t numnegbetter = 0; // # frames the neg competitor is better
    size_t numposbetter = 0; // # frames the pos competitor is better
    foreach_column (t, errorsignal)
    {
        const size_t sref = uids[t];
        const char *srefname = hset.getsenonename(sref);
        const size_t sbest = bestpath[t];
        const char *sbestname = hset.getsenonename(sbest);
        if (sref == sbest)
            numcor++;
        // for each frame, print error signal for ground truth and runner up (second largest abs value)
        size_t sneg = SIZE_MAX; // competitor
        float eneg = 0.0f;
        size_t spos = SIZE_MAX; // best positive competitor
        float epos = 0.0f;
        foreach_row (s, errorsignal)
        {
            if (s == sref)
                continue;
            if (errorsignal(s, t) < eneg)
            {
                sneg = s;
                eneg = errorsignal(s, t);
            }
            if (errorsignal(s, t) > epos)
            {
                spos = s;
                epos = errorsignal(s, t);
            }
        }
        if (fabs(errorsignal(sref, t)) > 0.0001f && errorsignal(sref, t) < -eneg)
            numnegbetter++;
        if (fabs(errorsignal(sref, t)) > 0.0001f && errorsignal(sref, t) < epos)
            numposbetter++;
        const char *snegname = sneg == SIZE_MAX ? "-" : hset.getsenonename(sneg);
        const char *sposname = spos == SIZE_MAX ? "-" : hset.getsenonename(spos);
        fprintf(stderr, "e(%d): ref %s: %.6f / %s: %.6f / %s: %.6f / top %s: %.6f%s%s%s%s%s\n",
                (int) t, srefname, errorsignal(sref, t), snegname, eneg, sposname, epos, sbestname, errorsignal(sbest, t),
                sbest == sref ? "" : " ERR",
                fabs(errorsignal(sref, t)) > 0.0001f && errorsignal(sref, t) < 0 ? " INV!!" : "",
                fabs(errorsignal(sref, t)) > 0.0001f && errorsignal(sref, t) < -eneg ? " WEAK" : "",
                fabs(errorsignal(sref, t)) > 0.0001f && errorsignal(sref, t) < epos ? " 2ND" : "",
                refseen[t] ? "" : " NOREF");
    }
    // print this to validate our bestpath computation
    fprintf(stderr, "sMBRdiagnostics: %d frames correct out of %d, %.2f%%, neg better in %d, pos in %d\n",
            (int) numcor, (int) errorsignal.cols(), 100.0f * numcor / errorsignal.cols(),
            (int) numnegbetter, (int) numposbetter);
}

// static  // with 'static', compiler will complain if the function is not used (we only compile it in sometimes for diagnostics)
void sMBRsuppressweirdstuff(msra::math::ssematrixbase &errorsignal, const_array_ref<size_t> uids)
{
    size_t numweird = 0;
    foreach_column (t, errorsignal)
    {
        const size_t sref = uids[t];
        // for each frame, print error signal for ground truth and runner up (second largest abs value)
        const float eref = errorsignal(sref, t);
        bool isweird = eref < 0.0f; // negative for reference!?
        for (size_t s = 0; s < errorsignal.rows() && !isweird; s++)
        {
            if (s == sref)
                continue;
            if (fabs(errorsignal(s, t)) > eref)
                isweird = true;
        }
        if (isweird)
        {
            foreach_row (s, errorsignal)
                errorsignal(s, t) = 0.0f;
            numweird++;
        }
    }
    // print this to validate our bestpath computation
    fprintf(stderr, "sMBRsuppressweirdstuff: %d weird frames out of %d, %.2f%% were flattened\n",
            (int) numweird, (int) errorsignal.cols(), 100.0f * numweird / errorsignal.cols());
}

// ---------------------------------------------------------------------------
// forwardbackward() -- main function for MMI/sMBR
//
// This computes the lattice state-level statistics for sequence training using MMI or sMBR.
//
// Outputs, MMI mode:
//  - result = dengammas = denominator gammas (non-log form)
//  - returns log of sum over all paths' likelihoods (the denominator of the MMI objective)
//  - note: numgammas is not used/touched in MMI mode
//
// Outputs, sMBR mode:
// TODO: fix this comment
//  - result = errorsignal = (abs value of) negative contributions to error signal
//  - errorsignalbuf = for temporarily use to get errorsignal
//  - returns expected frames-correct count (the sMBR objective)
// ---------------------------------------------------------------------------

double lattice::forwardbackward(parallelstate &parallelstate, const msra::math::ssematrixbase &logLLs, const msra::asr::simplesenonehmm &hset,
                                msra::math::ssematrixbase &result, msra::math::ssematrixbase &errorsignalbuf,
                                const float lmf, const float wp, const float amf, const float boostingfactor,
                                const bool sMBRmode, array_ref<size_t> uids, const_array_ref<size_t> bounds,
                                const_array_ref<htkmlfwordsequence::word> transcript, const std::vector<float> &transcriptunigrams) const
{
    bool softalign = true;
    bool softalignstates = false;      // true if soft alignment within edges, currently we only support soft within edge in cpu mode
    bool softalignlattice = softalign; // w.r.t. whole lattice

    edgealignments thisedgealignments(*this);   // alignments memory allocate for this lattice
    backpointers thisbackpointers(*this, hset); // memory for forwardbackward

    if (info.numframes != logLLs.cols())
        LogicError("forwardbackward: #frames mismatch between lattice (%d) and LLs (%d)", (int) info.numframes, (int) logLLs.cols());
    // TODO: the following checks should throw, but I don't dare in case this will crash a critical job... if we never see this warning, then
    if (info.numframes != uids.size())
        fprintf(stderr, "forwardbackward: #frames mismatch between lattice (%d) and uids (%d)\n", (int) info.numframes, (int) uids.size());
    if (info.numframes != result.cols())
        fprintf(stderr, "forwardbackward: #frames mismatch between lattice (%d) and result (%d)\n", (int) info.numframes, (int) result.cols());

    littlematrixheap matrixheap(info.numedges); // for abcs

    // PHASE 0: fake word level forward backwards --only used when pruning enabled
    const double minlogpp = LOGZERO; // pruning threshold  --LOGZERO means disabled
    std::vector<double> origlogpps;  // word posterior from original lattice, for pruning decision

    // PHASE 1: per-edge forward backwards (="time alignments")

    // score the ground truth  --only if a transcript is provided, which happens if the user provides a language model
    // TODO: no longer used, remove this. 'transcript' parameter is no longer used in this function.
    transcript;
    transcriptunigrams;

    // allocate alpha/beta/gamma matrices (all are sharing the same memory in-place)
    std::vector<msra::math::ssematrixbase *> abcs;
    std::vector<float> edgeacscores; // [edge index] acoustic scores
    // funcation call for forwardbackward on edge level
    forwardbackwardalign(parallelstate, hset, softalignstates, minlogpp, origlogpps, abcs, matrixheap, sMBRmode /*returnsenoneids*/, edgeacscores, logLLs, thisedgealignments, thisbackpointers, uids, bounds);

// PHASE 2: lattice-level forward backward

// we exploit that the lattice is sorted by (end node, start node) for in-place processing

// checklattice();      // comment out by v-hansu to save time
#ifdef PRINT_TIME_MEASUREMENT
    auto_timer latlevelfwbw;
#endif
    std::vector<double> logpps;
    std::vector<double> Eframescorrectbuf; // this is used for compute the Eframescorrectdiff
    std::vector<double> logEframescorrect; // this is the final output of PHASE 2
    std::vector<double> logalphas;
    std::vector<double> logbetas;
    double totalfwscore = 0; // TODO: name no longer precise in sMBRmode
    double logEframescorrecttotal = LOGZERO;

    bool returnEframescorrect = sMBRmode;
    if (softalignlattice)
    {
        totalfwscore = forwardbackwardlattice(edgeacscores, parallelstate, logpps, logalphas, logbetas, lmf, wp, amf, boostingfactor, returnEframescorrect, (const_array_ref<size_t> &) uids, thisedgealignments, logEframescorrect, Eframescorrectbuf, logEframescorrecttotal);
        if (sMBRmode && !returnEframescorrect)
            logEframescorrecttotal = forwardbackwardlatticesMBR(edgeacscores, hset, logalphas, logbetas, lmf, wp, amf, (const_array_ref<size_t> &) uids, thisedgealignments, Eframescorrectbuf);
        // ^^ BUGBUG not tested
    }
    else
        totalfwscore = bestpathlattice(edgeacscores, logpps, lmf, wp, amf);
#ifdef PRINT_TIME_MEASUREMENT
    latlevelfwbw.show("latlevelfwbw"); // 68.395682 ms
#endif
    if (islogzero(totalfwscore))
    {
        fprintf(stderr, "forwardbackward: WARNING: no path found in lattice (%d nodes/%d edges)\n", (int) nodes.size(), (int) edges.size());
        return LOGZERO; // failed, do not use resulting matrix
    }

    // PHASE 3: compute final state-level posteriors (MMI mode)

    // compute expected frames correct in sMBRmode

    const size_t numframes = logLLs.cols();
    assert(numframes == info.numframes);
    // fprintf (stderr, "forwardbackward: total forward score %.6f (%d frames)\n", totalfwscore, (int) numframes);   // for now--while we are debugging the GPU port

    // MMI mode
    if (!sMBRmode)
    {
        // we first take the sum in log domain to avoid numerical issues
        auto &dengammas = result; // result is denominator gammas
        mmierrorsignal(parallelstate, minlogpp, origlogpps, abcs, softalignstates, logpps, hset, thisedgealignments, dengammas);
        return totalfwscore / numframes; // return value is av. posterior
    }
    // sMBR mode
    else
    {
        auto &errorsignal = result;
        sMBRerrorsignal(parallelstate, errorsignal, errorsignalbuf, logpps, amf, minlogpp, origlogpps, logEframescorrect, logEframescorrecttotal, thisedgealignments);

        static bool dummyvariable = (fprintf(stderr, "note: new version with kappa adjustment, kappa = %.2f\n", 1 / amf), true); // we only print once
        return exp(logEframescorrecttotal) / numframes;                                                                          // return value is av. expected frame-correct count
    }
}
};
};