{

std::map<faststring,unsigned>::iterator it;

it = m.find(x);

if (it == m.end() )

{

m[x] = 1;

}

else

{

++it->second;

}

{

std::map<faststring,unsigned>::iterator it;

it = m.find(x);

if (it == m.end() )

{

return m[x] = 1;

}

else

{

return (++it->second);

}

{

size_t N = v.size();

for (size_t i=0; i< N; ++i)

{

v[i].shorten(len);

}

{

size_t N = v.size();

for (size_t i=0; i< N; ++i)

{

v[i].resize(len, fill);

}

{

unsigned char all_symbols_lookup[256];

for (unsigned j=0; j<256; ++j)

{

all_symbols_lookup[j] = (unsigned char)j;

}

faststring::size_type n = string_of_replacement_symbols.size();

for (faststring::size_type i=0; i<n; i+=2)

{

all_symbols_lookup[(unsigned char)string_of_replacement_symbols[i]] = string_of_replacement_symbols[i+1];

}

for (size_t i=0; i < v.size(); ++i)

{

char *pos = v[i].begin();

char *pos_end = v[i].end();

char *pos_write;

char c;

unsigned removed = 0;

size_t N = v[i].length();

pos_write = pos;

while (pos < pos_end)

{

c = all_symbols_lookup[(unsigned char)*pos];

if (c != skip_char)

{

*pos_write = c;

++pos_write;

}

else

{

++removed;

}

++pos;

}

v[i].shorten(N-removed);

}

{

unsigned digits = (unsigned)log10(v.size()) + 1;

unsigned mlen;

if (digits+1 >= len)

mlen = 0;

else

mlen = (unsigned)(len-digits-1);

std::map<faststring, unsigned> m_names, m_names_short;

std::map<faststring, unsigned>::iterator f_it, f_it2;

faststring tmp;

for (size_t i=0; i < v.size(); ++i)

{

add_or_count(m_names, v[i]);

}

for (size_t i=0; i < v.size(); ++i)

{

if (m_names[v[i]] > 1)

{

v[i].shorten(mlen);

}

add_or_count(m_names_short, v[i]);

}

m_names.clear();

char s = '_';

if (mlen == 0)

s = 'S';

unsigned n;

for (unsigned i=0; i < v.size(); ++i)

{

if (m_names_short[v[i]] > 1)

{

n = add_or_count_return(m_names, v[i]);

v[i] += s + faststring(n, '0', digits);

}

}

{

unsigned N = (unsigned)v.size();

for (unsigned i=0; i< N; ++i)

{

v[i].removeSpacesFront();

v[i].removeSpacesBack();

}

{

typedef char chartype;

SequenceDataTypesEnum datatype; // DNA, Protein, or other

unsigned taxaNum;

char ambig_char;

bool originalPosNumbers_supplied;

faststring::size_type posNum; // Stores 0 if sequences have unequal lengths.

// Checked for fasta. Phylip input guarantees

// equal lengths.

std::vector<faststring::size_type> originalPosNumbers; // In case we excluded positions from the alignment, we want to remember the original positions.

std::vector<CSequence_Mol*> seqData; // The vector of pointers to sequences

// Map of short names. Only makes sense if the short names are unique!

std::map<faststring, CSequence_Mol*> sn_map; // Obtain pointer to the sequence by sequence name. Is map of short names.

// Short names are normal names

bool short_names_unique; // The short_name is the name before the first space or the full name if the name has no spaces.

void add_seq(CSequence_Mol* seq)

{

// std::cerr << "Adding: " << seq->getName() << '\n';

if (sn_map.find(seq->getName()) != sn_map.end() )

short_names_unique = false;

seqData.push_back(seq);

sn_map[seq->getName()] = seq;

}

void determine_map_of_sequence_names()

{

short_names_unique = true;

// Creates/recreates the map of short names.

if (sn_map.size() == 0)

{

std::vector<CSequence_Mol*>::size_type N=seqData.size();

for(std::vector<CSequence_Mol*>::size_type i=0; i<N; ++i)

{

if (sn_map.find(seqData[i]->getName()) != sn_map.end() )

short_names_unique = false;

sn_map[seqData[i]->getName()] = seqData[i];

}

}

}

void recompute_are_shortnames_unique()

{

sn_map.clear();

determine_map_of_sequence_names();

}

// Minimal constructor: Empty sequences object

CSequences3(SequenceDataTypesEnum Dt, char set_ambig_char='?'):

datatype(Dt), taxaNum(0), ambig_char(set_ambig_char), originalPosNumbers_supplied(false),posNum(0), short_names_unique(true)

{}

// Constructor for a set of empty sequences with names and length.

CSequences3(SequenceDataTypesEnum Dt, std::vector<faststring> names, faststring::size_type len):

datatype(Dt), taxaNum((unsigned)names.size()), ambig_char('?'), originalPosNumbers_supplied(false),posNum(0), short_names_unique(true)

{

seqData.reserve(taxaNum);

CSequence_Mol *seq;

for (unsigned i=0; i<taxaNum; ++i)

{

seq = new CSequence_Mol (SeqType_dna, names[i], len, ambig_char);

add_seq(seq);

}

}

~CSequences3()

{

int n=taxaNum;

for (int i=0; i<n; ++i)

{

delete seqData[i];

}

}

// This constructor can be used as/to

// - general copy constructor

// - extract a range of sites

// Coordinates:

// pos1 must be the first column. Is 0 based index.

// pos2 must be the index after the last column. Is 0 based index.

// pos2-pos1 must be the number of bases that are copied to this sequence.

CSequences3(const CSequences3 &s, faststring::size_type pos1 = 0, faststring::size_type pos2 = faststring::npos):

datatype(s.datatype), taxaNum(s.taxaNum), ambig_char(s.ambig_char),

originalPosNumbers_supplied(false), posNum(0), short_names_unique(true)

{

unsigned n=taxaNum;

CSequence_Mol *seq;

if (n > 0)

{

// This is taken out of the loop to silence the warning that seq is used uninitialised below.

seq = new CSequence_Mol ( *(s.seqData[0]), pos1, pos2);

add_seq(seq);

for (unsigned i=1; i<n; ++i)

{

seq = new CSequence_Mol ( *(s.seqData[i]), pos1, pos2);

add_seq(seq);

}

// The following code does not check anything - remove later.

// if (i != n)

// {

// std::cerr << "Critical error in CSequences3 constructor: taxaNum and number of sequences found disagree.\n";

// }

posNum = seq->length();

}

}

bool are_short_names_unique()

{

return short_names_unique;

}

bool equal_length_of_all_sequences()

{

faststring::size_type len;

if (taxaNum==0)

return true;

len = seqData[0]->length();

for (unsigned i=1; i<taxaNum; ++i)

{

if (len != seqData[i]->length())

return false;

}

return true;

}

void trim_seq_names(faststring::size_type max_len)

{

for (unsigned i=0; i<taxaNum; ++i)

{

seqData[i]->trim_seq_name(max_len);

}

recompute_are_shortnames_unique();

}

void trim_seq_names(const char *trim_at_these_symbols)

{

for (unsigned i=0; i<taxaNum; ++i)

{

seqData[i]->trim_seq_name(trim_at_these_symbols);

}

recompute_are_shortnames_unique();

}

unsigned long memory_usage(unsigned long &value1, unsigned long &value2, unsigned long &value3, float &mean2)

{

unsigned long tmp1=sizeof(*this);

unsigned long tmp2=0, tmp3=0;

unsigned long i,N;

for (i=0, N=seqData.size(); i<N; ++i)

tmp2 += seqData[i]->memory_usage();

std::map<faststring, CSequence_Mol*>::iterator it;

it = sn_map.begin();

while (it != sn_map.end() )

{

tmp3 += it->first.memory_usage();

tmp3 += sizeof(it->first);

tmp3 += sizeof(it->second);

++it;

}

value1 = tmp1;

value2 = tmp2;

value3 = tmp3;

mean2 = (float)tmp2/N;

return tmp1 + tmp2 + tmp3;

}

CSequence_Mol* get_seq_by_name(faststring name) // Must be short name.

{

// Test code:

{

// print_DEBUG(cerr, 0);

}

std::map<faststring, CSequence_Mol*>::iterator find_it = sn_map.find(name);

if (find_it != sn_map.end() )

return find_it->second;

else

return NULL;

}

CSequence_Mol* get_seq_by_index(unsigned id)

{

if (id >= seqData.size())

{

return NULL;

}

else

{

return seqData[id];

}

}

// Depreciated: Use get_seq_by_index instead

/*

const char* get_Seq_Data(unsigned id)

{

if (id >= seqData.size())

{

return NULL;

}

else

{

return seqData[id]->getSeqStr();

}

}

const char* get_Seq_Name(unsigned id)

{

if (id >= seqData.size())

{

return NULL;

}

else

{

return seqData[id]->getName();

}

}

const char* get_Seq_FullName(unsigned id)

{

if (id >= seqData.size())

{

return NULL;

}

else

{

return seqData[id]->getFullName();

}

}

unsigned GetTaxaNum() { return taxaNum;}

faststring::size_type GetPosNum() { return posNum;}

SequenceDataTypesEnum get_datatype()

{

return datatype;

}

chartype GetChar(unsigned TaxaIndex, faststring::size_type PosIndex) const

{

assert(TaxaIndex < taxaNum);

assert(PosIndex < posNum);

return seqData[TaxaIndex]->get_pos(PosIndex);

}

void print_DEBUG(std::ostream &os, unsigned flag=1)

{

if (flag & 1) // scalars:

{

os << "Debug output CSequences3 object, flag==0\n";

os << "Data type: " << (int)datatype<< '\n';

os << "taxaNum: " << taxaNum << '\n';

os << "posNum: " << posNum << '\n';

os << "ambig_char: " << ambig_char << '\n';

os << "originalPosNumbers_supplied " << (int)originalPosNumbers_supplied << '\n';

os << "size of originalPosNumbers vector: " << originalPosNumbers.size() << '\n';

os << "size of seqData : " << seqData.size() << '\n';

os << "size of sn_map : " << sn_map.size() << '\n';

}

if (flag & 2)

{

faststring::size_type n=seqData.size();

os << "Data types of sequences:\n";

for (unsigned i=0; i<n; ++i)

{

os << i << ": " << seqData[i]->get_datatype() << " " << seqData[i]->type_as_string() << '\n';

}

}

}

/* void SetChar(unsigned TaxaIndex, unsigned PosIndex, chartype mychar) */

/* { */

/* assert(TaxaIndex < taxaNum); */

/* assert(PosIndex < posNum); */

/* seqData[TaxaIndex]->set_pos(PosIndex, mychar); */

/* } */

void SetOriginalPosNumber(faststring::size_type index, faststring::size_type theOriginalPosNumber)

{

originalPosNumbers_supplied = true;

assert(index < posNum);

originalPosNumbers[index] = theOriginalPosNumber;

}

faststring::size_type GetOriginalPosNumber(faststring::size_type index) const

{

assert(index < posNum);

if (originalPosNumbers_supplied)

return originalPosNumbers[index];

else

return UINT_MAX;

}

// Moves taxon with index index to top, preserving the order of all other taxa.

// Of course this changes all indices with the only exception that this taxon is

// already at the top of the vector.

// The index is of course 0 based.

void reorder_move_seq_to_top(unsigned index)

{

if (index < taxaNum && index > 0)

{

seqData.insert(seqData.begin(), seqData[index]);

seqData.erase(seqData.begin()+index+1); // We have to add 1 since we have one additional entry at the beginning

}

}

// TODO: Not very efficient!!

// Return true if taxon_name has been found, false otherwise

bool reorder_move_seq_to_top(faststring &taxon_name)

{

for (unsigned i=0; i< taxaNum; ++i)

{

faststring iname = seqData[i]->getFullName();

if (seqData[i]->getFullName() == taxon_name)

{

std::cerr << "Move to top: " << i << '\n';

reorder_move_seq_to_top(i);

return true;

}

}

return false;

}

// Reorder the sequences such the those in the given vector are at the top

// in the order given by the vector.

// TODO: Not very efficient! - Very simple and thus secure implementation.

bool reorder_sort_by(std::vector<faststring> &names)

{

std::vector<faststring*>::size_type n=names.size();

bool success = true;

for (int i=(int)n-1; i>=0; --i)

{

success = success & reorder_move_seq_to_top(names[i]);

}

return success;

}

// Reorder the sequences according to the

// in the order given by the vector.

// TODO: Not very efficient! - Very simple and thus secure implementation.

void reorder_sort_by_first_number_in_seqquence_names()

{

sort(seqData.begin(), seqData.end(), less_than_sequence_name_first_unsigned_in_fullname_pointer);

}

bool get_Originalposnumbers_Supplied()

{

return originalPosNumbers_supplied;

}

char get_ambiguity_character()

{

return ambig_char;

}

void trimSeqNamesAtSymbolsInString(const char *symbol_list)

{

for (unsigned i=0; i < taxaNum; ++i)

{

seqData[i]->trim_seq_name(symbol_list);

}

}

// Returns the sequence names upto the first space

void get_short_sequence_names(std::vector<faststring> &snames)

{

snames.clear();

snames.reserve(taxaNum);

// unsigned max=0;

for (unsigned i=0; i < taxaNum; ++i)

{

snames.push_back(seqData[i]->getName_faststring());

}

// return max;

}

// Returns the maximum sequence, length

void get_full_sequence_names(std::vector<faststring> &snames)

{

snames.clear();

snames.reserve(taxaNum);

// unsigned max=0;

for (unsigned i=0; i < taxaNum; ++i)

{

snames.push_back(seqData[i]->getFullName());

/* if (snames[i].size() > max) */

/* { */

/* max = snames[i].size(); */

/* } */

}

// return max;

}

// Only works for recoded sequences!!!!

void ConsensusSequence(faststring& conSeq,

const CSplit& setOfTaxa,

unsigned char numDistCharacters, // Distinguishable characters. Recoded characters have indices 0, ..., numDistCharacters-1

unsigned char numSymbols,

double consensusThreshold )

{

unsigned char maxindex, equalindex;

unsigned taxaCount = 0;

unsigned consensusSymMinimum;

std::vector<unsigned> counterSymbolsVec(numSymbols,0);

for (faststring::size_type pos=0; pos < posNum; ++pos)

{

// Initialise variable

for (unsigned char i=0; i < numDistCharacters; ++i)

counterSymbolsVec[i] = 0;

// Count number of occuring symbols for this position over all taxa

for (unsigned taxon = 0; taxon < taxaNum; ++taxon)

{

if (setOfTaxa.test(taxon))

{

++taxaCount;

++counterSymbolsVec[seqData[taxon]->get_pos(pos)];

}

}

consensusSymMinimum = (unsigned) std::ceil(consensusThreshold * taxaCount);

maxindex = 0;

equalindex = numDistCharacters;

for (unsigned char i = 1; i < numDistCharacters; ++i)

{

if (counterSymbolsVec[i] >= counterSymbolsVec[maxindex])

{

if (counterSymbolsVec[i] == counterSymbolsVec[maxindex])

equalindex = i;

maxindex = i;

}

}

if (counterSymbolsVec[maxindex] >= consensusSymMinimum &&

maxindex != equalindex)

conSeq[pos] = maxindex;

else // Default value, in case consensus cannot be determined

conSeq[pos] = numDistCharacters;

}

}

// This version works for any kind of sequence, not only for recoded sequences.

// See also special version for DNA sequences with additional features.

void ConsensusSequence(faststring& conSeq,

double consensusThreshold,

char default_char = '\0')

{

unsigned maxindex, equalindex;

unsigned taxaCount = 0;

unsigned consensusSymMinimum;

const char **seq_strs;

seq_strs = new const char* [taxaNum];

for (unsigned seq_index=0; seq_index<taxaNum; ++seq_index)

{

seq_strs[seq_index] = get_Seq_Data(seq_index);

}

conSeq.clear();

if (default_char == '\0')

{

if (get_datatype() == SeqType_dna && default_char == '\0')

{

default_char = 'N';

}

else if (get_datatype() == SeqType_protein)

{

default_char = 'X';

}

}

// std::cout << "Default char is: " << default_char << '\n';

unsigned *counterSymbolsVec = new unsigned[256];

for (faststring::size_type pos=0; pos < posNum; ++pos)

{

// std::cout << "--- \n";

// Initialise variable

taxaCount = 0;

std::memset(counterSymbolsVec, 0, 256*sizeof(unsigned));

unsigned char c;

// Count number of occuring symbols for this position over all taxa

for (unsigned taxon = 0; taxon < taxaNum; ++taxon)

{

// char c = toupper_lookup[(unsigned char)seqData[taxon]->get_pos(pos)];

c = (unsigned char)seq_strs[taxon][pos];

++counterSymbolsVec[(unsigned char)c];

// std::cout << "Counting << " << c << " in this column " << counterSymbolsVec[c] << '\n';

}

// All gaps

if (taxaNum == counterSymbolsVec[(unsigned char)'-'])

{

conSeq.push_back('-');

}

else

{

taxaCount = taxaNum - counterSymbolsVec[(unsigned char)'-'];

consensusSymMinimum = (unsigned) std::ceil(consensusThreshold * taxaCount);

maxindex = 0;

equalindex = UINT_MAX;

for (unsigned char i = 0; i < 254; ++i)

{

if (counterSymbolsVec[i] >= counterSymbolsVec[maxindex] && i != '-')

{

if (counterSymbolsVec[i] == counterSymbolsVec[maxindex])

equalindex = i;

maxindex = i;

}

}

// std::cout << "Pos: " << pos << " maxindex: " << (char) maxindex << " equalindex " << (char)equalindex << '\n';

if (counterSymbolsVec[maxindex] >= consensusSymMinimum && maxindex != equalindex)

conSeq.push_back(maxindex);

else // Default value, in case consensus cannot be determined

conSeq.push_back(default_char);

}

// std::cout << "conSeq: " << conSeq << '\n';

}

delete [] counterSymbolsVec;

}

void ConsensusSequence_DNA(faststring& conSeq,

double consensusThreshold,

unsigned minimum_uppercase_coverage,

unsigned minimum_total_coverage,

char default_char = '\0')

{

ConsensusSequence_DNA(conSeq, consensusThreshold, minimum_uppercase_coverage, minimum_total_coverage, NULL, default_char='\0');

}

void ConsensusSequence_DNA(faststring& conSeq,

double consensusThreshold,

unsigned minimum_uppercase_coverage,

unsigned minimum_total_coverage,

std::vector<unsigned> *Pointer_coverage_vec,

char default_char = '\0')

{

unsigned taxaCount = 0;

unsigned consensusSymMinimum;

const char **seq_strs;

seq_strs = new const char* [taxaNum];

for (unsigned seq_index=0; seq_index<taxaNum; ++seq_index)

{

seq_strs[seq_index] = get_Seq_Data(seq_index);

}

conSeq.clear();

if (Pointer_coverage_vec)

Pointer_coverage_vec->clear();

if (default_char == '\0')

{

default_char = 'N';

}

// std::cout << "Default char is: " << default_char << '\n';

unsigned *counterSymbolsVec = new unsigned[256];

std::memset(counterSymbolsVec, 0, 256*sizeof(unsigned));

for (faststring::size_type pos=0; pos < posNum; ++pos)

{

// std::cout << "--- \n";

// Initialise variable

taxaCount = 0;

unsigned char c;

// Count number of occuring symbols for this position over all taxa

for (unsigned taxon = 0; taxon < taxaNum; ++taxon)

{

// char c = toupper_lookup[(unsigned char)seqData[taxon]->get_pos(pos)];

c = (unsigned char)seq_strs[taxon][pos];

++counterSymbolsVec[(unsigned char)c];

// std::cout << "Counting << " << c << " in this column " << counterSymbolsVec[c] << '\n';

}

// All gaps

if (taxaNum == counterSymbolsVec[(unsigned char)'-'])

{

conSeq.push_back('-');

}

else

{

unsigned cA = counterSymbolsVec[(unsigned char)'A']+counterSymbolsVec[(unsigned char)'a'];

unsigned cC = counterSymbolsVec[(unsigned char)'C']+counterSymbolsVec[(unsigned char)'c'];

unsigned cG = counterSymbolsVec[(unsigned char)'G']+counterSymbolsVec[(unsigned char)'g'];

unsigned cT = counterSymbolsVec[(unsigned char)'T']+counterSymbolsVec[(unsigned char)'t'];

unsigned ctilde = counterSymbolsVec[(unsigned char)'~'];

// unsigned cN = counterSymbolsVec[(unsigned char)'N']+counterSymbolsVec[(unsigned char)'n'];

/*

taxaCount = cA + cC + cG + cT + ctilde;

consensusSymMinimum = (unsigned) std::ceil(consensusThreshold * taxaCount);

if (consensusSymMinimum < minimum_total_coverage)

consensusSymMinimum = minimum_total_coverage;

if (cA > cC && cA > cG && cA > cT && cA > ctilde)

{

if (cA >= consensusSymMinimum && counterSymbolsVec[(unsigned char)'A'] >= minimum_uppercase_coverage)

conSeq.push_back('A');

else

conSeq.push_back(default_char);

}

else if (cC > cA && cC > cG && cC > cT && cC > ctilde)

{

if (cC >= consensusSymMinimum && counterSymbolsVec[(unsigned char)'C'] >= minimum_uppercase_coverage)

conSeq.push_back('C');

else

conSeq.push_back(default_char);

}

else if (cG > cA && cG > cC && cG > cT && cG > ctilde)

{

if (cG >= consensusSymMinimum && counterSymbolsVec[(unsigned char)'G'] >= minimum_uppercase_coverage)

conSeq.push_back('G');

else

conSeq.push_back(default_char);

}

else if (cT > cA && cT > cC && cT > cG && cT > ctilde)

{

if (cT >= consensusSymMinimum && counterSymbolsVec[(unsigned char)'T'] >= minimum_uppercase_coverage)

conSeq.push_back('T');

else

conSeq.push_back(default_char);

}

else if (ctilde > cA && ctilde > cC && ctilde > cG && ctilde > cT)

{

if (ctilde >= consensusSymMinimum && ctilde >= minimum_uppercase_coverage)

conSeq.push_back('~');

else

conSeq.push_back(default_char);

}

else // No consensus, the larger two or more are equal

{

conSeq.push_back(default_char);

}

/*

} // END else of all gaps

// std::cout << "conSeq: " << conSeq << '\n';

if (Pointer_coverage_vec)

Pointer_coverage_vec->push_back(taxaCount);

// Reset the counter array:

counterSymbolsVec[(unsigned char)'A'] = 0;

counterSymbolsVec[(unsigned char)'a'] = 0;

counterSymbolsVec[(unsigned char)'C'] = 0;

counterSymbolsVec[(unsigned char)'c'] = 0;

counterSymbolsVec[(unsigned char)'G'] = 0;

counterSymbolsVec[(unsigned char)'g'] = 0;

counterSymbolsVec[(unsigned char)'T'] = 0;

counterSymbolsVec[(unsigned char)'t'] = 0;