Estimation.cc

/*
 * $Id: Estimation.cc 1691 2011-08-03 13:38:08Z hahn $
 *
 * Copyright (c) 2004-2005  RWTH Aachen University
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License Version 2 (June
 * 1991) as published by the Free Software Foundation.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, you will find it at
 * http://www.gnu.org/licenses/gpl.html, or write to the Free Software
 * Foundation, Inc., 51 Franlin Street, Fifth Floor, Boston, MA 02110,
 * USA.
 *
 * Should a provision of no. 9 and 10 of the GNU General Public License
 * be invalid or become invalid, a valid provision is deemed to have been
 * agreed upon which comes closest to what the parties intended
 * commercially. In any case guarantee/warranty shall be limited to gross
 * negligent actions or intended actions or fraudulent concealment.
 */

#include "Multigram.hh"
#include "MultigramGraph.hh"
#include "Probability.hh"
#include "Python.hh"
#include "SequenceModel.hh"
#include <tr1/unordered_map>
#include <vector>


class SequenceModelEstimator;

class EstimationGraph : public MultigramGraph {
    friend class EstimationGraphBuilder;
    friend class Accumulator;
    friend class ViterbiAccumulator;
    friend class OneForAllAccumulator;

    NodeList nodesInTopologicalOrder_;
    NodeMap<SequenceModel::History> histories_;
public:
    EstimationGraph() :
	MultigramGraph(),
	histories_(&graph_)
	{}

    void yield() {
	graph_.yield();
	token_.yield();
	histories_.yield();
	probability_.yield();
	NodeList tmp(nodesInTopologicalOrder_);
	nodesInTopologicalOrder_.swap(tmp);
    }

    void updateHistories(const SequenceModel*);
    void updateProbabilities(const SequenceModel*);
#ifdef OBSOLETE
    void draw(FILE*, const StringInventory*, const SequenceModel*) const;
#endif // OBSOLETE

    size_t memoryUsed() const {
	return MultigramGraph::memoryUsed()
	    + nodesInTopologicalOrder_.capacity() * sizeof(NodeList::value_type)
	    + histories_.memoryUsed();
    }
};

#ifdef OBSOLETE
void EstimationGraph::draw(FILE *f, const StringInventory *si, const SequenceModel *sm) const {
    fprintf(f,
	    "digraph \"estimation graph\" {\n"
	    "ranksep = 1.0;\n"
	    "rankdir = LR;\n"
	    "size = \"8.5,11\";\n"
	    "center = 1;\n"
	    "orientation = Landscape\n"
	    "node [fontname=\"Helvetica\"]\n"
	    "edge [fontname=\"Helvetica\"]\n");

    for (Graph::NodeId n = graph_.nodesBegin(); n != graph_.nodesEnd(); ++n) {
	std::string label = (sm)? sm->formatHistory(histories_[n], si) : std::string("?");
	fprintf(f, "n%d [label=\"%s\"] \n", n, label.c_str());
    }

    for (Graph::EdgeId e = graph_.edgesBegin(); e != graph_.edgesEnd(); ++e) {
	std::string label = (si) ? si->symbol(token_[e]) : std::string("?");
	fprintf(f, "n%d -> n%d [label=\"%s %f\"]\n",
		graph_.source(e), graph_.target(e),
		label.c_str(), probability_[e].probability());
    }

    fprintf(f, "}\n");
    fflush(f);
}
#endif // OBSOLETE

void EstimationGraph::updateHistories(const SequenceModel *sm) {
    histories_.sync();
    histories_.fill(0);
    histories_[initial_] = sm->initial();
    histories_[final_] = sm->culDeSac();
    for (NodeList::const_iterator n = nodesInTopologicalOrder_.begin(); n != nodesInTopologicalOrder_.end(); ++n) {
	SequenceModel::History oldHistory = histories_[*n];
	for (Graph::OutgoingEdgeIterator e = graph_.outgoingEdges(*n); e; ++e) {
	    Graph::NodeId target = graph_.target(*e);
	    if (target == final_) continue;
	    SequenceModel::History newHistory = sm->advanced(oldHistory, token_[*e]);
	    verify(!histories_[target] || histories_[target] == newHistory);
	    histories_[target] = newHistory;
	}
    }
}

void EstimationGraph::updateProbabilities(const SequenceModel *sm) {
    probability_.sync();
    for (Graph::EdgeId e = graph_.edgesBegin(); e != graph_.edgesEnd(); ++e) {
	probability_[e] = sm->probability(token_[e], histories_[graph_.source(e)]);
    }
}

// ===========================================================================
class EvidenceStore {
public:
    struct Event {
	SequenceModel::History history;
	SequenceModel::Token token;

	bool operator==(const Event &other) const {
	    return (token   == other.token)
		&& (history == other.history);
	}

	struct Hash {
	    size_t operator() (const Event &e) const {
		size_t h = size_t(e.history);
		h = (h << 4) ^ size_t(e.token);
		return h;
	    }
	};
    };

private:
    typedef std::unordered_map<Event, Probability, Event::Hash> Store;
    Store evidence_;

    const SequenceModel *sequenceModel_;

public:
    EvidenceStore() : sequenceModel_(0) {}

    void setSequenceModel(SequenceModel *sm) {
	sequenceModel_ = sm;
    }

    void accumulate(
	SequenceModel::History history,
	SequenceModel::Token token,
	LogProbability evidence)
    {
	require(token != MultigramInventory::voidIndex()); // don't accumulate with anonymized unknowns!
	Event ev;
	ev.history = history;
	ev.token   = token;
	evidence_[ev] += evidence;
    }

    PyObject *asList() const {
	PyObject *result = PyList_New(evidence_.size());
	u32 i = 0;
	for (Store::const_iterator ev = evidence_.begin(); ev != evidence_.end(); ++ev) {
	    PyObject *hte = Py_BuildValue("(Nif)",
		sequenceModel_->historyAsTuple(ev->first.history),
		int(ev->first.token),
		ev->second.probability());
	    PyList_SET_ITEM(result, i++, hte);
	}
	return result;
    }

    size_t size() const {
	return evidence_.size();
    }

    u32 maximumHistoryLength() const {
	u32 result = 0;
	for (Store::const_iterator ev = evidence_.begin(); ev != evidence_.end(); ++ev)
	    result = std::max(result, sequenceModel_->historyLength(ev->first.history));
	return result;
    }

    Probability maximum() const {
	Probability result = Probability::impossible();
	for (Store::const_iterator ev = evidence_.begin(); ev != evidence_.end(); ++ev)
	    result = std::max(result, ev->second);
	return result;
    }
    Probability total() const {
	Probability result = Probability::impossible();
	for (Store::const_iterator ev = evidence_.begin(); ev != evidence_.end(); ++ev)
	    result += ev->second;
	return result;
    }

    SequenceModelEstimator *makeSequenceModelEstimator() const;

    size_t memoryUsed() const {
#if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 3)
	typedef std::__detail::_Hash_node<Store::value_type, false> StoreNode;
#elif __GNUC__ == 4 && __GNUC_MINOR__ == 2
	typedef std::tr1::__detail::_Hash_node<Store::value_type, false> StoreNode;
#elif __GNUC__ == 4 && __GNUC_MINOR__ <= 1
	typedef Internal::hash_node<Store::value_type, false> StoreNode;
#endif
	return sizeof(EvidenceStore)
	    + evidence_.size() * sizeof(StoreNode)
	    + evidence_.bucket_count() * sizeof(StoreNode*);
    };
};

class Accumulator {
    EvidenceStore *target_;

    ProbabilityAccumulator accu_;
    NodeMap<LogProbability> forw_, bckw_;

    void forward(EstimationGraph *eg) {
	forw_[eg->initial_] = LogProbability::certain();
	for (EstimationGraph::NodeList::const_iterator n = eg->nodesInTopologicalOrder_.begin()+1; n != eg->nodesInTopologicalOrder_.end(); ++n) {
	    verify_(*n != eg->initial_);
	    accu_.clear() ;
	    for (Graph::IncomingEdgeIterator e = eg->graph_.incomingEdges(*n); e; ++e)
		accu_.add(forw_[eg->graph_.source(*e)] * eg->probability_[*e]);
	    forw_[*n] = accu_.sum();
	}
    }

    void backward(EstimationGraph *eg) {
	bckw_[eg->final_] = LogProbability::certain();
	for (EstimationGraph::NodeList::reverse_iterator n = eg->nodesInTopologicalOrder_.rbegin()+1; n != eg->nodesInTopologicalOrder_.rend(); ++n) {
	    verify_(*n != eg->final_);
	    accu_.clear() ;
	    for (Graph::OutgoingEdgeIterator e = eg->graph_.outgoingEdges(*n); e; ++e)
		accu_.add(bckw_[eg->graph_.target(*e)] * eg->probability_[*e]);
	    bckw_[*n] = accu_.sum();
	}
    }

public:
    Accumulator() : target_(0) {}

    void setTarget(EvidenceStore *es) {
	target_ = es;
    }

    LogProbability accumulate(EstimationGraph *eg, LogProbability weight) {
	forw_.sync(&eg->graph_); forward(eg);
	bckw_.sync(&eg->graph_); backward(eg);
#if 1 // DEBUG
	if (!isNearlyEqual(forw_[eg->final_].score(), bckw_[eg->initial_].score(), 100)) {
		std::cerr << __FILE__ << ":" << __LINE__ << "\t"
			  << forw_[eg->final_].score() << "\t"
			  << bckw_[eg->initial_].score() << std::endl;
	}
#endif
	LogProbability total = (forw_[eg->final_] * bckw_[eg->initial_]).pow(0.5);
	for (Graph::EdgeId e = eg->graph_.edgesBegin(); e != eg->graph_.edgesEnd(); ++e) {
	    Graph::NodeId source = eg->graph_.source(e);
	    LogProbability post
		= forw_[source]
		* eg->probability_[e]
		* bckw_[eg->graph_.target(e)]
		/ total;
#if 1 // DEBUG
	    if (post > LogProbability::certain() &&
		!isNearlyEqual(post.probability(), 1.0, 100)) {
		std::cerr << __FILE__ << ":" << __LINE__ << "\t"
			  << forw_[eg->final_].score() << "\t"
			  << bckw_[eg->initial_].score() << "\t"
			  << total.score() << "\t"
			  << forw_[source].score() << "\t"
			  << eg->probability_[e].score() << "\t"
			  << bckw_[eg->graph_.target(e)].score() << "\t"
			  << post.score() << std::endl;
	    }
#endif
	    target_->accumulate(eg->histories_[source], eg->token_[e], weight * post);
	}
	return total;
    }

    LogProbability logLik(EstimationGraph *eg) {
	forw_.sync(&eg->graph_); forward(eg);
#if 1
	return forw_[eg->final_];
#else
	bckw_.sync(&eg->graph_); backward(eg);
	return (forw_[eg->final_] * bckw_[eg->initial_]).pow(0.5);
#endif
    }
};

class ViterbiAccumulator {
    EvidenceStore *target_;

    NodeMap<LogProbability> forw_;
    NodeMap<Graph::EdgeId> back_;

    void forward(EstimationGraph *eg) {
	forw_[eg->initial_] = LogProbability::certain();
	for (EstimationGraph::NodeList::const_iterator n = eg->nodesInTopologicalOrder_.begin()+1; n != eg->nodesInTopologicalOrder_.end(); ++n) {
	    verify_(*n != eg->initial_);
	    LogProbability bestProb = LogProbability::impossible();
	    Graph::EdgeId bestBack = 0;
	    for (Graph::IncomingEdgeIterator e = eg->graph_.incomingEdges(*n); e; ++e) {
		LogProbability f = forw_[eg->graph_.source(*e)] * eg->probability_[*e];
		if (bestProb < f) {
		    bestProb = f;
		    bestBack = *e;
		}
	    }
	    forw_[*n] = bestProb;
	    back_[*n] = bestBack;
	}
    }

public:
    ViterbiAccumulator() : target_(0) {}

    void setTarget(EvidenceStore *es) {
	target_ = es;
    }

    LogProbability accumulate(EstimationGraph *eg, LogProbability weight) {
	forw_.sync(&eg->graph_); back_.sync(&eg->graph_); forward(eg);

	for (Graph::NodeId n = eg->final_; n != eg->initial_;) {
	    Graph::EdgeId e = back_[n];
	    Graph::NodeId source = eg->graph_.source(e);
	    target_->accumulate(eg->histories_[source], eg->token_[e], weight);
	    n = source;
	}
	return forw_[eg->final_];
    }

    LogProbability logLik(EstimationGraph *eg) {
	forw_.sync(&eg->graph_); back_.sync(&eg->graph_); forward(eg);
	return forw_[eg->final_];
    }

    LogProbability segment(EstimationGraph *eg, std::vector<MultigramIndex> &result) {
	forw_.sync(&eg->graph_); back_.sync(&eg->graph_); forward(eg);
	result.clear();
	for (Graph::NodeId n = eg->final_; n != eg->initial_;) {
	    Graph::EdgeId e = back_[n];
	    result.push_back(eg->token_[e]);
	    n = eg->graph_.source(e);
	}
	std::reverse(result.begin(), result.end());
	return forw_[eg->final_];
    }
};

class OneForAllAccumulator {
    EvidenceStore *target_;

public:
    OneForAllAccumulator() : target_(0) {}

    void setTarget(EvidenceStore *es) {
	target_ = es;
    }

    void accumulate(EstimationGraph *eg, LogProbability weight) {
	for (Graph::EdgeId e = eg->graph_.edgesBegin(); e != eg->graph_.edgesEnd(); ++e) {
	    target_->accumulate(eg->histories_[eg->graph_.source(e)], eg->token_[e], weight);
	}
    }
};

// ===========================================================================
class EstimationGraphBuilder :
    public GraphSorter
{
public:
    enum MultigramEmergenceMode {
	emergeNewMultigrams,
	suppressNewMultigrams,
	anonymizeNewMultigrams
    };

private:
    struct SizeTemplate {
	u32 left, right;
    };
    typedef std::vector<SizeTemplate> SizeTemplateList;
    SizeTemplateList sizeTemplates_;
    MultigramEmergenceMode multigramEmergence_;

public:
    void clearSizeTemplates() {
	sizeTemplates_.clear();
    }
    void addSizeTemplate(int left, int right) {
	require(left >= 0);
	require(right >= 0);
	require(left > 0 || right > 0);
	require(left  <= Multigram::maximumLength);
	require(right <= Multigram::maximumLength);
	SizeTemplate st;
	st.left = left;
	st.right = right;
	sizeTemplates_.push_back(st);
    }

    void setEmergenceMode(MultigramEmergenceMode mode) {
	multigramEmergence_ = mode;
    }

private:
    MultigramInventory *inventory_;
    SequenceModel *sequenceModel_;
public:
    void setSequenceModel(MultigramInventory *mi, SequenceModel *sm) {
	inventory_ = mi;
	sequenceModel_ = sm;
    }

private:
    Sequence left_, right_;
    EstimationGraph *target_;

    struct NodeDesc {
	struct {
	    u32 left, right;
	} position;
	SequenceModel::History history;

	bool operator==(const NodeDesc &other) const {
	    return (position.left  == other.position.left)
		&& (position.right == other.position.right)
		&& (history        == other.history);
	}

	struct Hash {
	    size_t operator() (const NodeDesc &s) const {
		size_t h = size_t(s.history);
		h = (h << 0) ^ s.position.left;
		h = (h << 4) ^ s.position.right;
		return h;
	    }
	};
    };

    static const Graph::NodeId newNode = 0;
    static const Graph::NodeId greyNode = 0xfffffff;
    static const Graph::NodeId deadNode = 0xffffffe;
    typedef std::unordered_map<NodeDesc, Graph::NodeId, NodeDesc::Hash> NodeStateMap;
    NodeStateMap nodeStates_;
    typedef std::pair<NodeDesc, SizeTemplateList::const_iterator> DfsStackItem;
    typedef std::vector<DfsStackItem> DfsStack;
    DfsStack stack_;

    bool expand(
	const NodeDesc &current,
	SizeTemplateList::const_iterator st,
	NodeDesc &next,
	SequenceModel::Token &token)
    {
	next.position.left = current.position.left + st->left;
	if (next.position.left > left_.size())
	    return false;
	next.position.right = current.position.right + st->right;
	if (next.position.right > right_.size())
	    return false;

	JointMultigram jmg(
	    &left_ [current.position.left ], &left_ [next.position.left ],
	    &right_[current.position.right], &right_[next.position.right]);

	switch (multigramEmergence_) {
	case emergeNewMultigrams:
	    token = inventory_->index(jmg);
	    break;
	case suppressNewMultigrams:
	    token = inventory_->testIndex(jmg);
	    if (token == MultigramInventory::voidIndex())
		return false;
	    break;
	case anonymizeNewMultigrams:
	    token = inventory_->testIndex(jmg);
	    break;
	default: defect();
	}

	next.history = sequenceModel_->advanced(current.history, token);

	return true;
    }

    bool isFinal(const NodeDesc &current) {
	return (current.position.left  == left_ .size() &&
		current.position.right == right_.size());
    }

    void explore() {
	while (!stack_.empty()) {
	    DfsStackItem &top(stack_.back());
	    NodeDesc &current(top.first);
	    SizeTemplateList::const_iterator &st(top.second);
	    Graph::NodeId currentState = nodeStates_[current];
	    verify(currentState != newNode);
	    verify(currentState != deadNode);

	    if (isFinal(current)) {
		verify(nodeStates_[current] == greyNode);
		if (!target_->final_) {
		    target_->final_ = target_->graph_.newNode();
		    target_->nodesInTopologicalOrder_.push_back(target_->final_);
		}
		Graph::NodeId currentState = nodeStates_[current] = target_->graph_.newNode();
		target_->nodesInTopologicalOrder_.push_back(currentState);
		Graph::EdgeId edge = target_->graph_.newEdge(currentState, target_->final_);
		target_->token_.set(edge, sequenceModel_->term());
		stack_.pop_back();
	    } else if (st != sizeTemplates_.end()) {
		NodeDesc next;
		SequenceModel::Token token;
		if (expand(current, st++, next, token)) {
		    Graph::NodeId nextState = nodeStates_[next];
		    if (nextState == newNode) {
			--st;
			nodeStates_[next] = greyNode;
			stack_.push_back(DfsStackItem(next, sizeTemplates_.begin()));
		    } else if (nextState == greyNode) {
			defect(); // cycle detected!
		    } else if (nextState != deadNode) {
			if (currentState == greyNode) {
			    currentState = nodeStates_[current] = target_->graph_.newNode();
			}
			Graph::EdgeId edge = target_->graph_.newEdge(currentState, nextState);
			target_->token_.set(edge, token);
		    }
		}
	    } else {
		if (currentState == greyNode)
		    nodeStates_[current] = deadNode;
		else
		    target_->nodesInTopologicalOrder_.push_back(currentState);
		stack_.pop_back();
	    }
	}
    }

public:
    EstimationGraphBuilder() :
	multigramEmergence_(emergeNewMultigrams),
	inventory_(0),
	sequenceModel_(0),
	target_(0)
    {}

    void build(EstimationGraph *eg, const Sequence &left, const Sequence &right) {
	left_   = left;
	right_  = right;
	target_ = eg;

	target_->graph_.clear();
	target_->initial_ = target_->final_ = 0;
	NodeDesc initial;
	initial.position.left = initial.position.right = 0;
	initial.history = sequenceModel_->initial();
	nodeStates_[initial] = greyNode;
	stack_.push_back(DfsStackItem(initial, sizeTemplates_.begin()));
	target_->nodesInTopologicalOrder_.clear();
	explore();
	target_->initial_ = nodeStates_[initial];
	nodeStates_.clear();
	std::reverse(target_->nodesInTopologicalOrder_.begin(),
		     target_->nodesInTopologicalOrder_.end());

	verify(target_->initial_ != greyNode);
	verify(target_->initial_ != newNode);
	if (target_->initial_ == deadNode)
	    throw std::runtime_error("final node not reachable");

	verify(target_->nodesInTopologicalOrder_.size() == target_->graph_.nNodes() - 1);
	verify(target_->nodesInTopologicalOrder_.front() == target_->initial_);
	verify(target_->nodesInTopologicalOrder_.back() == target_->final_);

	target_->updateHistories(sequenceModel_);
	target_->updateProbabilities(sequenceModel_);
    }

    EstimationGraph *create(const Sequence &left, const Sequence &right) {
	EstimationGraph *result = new EstimationGraph;
	try {
	    build(result, left, right);
	} catch (...) {
	    delete result; result = 0;
	    throw;
	}
	result->yield();
	return result;
    }

    void update(EstimationGraph *eg) {
	eg->updateHistories(sequenceModel_);
	eg->updateProbabilities(sequenceModel_);
    }

    size_t memoryUsed() const {
#if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 3)
	typedef std::__detail::_Hash_node<NodeStateMap::value_type, false> NodeStateMapNode;
#elif __GNUC__ == 4 && __GNUC_MINOR__ == 2
	typedef std::tr1::__detail::_Hash_node<NodeStateMap::value_type, false> NodeStateMapNode;
#elif __GNUC__ == 4 && __GNUC_MINOR__ <= 1
	typedef Internal::hash_node<NodeStateMap::value_type, false> NodeStateMapNode;
#endif
	return sizeof(EstimationGraphBuilder)
	    + GraphSorter::memoryUsed() - sizeof(GraphSorter)
	    + sizeTemplates_.capacity() * sizeof(SizeTemplateList::value_type)
	    + left_.capacity() * sizeof(Sequence::value_type)
	    + right_.capacity() * sizeof(Sequence::value_type)
	    + nodeStates_.size() * sizeof(NodeStateMapNode)
	    + nodeStates_.bucket_count() * sizeof(NodeStateMapNode*)
	    + stack_.capacity() * sizeof(DfsStack::value_type);
    }

};

// ===========================================================================
class SequenceModelEstimator {
private:
    friend class EvidenceStore;

    struct Item : public EvidenceStore::Event {
	Probability evidence;
	Probability probability;

	struct Ordering {
	    bool operator() (const Item &a, const Item &b) const {
		if (a.history != b.history) return a.history < b.history;
		return a.token < b.token;
	    }
	};
	struct TokenOrdering {
	    bool operator() (const Item &a, const Item &b) const {
		return a.token < b.token;
	    }
	};
    };
    typedef std::vector<Item> ItemList;

    struct Group {
	struct { ItemList::iterator begin, end; } items;
	Probability total;
	Probability backOffWeight;
	Group() : total() {}
    };
    typedef std::unordered_map<SequenceModel::History, Group, Core::conversion<SequenceModel::History, size_t> > GroupStore;

    const SequenceModel *sequenceModel_;
    ItemList items;
    Item::Ordering itemOrdering;
    Item::TokenOrdering itemTokenOrdering;
    GroupStore groups;
    std::vector<std::vector<SequenceModel::History> > historiesByLength;

    void init(const SequenceModel*);
    void reset();
    void doKneserNeyDiscounting(const std::vector<double> &discounts);
    void computeProbabilities(double vocabularySize);
public:
    void makeSequenceModel(
	SequenceModel *target,
	double vocabularySize,
	const std::vector<double> &discounts);
};

SequenceModelEstimator *EvidenceStore::makeSequenceModelEstimator() const {
    SequenceModelEstimator *sme = new SequenceModelEstimator();

    sme->items.clear();
    SequenceModelEstimator::Item item;
    std::unordered_map<Event, size_t, Event::Hash> pos;
    for (Store::const_iterator ev = evidence_.begin(); ev != evidence_.end(); ++ev) {
	item.history     = ev->first.history;
	item.token       = ev->first.token;
	item.evidence    = ev->second;
	pos[item] = sme->items.size();
	sme->items.push_back(item);
    }

    // ensure all events we are going to discout to are present
    for (size_t i = 0; i < sme->items.size(); ++i) {
	item = sme->items[i];
	item.evidence = Probability(0.0);
	item.history = sequenceModel_->shortened(item.history);
	if (item.history && pos.find(item) == pos.end()) {
	    pos[item] = sme->items.size();
	    sme->items.push_back(item);
	}
    }

    sme->init(sequenceModel_);
    return sme;
}

void SequenceModelEstimator::makeSequenceModel(
    SequenceModel *target,
    double vocabularySize,
    const std::vector<double> &discounts)
{
    reset();
    doKneserNeyDiscounting(discounts);
    computeProbabilities(vocabularySize);

    std::auto_ptr<SequenceModel::InitData> data(new SequenceModel::InitData);
    std::vector<Token> history;
    for (GroupStore::const_iterator g = groups.begin(); g != groups.end(); ++g) {
	sequenceModel_->historyAsVector(g->first, history);
	std::reverse(history.begin(), history.end());
	data->setHistory(&*history.begin(), &*history.end());

	if (g->second.backOffWeight != Probability(1.0))
	    data->addBackOffWeight(g->second.backOffWeight);
	for (ItemList::const_iterator i = g->second.items.begin; i != g->second.items.end; ++i) {
	    if (i->probability > Probability(0.0))
		data->addProbability(i->token, i->probability);
	}
    }
    target->setInitAndTerm(sequenceModel_->init(), sequenceModel_->term());
    target->set(&*data);
}

void SequenceModelEstimator::init(const SequenceModel *sm) {
    require(items.size() > 0);

    sequenceModel_ = sm;

    std::sort(items.begin(), items.end(), itemOrdering);

    groups.clear();
    historiesByLength.clear();
    Group ng;
    ItemList::iterator &i(ng.items.begin);
    i = items.begin();
    GroupStore::iterator g = groups.insert(std::make_pair(i->history, ng)).first;
    for (++i; i != items.end(); ++i) {
	if (i->history != g->first) {
	    g->second.items.end = i;
	    g = groups.insert(std::make_pair(i->history, ng)).first;
	}
    }
    g->second.items.end = items.end();

    for (GroupStore::const_iterator g = groups.begin(); g != groups.end(); ++g) {
	u32 hl = sequenceModel_->historyLength(g->first);
	if (hl >= historiesByLength.size()) historiesByLength.resize(hl+1);
	historiesByLength[hl].push_back(g->first);
    }
}

void SequenceModelEstimator::reset() {
    for (ItemList::iterator i = items.begin(); i != items.end(); ++i)
	i->probability = i->evidence;
}

void SequenceModelEstimator::doKneserNeyDiscounting(const std::vector<double> &discounts) {
    require(historiesByLength.size() > 0);
    require(discounts.size() >= historiesByLength.size());
    for (u32 level = historiesByLength.size()-1; level > 0; --level) {
	Probability discount(discounts[level]);
	std::vector<SequenceModel::History>::const_iterator h, h_end = historiesByLength[level].end();
	for (h = historiesByLength[level].begin(); h != h_end; ++h) {
	    SequenceModel::History history = *h;
	    SequenceModel::History shorterHistory = sequenceModel_->shortened(history);
	    Group &g(groups[history]);
	    verify_(groups.find(shorterHistory) != groups.end());
	    Group &sg(groups[shorterHistory]);
	    Probability sum;
	    ItemList::iterator si = sg.items.begin;
	    for (ItemList::iterator i = g.items.begin; i != g.items.end; ++i) {
#if 0 // DEBUG
		std::cout << __FILE__ << ":" << __LINE__ << "\t"
			  << i - items.begin() << "\t"
			  << "(" << sequenceModel_->formatHistory(i->history) << ")\t"
			  << i->token << "\t"
			  << i->probability.probability() << std::endl;
#endif
		Item sItem = *i;
		sum += i->probability;
		if (i->probability > discount) {
		    i->probability -= discount;
		    sItem.probability = discount;
		} else {
		    i->probability = Probability(0.0);
		}
		verify_(si != sg.items.end); while (itemTokenOrdering(*si, sItem)) {
		    ++si;
#if 0 // DEBUG
		    std::cout << __FILE__ << ":" << __LINE__ << "\t"
			      << si - items.begin() << "\t"
			      << "(" << sequenceModel_->formatHistory(si->history) << ")\t"
			      << si->token << "\t"
			      << si->probability.probability() << std::endl;
#endif
		    verify_(si != sg.items.end);
		}
		verify(si->token == sItem.token);
		si->probability += sItem.probability;
	    }
	    g.total = sum;
	}
    }
    u32 level = 0;
    Probability discount(discounts[level]);
    std::vector<SequenceModel::History>::const_iterator h, h_end = historiesByLength[level].end();
    for (h = historiesByLength[level].begin(); h != h_end; ++h) {
	SequenceModel::History history = *h;
	Group &g(groups[history]);
	Probability sum;
	for (ItemList::iterator i = g.items.begin; i != g.items.end; ++i) {
	    sum += i->probability;
	    if (i->probability > discount) {
		i->probability -= discount;
	    } else {
		i->probability = Probability(0.0);
	    }
	}
	g.total = sum;
    }
}

void SequenceModelEstimator::computeProbabilities(double vocabularySize) {
    // compute probabilities with interpolation
    Probability zeroGramProbability(1.0 / vocabularySize);

    for (u32 level = 0; level < historiesByLength.size(); ++level) {
	std::vector<SequenceModel::History>::const_iterator h, h_end = historiesByLength[level].end();
	for (h = historiesByLength[level].begin(); h != h_end; ++h) {
	    SequenceModel::History history = *h;
	    Group &g(groups[history]);

	    Probability sum(0.0);
	    for (ItemList::const_iterator i = g.items.begin; i != g.items.end; ++i)
		sum += i->probability;
	    if (sum > g.total) {
		g.backOffWeight = Probability(0.0);
	    } else if (sum <= Probability(0.0)) {
		g.backOffWeight = Probability(1.0);
	    } else {
		g.backOffWeight = (sum / g.total).complement();
	    }

	    SequenceModel::History shorterHistory = sequenceModel_->shortened(history);
	    if (shorterHistory == 0) {
		g.backOffWeight *= zeroGramProbability;
		for (ItemList::iterator i = g.items.begin; i != g.items.end; ++i) {
		    if (i->probability == Probability(0.0)) continue;
		    i->probability = i->probability / g.total + g.backOffWeight;
		}
	    } else {
		const Group &sg(groups[shorterHistory]);
		ItemList::const_iterator si = sg.items.begin;
		for (ItemList::iterator i = g.items.begin; i != g.items.end; ++i) {
		    if (i->probability == Probability(0.0)) continue;
		    Probability pLowerOrder;
		    verify_(si != sg.items.end); while (itemTokenOrdering(*si, *i)) { ++si; verify_(si != sg.items.end); }
		    verify_(si->token == i->token);
		    if (si->probability > Probability(0.0)) {
			pLowerOrder = si->probability;
		    } else {
			pLowerOrder = sg.backOffWeight;
			for (SequenceModel::History boh = sequenceModel_->shortened(shorterHistory); boh; boh = sequenceModel_->shortened(boh)) {
			    Group &bog(groups[boh]);
			    ItemList::const_iterator boi = std::lower_bound(bog.items.begin, bog.items.end, *i, itemTokenOrdering);
			    verify_(boi != bog.items.end && boi->token == i->token);
			    if (boi->probability > Probability(0.0)) {
				pLowerOrder *= boi->probability;
				break;
			    }
			    pLowerOrder *= bog.backOffWeight;
			}
		    }

		    i->probability = i->probability / g.total + g.backOffWeight * pLowerOrder;
		}
	    }
	} // for (h)
    } // for (level)
}