/*****************************************************************************\
 * Filename : pp_profile.hh
 * Author   : Oliver Keller
 * Project  : Gene Prediction with Protein Family Patterns
 *
 * Description: Data structures containing the ProteinProfile /
 *              IntronProfile 
 *
 *
 * Date       |   Author              |  Changes
 *------------|-----------------------|---------------------------------
 * 22.03.07   | Oliver Keller         | creation of the file
 * 06.11.08   | Oliver Keller         | first version in SVN
 * 28.01.09   | Oliver Keller         | rev. 157, final forward-only version
 \******************************************************************************/

#ifndef __PP_PROFILE_HH
#define __PP_PROFILE_HH

// project includes
#include "geneticcode.hh" // needed to score codons
#include "vitmatrix.hh"   // for SubstateId's

// standard C/C++ includes
#include <deque>
#include <cmath>  // for sqrt
#include <map>
#include <vector>

using namespace std;


// constants used in profile classes
// how many bits are available for i when coding PP::Position as int
#define INTERBLOCKBITS 15
#define MAXINTERBLOCKDIST ((1<<INTERBLOCKBITS)-1)
// this might well be reduced in the future to enable more profile-related
// substates
#define MAX_BLOCKCOUNT 64


// fraction by which the distance range is relaxed
// PP::DistanceType::makeTolerant allows for a range of (1-x)min,(1+x)max
#define RELAXATION 0.05

// minimum size of a partial block that can be discarded
#define MIN_CHECKCOUNT 3
#define MIN_BLOCKSIZE 6

// default values for quantiles (in units of stddev)
// that are transformed into scoring thresholds
#define MIN_SPEC 4.0
#define MIN_SENS 0.4
#define MIN_ANCHOR_SPEC 4.0
#define MIN_ANCHOR_COUNT 1
#define PARTIAL_SPEC 4.5
#define PARTIAL_SENS 2.0
#define GLOBAL_THRESH 2.5


struct Dist {
    double mu;
    double var;
    Dist(double m=0, double v=0) : mu(m), var(v) {}
    Dist& add(double m=0, double v=0) { 
	mu+=m; var+=v; return *this;
    }
    Dist& operator+=(const Dist& other) {
	return add(other.mu, other.var); 
    }
    Dist& operator-=(const Dist& other) {
	return add(-other.mu, -other.var); 
    }
    Dist operator+(const Dist& other) const {
	return Dist(*this) += other;
    }
    Dist operator-(const Dist& other) const {
	return Dist(*this) -= other;
    }

    double stddev() const {
	return sqrt(var>=0?var:-var);
    }

    double normed(double abs) const {
	return (abs - mu)/stddev();
    }
    double abs(double normed) const {
	return normed * stddev() + mu;
    }

};

inline Dist operator-(const Dist& d) {
    return Dist(-d.mu, -d.var);
}

/*
 * A data structure representing the intron profile of a block profile
 *
 */
struct ProfileMap : private map < pair<int, int>, double > {
    void add(int n, int f, double p) {
	insert(make_pair(make_pair(n,f), p));
    }
    double get(int n, int f) const {
	const_iterator it = find(make_pair(n,f));
	return (it == end()) ? 0 : it->second;
    }
}; // class ProfileMap


namespace PP {
    /*
     * Exception class thrown from PP:: classes
     */
    struct ProfileNotFoundError : ProjectError {
 	ProfileNotFoundError(string filename) :
	    ProjectError("PP::Profile: Could not open file \"" + filename + "\".")
	{}
    };
    struct ProfileReadError : ProjectError {
	ProfileReadError(string filename, int lineno) :
	    ProjectError("PP::Profile: Error parsing pattern file\"" + filename + 
			 "\", line " + itoa(lineno) + ".")
	{}
    };

    struct ProfileParseError {
	int lineno; 
	ProfileParseError(int n) : lineno(n) {}
    };

    struct PartParseError {
	int offset;
	PartParseError(int n) : offset(n) {}
    };

    struct ProfileInsigError : ProjectError {
	ProfileInsigError(string message) :
	    ProjectError("PP::Profile: " + message) {}
    };
    
    
    /*
     * A data structure representing an integer interval
     */
    struct Range {
	Range(int x=0) : min(x), max(x) {}
	Range(int l, int h) : min(l), max(h) { 
	    if (min>max) { min = h; max = l; }
	}
	Range operator+= (Range other) {
            min+=other.min;
	    max+=other.max;
            return *this;
        }
        Range operator+ (Range other) const {
	    return Range(min,max)+=other;
        }
	bool has(int elem) const {
	    return min <= elem && elem <= max;
	}
	int size() const { // <= 0 is equivalent to
	    return max-min+1;
	}

	// move interval such that max equals newmax
	Range alignRight(int newmax) const {
	    return Range(newmax-max+min, newmax);
	}
	// move interval such that min equals newmin
	Range alignLeft(int newmin) const {
	    return Range(newmin, newmin+max-min);
	}
	
	// data fields
	int min, max;
    }; // class PP::Range

    inline Range operator- (Range r) {
	return Range(-r.max, -r.min);
    }

    struct DistanceType {
	bool has_max;
	Range r;
	DistanceType() : has_max(true), r() {}
	DistanceType(Range other) : has_max(true), r(other) {}
	friend ostream& operator<< (ostream& strm, DistanceType d) {
	    if (d.has_max)
		return strm <<  d.r.min << "\t" << d.r.max;
  	    else
 		return strm << d.r.min << "\t*";
	}
        friend istream& operator>> (istream& strm, DistanceType& d) {
	    d.has_max = true;
	    strm >> d.r.min >> ws;
	    if (strm.peek() == '*') {
		strm.get();
		d.has_max = false;
		d.r.max = d.r.min;
		return strm;
	    }
	    return strm >> d.r.max;
        }
	DistanceType& operator+= (DistanceType other) {
	    r += other.r;
	    has_max = has_max && other.has_max;
	    return *this;
	}
	DistanceType operator+ (DistanceType other) const {
	    return DistanceType(*this) += other;
	}
	void setInfMax () {
	    r.max = r.min;
	    has_max = false;
	}
	void makeTolerant() {
	    r.min = int(r.min * (1-RELAXATION) + 0.5);
	    if (has_max) {
		r.max = int(r.max * (1+RELAXATION) + 0.5);
		if (r.max >= MAXINTERBLOCKDIST) 
		    setInfMax();
	    } else
		r.max = r.min;
	}
	bool has(int elem) const {
	    return has_max ? r.has(elem) : r.min <= elem;
	}
    };

    /*
     * A data structure representing a distribution of amino acids
     */
    class Column {
    public:
	Column() { values[0]=1.0; }
 	Column(const double* val) {
	    *this = val;
	}
	Column& operator= (const double*);

	// initialize odds-ratios factors 
	void initRatios();

	// return expectation and variance of odds-ratios
	Dist getDist(const Column& model) const;
	Dist getOwnDist() const {
	    return getDist(*this);
	}
	Dist getBackDist() const {
	    return getDist(background);
	}

	// accessing the Column entries
	// by value
	double operator[](int n) const { 
	    return (0 <= n && n<NUM_AA) ? values[n] : -1.0; 
	}
	// by oddRatio
	Double Q(int n) const {
	    return (0 <= n && n<NUM_AA) ? oddRatios[n] : stopCodonScore; 
	}
	Double Q(char c) const {
	    return Q(GeneticCode::get_aa_from_symbol(c));
	}
	// by LogOddRatio
	double L(int n) const { 
	    return Q(n).log(); 
	}

	// stream operations
	friend ostream& operator<<(ostream&, const Column&);
	friend istream& operator>>(istream&, Column& c);
	
	static const Double stopCodonScore;
	static double invalidScore;
	static double weight;
    private:
 	// internal data fields
        double values[NUM_AA];  // these numbers will be positive and
	                        // add up to one; class takes care of the normalisation 
	double oddRatios[NUM_AA]; // the odds against the background distribution

	static const Column background;
	static const double minFreq;   // a frequency of 0.0001 is always ensured by the class;
	                               // but we don't deal here with pseudocounts
    }; // class PP::Column



    struct Position {
	Position() :  b(-1), i(0) {}

	Position(int blockno, int columnNr) :
	    b(blockno), i(columnNr) {}
	Position(const SubstateId& id) : 
	    b(id.slot % MAX_BLOCKCOUNT),
	    i(id.value) {}
	Position operator+ (int l) {
	    return Position(b, i + l);
	}
	Position operator- (int l) {
	    return operator+(-l);
	}
	Position operator++(int) {
	    return Position(b,i++);
	}
	bool operator< (Position other) const {
	    return (b < other.b) || (b == other.b && i < other.i);
	}
	bool operator== (Position other) const {
	    return i==other.i && b==other.b;
	}
	bool operator!= (Position other) const {
	    return !(other == (*this));
	}
	bool operator<= (Position other) const {
	    return !(other < (*this));
	}
	void nextB() { b++; i=0; }

	// static int getB(int id) {
	//     return id >> (INTERBLOCKBITS+1);
	// }
	// static int getI(int id) {
	//     static const int mask = (-1) << (INTERBLOCKBITS+1);
	//     static const int offset = 1 << INTERBLOCKBITS;
	//     return (id | mask) + offset;
	// }
	// static Position fromId(int id) {
	//     if (id < 0) return Position();
	//     return Position(getB(id), getI(id));
	// }

	Position operator= (const SubstateId& id) {
	    b = id.slot % MAX_BLOCKCOUNT;
	    i = id.value;
	    return *this;
	}
	SubstateId id() const {
	    return b<0 ? SubstateId() : SubstateId(b,i);
	}
	SubstateId copyId() const {
	    return b<0 ? SubstateId() : SubstateId(b + MAX_BLOCKCOUNT, i);
	}

     	short b, i;
   }; // class PP::Position
    
    class BlockScoreType;
    
    struct PartScoreType {
	int from, to;
	double score;
    };

    class Block {
//	friend class BlockScoreType;
    public:
	Block(DistanceType, const vector<string>&, string);
	~Block() {
	    delete iP;
	}
	// score(x,b,l) is the emission probability factor for the sequence starting 
	// at dna position x, from the partial block of length l (codons), starting at 
	// offset b
	// BlockScoreType score;
	void initDistributions();
  	bool initThresholds();
	void setIntronProfile(const vector<string>&);
	int size() const { return columns.size(); }
	const Column& operator[] (int columnNr) const { return columns[columnNr]; }
	// bool addBlocksUntil(int newbase, map<int,Double>* result=NULL);
	// void removeBlocksUntil(int newbase) {
	//     newbase -= 3 * size();
	//     while (score.begin() < newbase)
	// 	score.pop_front();
	// }
//	Double saveNewScore();
	// Double savedScore(int offset) {
	//     map<int,Double>::iterator it = hits[offset%3].find(offset);
	//     return it == hits[offset%3].end() ? 0 : it->second;
	// }
	Double scoreFromScratch(bool complement, int dna_offset, int block_offset, int len) const;
	Double scoreFromScratch(bool complement, int dna_offset, int block_offset) const {
	    return scoreFromScratch(complement, dna_offset, block_offset, size()-block_offset);
	}
	Double scoreFromScratch(bool complement, int dna_offset) const {
	    return scoreFromScratch(complement, dna_offset, 0);
	}
	Double checkedSuffixScore(bool complement, int dna_offset, int block_offset) const;
	void bestPartialLogScore(bool complement, int dna_offset, PartScoreType&) const;

	// Double checkedPrefixScore(bool complement, int dna_offset, int len) const {
	//     Double result = scoreFromScratch(complement, dna_offset, 0,len);
	//     return result > getPrefixThresh(complement, len) ? result : 0;
	// }
	
	// const map<int,Double>& savedHits(int f) const {
	//     return hits[f];
	// }
#ifdef DEBUG
	void show_score(int n);
#endif       

	Block reverse(DistanceType dist) {
	    Block result;
	    result.distance = dist;
	    for (int i=size()-1; i>=0; i--) 
		result.columns.push_back(columns[i]);
	    return result;
	}

	friend void initConstants();

	Double getPrefixThresh(bool complement, int i) const;
	Double getSuffixThresh(bool complement, int i) const;
	Double getPartialThresh(bool complement, int from, int to) const {
	    return complement ? 
		thresholdMatrix[size()-from][size()-to] :
		thresholdMatrix[to][from];
	}
	Double getPartialThresh(bool complement, int from) const {
	    return getPartialThresh(complement, from, size());
	}
	Double getThreshold() const {
	    return threshold;
	}
	bool hasIP() const {
	    return iP != 0;
	}
	double getIntronFreq(int i, int frame) const {
	    return iP->get(i,frame);
	}
	bool isAnchor() {
	    return threshold > getSpecThresh(min_anchor_spec); //getBackDist().abs(min_anchor_spec);
	}
// 	Dist getBlockscoreBack() {
// 	    return backDists[0];
// 	}
	// the following methods use logscores
	Dist getOwnDist(int from, int to, bool complement = false) const {
	    return complement ? 
		ownDists[size()-to] - ownDists[size()-from] : 
		ownDists[from] - ownDists[to];
	}
	Dist getOwnDist(int from=0) const {
	    return ownDists[from];
	}
	Dist getBackDist(int from, int to, bool complement = false ) const {
	    return complement ? 
		backDists[size()-to] - ownDists[size()-from]:
		backDists[from] - backDists[to];
	}
	Dist getBackDist(int from=0) const {
	    return backDists[from];
	}
	double ownNormalize(double logscore, int from, int to, bool complement = false) const {
	    return getOwnDist(complement, from, to).normed(logscore);
	}
	double backNormalize(double logscore, int from, int to, bool complement = false) const {
	    return getBackDist(complement, from, to).normed(logscore);
	}
	double getSensThresh(double stddev, int from, int to) const {
	    return getOwnDist(from, to).abs(-stddev);
	}
	double getSensThresh(double stddev, int from=0) const {
	    return getOwnDist(from).abs(-stddev);
	}
	double getSpecThresh(double stddev, int from, int to) const {
	    return getBackDist(from, to).abs(stddev);
	}
	double getSpecThresh(double stddev, int from=0) const {
	    return getBackDist(from).abs(stddev);
	}
	const vector<Dist>& getAllOwnDists() const {
	    return ownDists;
	}

	string id;
	DistanceType distance; // distance range to previous block

	static const Double almostZero;
	static const double invalidScore;

    private:
	Block() : iP(0) {}

	vector<Column> columns;
	ProfileMap* iP;
//	map<int,Double> hits[3]; // will be updated by addBlocksUntil
 
#ifdef DEBUG
	vector<Double> prefixThresh, suffixThresh;
#endif
	vector< vector<Double> > thresholdMatrix;
	vector<Dist> backDists, ownDists;
	// Dist blockscoreBack;
	Double threshold; // threshold for full block scores to be inserted in the hitlist
	static double min_spec;
	static double min_sens;
	static double min_anchor_spec;
	static double partial_spec;
	static double partial_sens;
    }; // class PP::Block

    
    class BlockScoreType {
    public: 
	typedef vector<Double> aligned_type;    // partial (prefix) scores for the block at one dna position
	                                        // back() is the full block score

	typedef deque<aligned_type> array_type; // the aligned_types for all dna positions

	BlockScoreType(const Block& block) : 
	    offset(0), partner(&block) {}
	bool addBlocksUntil(bool complement, int newbase, map<int,Double>* result=NULL);
	void removeBlocksUntil(int newbase) {
	    newbase -= 3 * size();
	    while (begin() < newbase)
		pop_front();
	}
	Double savedScore(int offset) {
	    map<int,Double>::iterator it = hits[offset%3].find(offset);
	    return it == hits[offset%3].end() ? 0 : it->second;
	}
	const map<int,Double>& savedHits(int f) const {
	    return hits[f];
	}
	void init(int newoffset = 0) { 
	    offset = newoffset;
	    values.clear();
	    for (int f=0; f<3; f++) hits[f].clear();
	}
	void pop_front() { // removes the first block from the array and returns its score
	    values.pop_front();
	    ++offset;
 	}
	Double front() {
	    return values.front().back();
	}
// 	void push_back(const aligned_type& new_vals) {
// 	    values.push_back(new_vals);
// 	}
// 	void push_front(const aligned_type& new_vals) {
// 	    values.push_front(new_vals);
// 	    offset--;
// 	}
	const aligned_type& getAlignedScores(int  x) const {
	    return values.at(x-offset);
	}
	Double get(int dna_offset, int block_offset, int len)  const { 
	    if (len==0)
		return 1;
	    if (block_offset == 0)
		return getAlignedScores(dna_offset).at(len-1);
	    const aligned_type& alignedScores = getAlignedScores(dna_offset - 3*block_offset);
	    return alignedScores.at(len+block_offset-1)/alignedScores.at(block_offset-1);
 	}
	Double get(int dna_offset, int block_offset)  const { 
	    if (block_offset == 0)
		return get(dna_offset);
	    const aligned_type& alignedScores = getAlignedScores(dna_offset - 3*block_offset);
	    return alignedScores.back()/alignedScores.at(block_offset-1);
 	}
	Double get(int dna_offset) const {
	    return getAlignedScores(dna_offset).back();
	}
	const int& begin() const {
	    return offset;
	}
	int count() const {
	    return values.size();
	}
	int end() const {
	    return offset + count();
	}
 	int size() const {
	    return partner->size();
	}
	// Double threshold() const {
	//     return partner->threshold;
	// }
	Double prefixThresh(bool complement, int i) const {
	    return partner->getPrefixThresh(complement, i);
	}
	Double suffixThresh(bool complement, int i) const {
	    return partner->getSuffixThresh(complement, i);
	}
   private:
	aligned_type& new_back() {
	    values.push_back(aligned_type());
	    return values.back();
	}
		

	// internal data fields
	array_type values;
	int offset; // dna pos of first base of first evaluated block
	const Block* partner;
	map<int,Double> hits[3]; // will be updated by addBlocksUntil
    }; // class PP::BlockScoreType


	

    class Profile {
	static Column background;
	static double lowest_odds, exp_lowest_odds;

    public:
	// methods of class ProteinProfile
	// initialise pattern from file <filename>
	Profile(string filename);
	
    

// 	static bool idInBlock(int id) {
// 	    return ((id >> INTERBLOCKBITS)) % 2 > 0;
// 	}

// 	Position getPosOfId(int id) const {
// 	    if (id < 0) return Position();
// 	    Position result(getB(id), getI(id));
// #ifdef DEBUG
// 	    if (result.b > blockCount()) {
// 		cout << "getB() > #blocks\n";
// 	    } else if (result.i > blockSize(result.b)) {
// 		cout << "getI() > blocksize\n";
// 	    }
// #endif
// 	    if (result.i < 0) result.i <<= skipBits[result.b];
// 	    return result;
// 	}
// 	Position getPosOfCopyId(int id, bool comp=false) const {
// 	    return getPosOfId(id + copyOffset);
// 	}
	Position firstInRange(int blockno=0) const {
	    return Position(blockno, -interBlockRange(blockno).max);
	}

	int blockCount() const {
	    return blocks.size();
	}
	int blockSize(int blockno) const {
#ifdef DEBUG
	    if (blockno < 0 || blockno >= blocks.size())
		throw ProjectError("prfl.blockSize() called with index out of range");
#endif
	    return blocks[blockno].size();
	}
	DistanceType interBlockDist(int blockno) const {
	    return blockno == blockCount() ? finalDist : blocks[blockno].distance;
	}
	Range interBlockRange(int blockno) const {
	    return interBlockDist(blockno).r;
	}
	// bool has_max(int blockno) const {
	//     return interBlockDist(blockno).has_max;
	// }
	// DistanceType interBlockFullDist(int blockno) const {
	//     if (blockno == 0)
	// 	return blocks[0].distance;
	//     DistanceType result = interBlockDist(blockno);
	//     result.r += blocks[blockno-1].size();
	//     return result;
	// }
// 	int interBlockMaxDist(int blockno) const {
// 	    return interBlockFullRange(blockno).max;
// 	}
// 	bool isInRange(int blockno, int len) const {
// 	    return interBlockRange(blockno).has(len);
// 	}
	Column getColumn(Position pos) const {
	    return blocks[pos.b][pos.i];
	}
	string getName() const {
	    return name;
	}
	string getInfoLine() const {
	    ostringstream result;
	    int blockno=0;
	    while (true) {
		result << "--";
		DistanceType d = interBlockDist(blockno);
		if (d.has_max)
		    result << "[" << d.r.min << ".." << d.r.max << "]";
		else if (d.r.min>0)
		    
		    result << "[" << d.r.min << ",...]";
		result << "--";
		if (blockno == blockCount())
		    return result.str();
		result << "> " << blocks[blockno].id << " (" << blockSize(blockno) << ") <";
		blockno++;
	    }
	}
	ostream& write(ostream&) const;
	
	friend ostream& operator<< (ostream& strm, const Profile& P) {
	    return P.write(strm);
	}
	
	// Position finalPos() const {
	//     return Position(blockCount(),0);
	// }
	// SubstateId getTerminalId() const {
	//     return Position(blockCount(),0).id();
	// }

	Block& operator[] (int blockno) {
	    return blocks[blockno];
	}
	const Block& operator[] (int blockno) const {
	    return blocks[blockno];
	}
	Double getGlobalThresh(bool complement, Position ppos) const {
	    return 
		ppos.b == blockCount() ? 1 :
		globalThresh[complement][ppos.b][ppos.i >= 0 ? ppos.i : 0];
	}
	
	friend void initConstants();

#ifdef DEBUG
	static bool vv_equal(const vector<Double>& v1, const vector<Double>& v2) {
	    return equal(v1.begin(), v1.end(), v2.begin(), almost_equal);
	}
#endif

#ifdef DEBUG
	Profile reverse() {
	    Profile result;
	    vector< vector<Dist> > ownDists;
	    result.name = name + "(rev.)";
	    for (int b=blockCount(); b>0; b--) {
		result.blocks.push_back(blocks[b-1].reverse(interBlockDist(b)));
		result.blocks.back().initThresholds();
		ownDists.push_back(result.blocks.back().getAllOwnDists());
	    }
	    result.calcGlobalThresh(ownDists);
	    if (!equal(result.globalThresh[0].begin(),
		       result.globalThresh[0].end(),
		       globalThresh[1].begin(),
		       vv_equal) ||
		!equal(result.globalThresh[1].begin(),
		       result.globalThresh[1].end(),
		       globalThresh[0].begin(),
		       vv_equal) )
		throw ProjectError("globalThreshs not equal!");
	    for (int b=0; b<blockCount(); b++) 
		for (int i=1; i<result.blockSize(b); i++)
		    if (!almost_equal(result[b].getPrefixThresh(false, i),
				      blocks[blockCount()-b-1].getPrefixThresh(true,i)) ||
			!almost_equal(result[b].getPrefixThresh(true, i),
				      blocks[blockCount()-b-1].getPrefixThresh(false,i)))
			throw ProjectError("prefixThreshs not equal!");
		    else if (!almost_equal(result[b].getSuffixThresh(false, i),
					   blocks[blockCount()-b-1].getSuffixThresh(true,i)) ||
			     !almost_equal(result[b].getSuffixThresh(true, i),
					   blocks[blockCount()-b-1].getSuffixThresh(false,i)))
			throw ProjectError("suffixThreshs not equal!");
	    result.finalDist = interBlockDist(0);
	    return result;
	}
#endif

    private:
	Profile() {}
	void calcGlobalThresh(const vector< vector<Dist> >&); 
	void parse_stream(istream&);

	// data fields 
	string name;
	vector<Block> blocks;  // the blocks
	vector< vector<Double> > globalThresh[2];
	DistanceType finalDist;       // final non-block range
	ProfileMap iP;      // not implemented

	static int    min_anchor_count;
	static double global_thresh;
	static double absolute_malus_threshold;
// 	// moved to SubstateModel
// 	void addSkipBits(Range r) {
// 	    int n=0; 
// 	    while (r.max > (1<<INTERBLOCKBITS)) {
// 		r.max >>= 1;
// 		n++;
// 	    }
// 	    skipBits[0].push_back(n);
// 	}
// 	int copyOffset;

    };  // class PP::Profile		
	    
    struct DNA {
	static const char* sequence;
	static int len;
	static void initSeq(const char* dna, int seqlen) {
	    sequence = dna;
	    len = seqlen;
	}
	static void initSeq(string s) {
	    sequence = s.c_str();
	    len = s.length();
	}
    };

    void initConstants();
} // namespace PP

#ifdef DEBUG
namespace std {
inline ostream& operator<< (ostream& strm, PP::Position ppos) {
    return strm << ppos.b << "/" << ppos.i;
}
}
#endif

inline Double PP::Block::getPrefixThresh(bool complement, int i) const {
#ifdef DEBUG
    if (i<=0 || i>=size()) 
	throw ProjectError("Out of range!");
    Double result = complement ? 
	suffixThresh[size() - i] :
	prefixThresh[i];
    if (!almost_equal(result,getPartialThresh(complement, 0, i)))
	throw ProjectError("partialThresh bug");
#endif
    return getPartialThresh(complement, 0, i);
}

inline Double PP::Block::getSuffixThresh(bool complement, int i) const {
#ifdef DEBUG
    if (i<=0 || i>=size()) 
	throw ProjectError("Out of range!");
    Double result=complement ?
	prefixThresh[size() - i] :
	suffixThresh[i];
    if (!almost_equal(result,getPartialThresh(complement, i)))
	throw ProjectError("partialThresh bug");
#endif
    return getPartialThresh(complement, i);
}
#endif