///CalcNumOfSymmetries calculates the number of symmetries in a Newick.
///From http://www.richelbilderbeek.nl/CppCalcNumOfCombinations.htm
int CalcNumOfSymmetries(std::vector<int> v)
{
assert(IsNewick(v));
if (v.size() == 4) return (v[1] > 0 && v[1]==v[2] ? 1 : 0);
const int n_reserved
= *std::max_element(v.begin(),v.end())
+ std::count_if(v.begin(), v.end(),
std::bind2nd(std::greater<int>(),0));
int n_symmetries = 0;
int id = n_reserved + 1;
std::map<std::pair<int,int>,int> ids;
while (1)
{
//Count number of symmetries
assert(!v.empty());
const std::size_t sz = v.size();
assert(sz >= 2);
const std::size_t j = sz - 1;
for (std::size_t i = 0; i!=j; ++i)
{
if (v[i] > 0 && v[i]==v[i+1]) ++n_symmetries;
}
//Collect all leaves and store new leaves
std::vector<std::pair<int,int> > leaves;
for (std::size_t i = 0; i!=j; ++i)
{
if (v[i] > 0 && v[i+1] > 0)
{
//Keep pair sorted
const std::pair<int,int> p
= (v[i] <= v[i+1]
? std::make_pair(v[i+0],v[i+1])
: std::make_pair(v[i+1],v[i+0]));
//If this leaf is new, store it
if (ids.find(p)==ids.end())
{
ids[p] = id;
++id;
}
}
}
//Generalize all leaves
for (std::size_t i = 0; i < v.size()-1; ++i)
{
assert(v.size()>2);
if (v[i] > 0 && v[i+1] > 0)
{
//Keep pair sorted
const std::pair<int,int> p
= (v[i] <= v[i+1]
? std::make_pair(v[i+0],v[i+1])
: std::make_pair(v[i+1],v[i+0]));
//If this leaf is new, store it
assert(ids.find(p)!=ids.end()
&& "Leaf should have been stored already");
assert(i > 0);
std::vector<int> v_new;
std::copy(v.begin(),v.begin() + i - 1,std::back_inserter(v_new));
const int id = ids[p];
v_new.push_back(id);
std::copy(v.begin() + i + 3,v.end(),std::back_inserter(v_new));
v = v_new;
i = (i-1 > 0 ? i-1 : 0);
//Check if there are more leaves to be generalized
if (v.size()<=4)
{
//Check if the last (X,Y) is symmetrical...
return n_symmetries + (v[1] > 0 && v[1]==v[2] ? 1 : 0);
break;
}
}
}
}
}
|
