nojhan/paradiseo
Watch
1
0
Fork
You've already forked paradiseo
 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

git-svn-id: svn://scm.gforge.inria.fr/svnroot/paradiseo@1493 331e1502-861f-0410-8da2-ba01fb791d7f

Browse source
This commit is contained in:
wcancino 2009-03-11 11:07:24 +00:00
commit 8a63886dbf
46 changed files with 8759 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

37
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/CMakeLists.txt Normal file
View file

@ -0,0 +1,37 @@
######################################################################################
### 1) Include the sources
######################################################################################
INCLUDE_DIRECTORIES(${CMAKE_INSTALL_PREFIX}/include/paradiseo/eo)
INCLUDE_DIRECTORIES(${CMAKE_INSTALL_PREFIX}/include/paradiseo/moeo)
INCLUDE_DIRECTORIES(${CMAKE_INSTALL_PREFIX}/include/paradiseo/peo)
INCLUDE_DIRECTORIES(.)
INCLUDE_DIRECTORIES(/usr/local/include/GTL)
LINK_DIRECTORIES(${CMAKE_INSTALL_PREFIX}/lib)
LINK_DIRECTORIES(/usr/local/lib)
######################################################################################
### 2) Define PhyloMOEA target
######################################################################################
SET( PHYLOMOEA_SOURCES eigensolver.cpp
likelihoodcalculator.cpp
likoptimizer.cpp
matrixutils.cpp
parsimonycalculator.cpp
PhyloMOEA.cpp
probmatrixcontainer.cpp
ProbMatrix.cpp
RandomNr.cpp
Sequences.cpp
SubsModel.cpp
phylotreeIND.cpp
treeIterator.cpp
utils.cpp )
ADD_EXECUTABLE( PhyloMOEA ${PHYLOMOEA_SOURCES} )
TARGET_LINK_LIBRARIES(PhyloMOEA gsl gslcblas GTL eo eoutils ga moeo cma)
INSTALL( TARGETS PhyloMOEA RUNTIME DESTINATION bin)

57
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/ExceptionManager.h Normal file
View file

@ -0,0 +1,57 @@
/***************************************************************************
* Copyright (C) 2005 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#ifndef _ExceptionManager_H_
#define _ExceptionManager_H_
using namespace std;
class ExceptionManager
{
string msg;
short number;
public:
ExceptionManager(int n) {
number = n;
switch(number){
case 10:
msg = "Could not open file\n";
break;
case 11:
msg ="Incorrect file format\n";
break;
case 12:
msg = "Could not write to file\n";
break;
}
};
void Report()
{
cout << endl << msg << endl;
cout << "PhyloMOEA was finished because errors." << endl;
exit(1);
}
};
#endif

54
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/Patterns.h Normal file
View file

@ -0,0 +1,54 @@
/***************************************************************************
* Copyright (C) 2005 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef PATTERNS_H
#define PATTERNS_H
#include <Sequences.h>
#include <vector>
#include <string>
class Patterns
{
private:
int num_patterns;
double freqs[4];
std::vector<unsigned char*> pattern;
std::vector<int> pattern_count;
Sequences &seqs;
public:
// return a pattern content
Patterns(Sequences &S);
~Patterns();
unsigned char* operator[](int n) const { return pattern[n]; }
inline int count(int n) const { return pattern_count[n]; }
inline int count() const { return num_patterns; }
inline double *frequences() { return freqs; }
inline double frequence(int i) const { return freqs[i]; }
inline std::string pattern_name(int i) const { return seqs.seqname(i); }
inline int search_taxon_name(std::string &s) const { return seqs.search_taxon_name(s); }
inline bool is_ambiguous(short int l) const { return seqs.is_ambiguous(l); }
inline bool is_gap(short int l) const { return seqs.is_gap(l); }
inline bool is_undefined(short int l) const { return seqs.is_undefined(l); }
inline unsigned char* ambiguos_meaning(short int l) const { return seqs.ambiguos_meaning(l); }
void calculate_frequences();
int num_sequences();
};
#endif

250
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/PhyloMOEA.cpp Normal file
View file

@ -0,0 +1,250 @@
#include <eo>
#include <moeo>
#include <PhyloMOEO.h>
#include <PhyloMOEO_operators.h>
#include <PhyloMOEO_init.h>
#include <PhyloMOEO_eval.h>
#include <PhyloMOEO_archive.h>
#include <PhyloMOEOProbMatrixContainerUpdater.h>
#include <moeoNSGAII2.h>
#include <utils/moeoBestObjVecStat.h>
#include <utils/moeoAverageObjVecStat.h>
#include <PhyloMOEOPartitionStat.h>
#include <eoCountedFileMonitor.h>
#include <eoSingleFileCountedStateSaver.h>
#include <vectorSortIndex.h>
#include <utils.h>
#include <ctime>
#include <apply.h>
gsl_rng *rn2;
RandomNr *rn;
//Sequences *seq;
long seed;
//vector<phylotreeIND> arbores;
string datafile,usertree, expid, path, algotype;
double pcrossover, pmutation, kappa, alpha;
unsigned int ngenerations, popsize, ncats;
ofstream exp_data,evolution_data, best_media_scores, final_trees, final_pareto_trees, clades_pareto, clades_final,final_scores,pareto_scores;
int main(int argc, char *argv[])
{
welcome_message();
eoParser parser(argc, argv);
seed = parser.createParam((unsigned int)(time(NULL)), "seed", "Random Seed", 's',"Param").value();
popsize = parser.createParam((unsigned int)(50), "popSize", "Population size", 'n',"Param").value();
pcrossover = parser.createParam((double)(0.8), "pcross", "Crossover Rate", 'c',"Param").value();
pmutation = parser.createParam((double)(0.1), "pmut", "Mutation Rate", 'm',"Param").value();
ngenerations = parser.createParam((unsigned int)(500), "nGen", "Number of Generations", 'g',"Param").value();
kappa = parser.createParam((double)(4.0), "kappa", "Kappa value", 'k',"Param").value();
alpha = parser.createParam((double)(2.0), "alpha", "Alpha value", 'a',"Param").value();
ncats = parser.createParam((unsigned int)(4), "nCat", "Number of Categories", 'r',"Param").value();
datafile = parser.createParam(string(), "data", "Datafile", 'd',"Param").value();
usertree = parser.createParam(string(), "treef", "Treefile", 't',"Param").value();
path = parser.createParam(string(), "path", "Treefile", 'p',"Param").value();
algotype = parser.createParam(string("nsgaii"), "algo", "Algorith, Type", 'b',"Param").value();
ostringstream convert;
convert << seed;
expid = parser.createParam(convert.str(), "expid", "Experiment ID", 'e',"Param").value();
if( datafile.size()==0 )
{
parser.printHelp( cout );
return(-1);
}
string filename = path + datafile;
cout << "\n\nReading Sequence Datafile...";
Sequences seq(filename.c_str());
cout << " done.\n";
// calculate datafile
seq.calculate_patterns();
seq.calculate_frequences();
rn2 = gsl_rng_alloc(gsl_rng_default);
rn = new RandomNr(seed);
phylotreeIND templatetree( rn, seq, rn2);
ParsimonyCalculator parsi_calc(templatetree);
SubstModel modelHKY( seq, SubstModel::HKY85);
modelHKY.init();
modelHKY.set_kappa(kappa); // banco_grande
ProbMatrixContainer probmatrixs(modelHKY);
LikelihoodCalculator lik_calc(templatetree, modelHKY, probmatrixs,ncats);
lik_calc.set_alpha(alpha);
modelHKY.init();
PhyloEval byobj( parsi_calc, lik_calc );
Phyloraninit initializer(templatetree);
eoState state;
//eoPop <PhyloMOEO> &population = state.takeOwnership(eoPop<PhyloMOEO>(popsize, initializer));
eoPop<PhyloMOEO> population(popsize, initializer);
//state.registerObject( population );
cout << "\n\nReading Initial Trees...";
if( usertree.size() >0)
{
filename = path + usertree;
readtrees(filename.c_str(), population);
}
cout << " done.\n";
cout << "\n\nCreating output files...";
try{
filename = path + datafile + "_exp_param_" + expid + ".txt";
exp_data.open(filename.c_str());
exp_data.precision(15);
exp_data.setf(ios::fixed);
if( !exp_data.is_open() )
{
throw( ExceptionManager(12) );
}
cout << " done.\n";
}
catch ( ExceptionManager e )
{
e.Report();
}
// create the moea
save_exp_params(exp_data);
seq.save_seq_data(exp_data);
moeoAverageObjVecStat <PhyloMOEO> bestfit;
moeoBestObjVecStat <PhyloMOEO> avgfit;
eoPopStat<PhyloMOEO> popstats;
eoCountedFileMonitor media_scores( 2, path + datafile + "_media_scores_" + expid + ".txt", "\t", true,true );
media_scores.add( bestfit);
media_scores.add( avgfit) ;
eoCountedFileMonitor evolution_scores( 2, path + datafile + "_evolution_data_" + expid + ".txt", "\n", true,true );
evolution_scores.add( popstats);
//cout << "\n\nRunning NSGA-II ..." << endl;
eoGenContinue<PhyloMOEO> continuator(ngenerations);
eoCheckPoint<PhyloMOEO> cp(continuator);
eoValueParam<unsigned> generationCounter(0, "Gen.");
eoIncrementor<unsigned> increment(generationCounter.value());
cp.add(increment);
eoStdoutMonitor monitor(false);
monitor.add(generationCounter);
cp.add(monitor);
Phylomutate mutator;
Phylocross crossover;
eoSequentialOp<PhyloMOEO> operadores;
operadores.add(crossover,pcrossover);
operadores.add(mutator,pmutation);
PhyloMOEOProbMatrixContainerUpdater probmatrixupdater(probmatrixs);
cp.add( bestfit);
cp.add( avgfit);
cp.add( media_scores);
cp.add( evolution_scores );
cp.add( popstats);
cp.add( probmatrixupdater );
// apply<PhyloMOEO> ( byobj, population );
// population.printOn(cout);
if(algotype == "ibea")
{
moeoAdditiveEpsilonBinaryMetric < ObjectiveVector > metric;
moeoIBEA < PhyloMOEO > algo (cp, byobj, operadores, metric);
cout << "\n\nRunning IBEA ..." << endl;
algo(population);
}
else
{
moeoNSGAII < PhyloMOEO > algo (cp, byobj, operadores);
cout << "\n\nRunning NSGA-II ..." << endl;
algo(population);
}
cout << "\nCalculating Final Solutions...";
cout << " done\n";
PhyloMOEOFinalSolutionsArchive finalsolutions;
finalsolutions.operator()(population);
//remove_final_solutions( population );
// optimize remaining solutions
cout << "\nOptimizing tree branch lenghts...\n";
optimize_solutions( finalsolutions, lik_calc );
cout << "\nReevaluating individuals \n";
apply<PhyloMOEO> ( byobj, finalsolutions );
finalsolutions.save_scores(path + datafile + "_final_scores_" + expid + ".txt","#Final Solutions Scores");
finalsolutions.save_trees(path + datafile + "_final_trees_" + expid + ".txt");
cout << "\ndone \n";
// print the optimized solutions
//print_scores_pop( -2, population, evolution_data);
//print_scores_pop( -2, population, final_scores);
//save_trees(finalsolutions, final_trees);
cout << "\n\nCalculating Final Solutions clade support...";
PhyloMOEOPartitionStat splitstats;
splitstats(finalsolutions);
eoFileMonitor finalsplitstatsaver(path+datafile+"_clades_final_"+expid+".txt");
finalsplitstatsaver.add(splitstats);
finalsplitstatsaver();
//cout << splitstats.value() << endl;
//partition_map split_frequences;
//calculate_frequence_splits(finalsolutions,split_frequences);
cout << " done\n";
//save_partitions(splitstats.value(), clades_final);
//split_frequences.clear();
// remove dominate solutions
cout << "\nCalculating Pareto-optimal Solutions...";
PhyloMOEOParetoSolutionsArchive paretosolutions;
paretosolutions.operator()(finalsolutions);
paretosolutions.save_scores(path + datafile + "_pareto_scores_" + expid + ".txt","#Pareto Solutions Scores");
paretosolutions.save_trees(path + datafile + "_pareto_trees_" + expid + ".txt");
cout << " done\n";
// print final pareto trees
//save_trees( paretosolutions, final_pareto_trees);
cout << "\nCalculating Pareto-optimal Solutions clade support...";
splitstats(paretosolutions);
//calculate_frequence_splits(paretosolutions,split_frequences);
eoFileMonitor paretosplitstatsaver(path+datafile+"_clades_pareto_"+expid+".txt");
paretosplitstatsaver.add(splitstats);
paretosplitstatsaver();
// save_partitions(splitstats.value(), clades_pareto);
//split_frequences.clear();
cout << " done\n";
exp_data.close();
evolution_data.close();
pareto_scores.close();
final_scores.close();
best_media_scores.close();
final_trees.close();
final_pareto_trees.close();
clades_pareto.close();
clades_final.close();
gsl_rng_free(rn2);
// delete probmatrixs;
delete rn;
cout << "\nPhyloMOEA execution finishes !\n";
return 0;
}

83
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/PhyloMOEO.h Normal file
View file

@ -0,0 +1,83 @@
/***************************************************************************
* Copyright (C) 2008 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef PHYLOMOEO_H_
#define PHYLOMOEO_H_
#include "phylotreeIND.h"
#include <iostream>
class ObjectiveVectorTraits : public moeoObjectiveVectorTraits
{
public:
static bool minimizing (int i)
{
return true;
}
static bool maximizing (int i)
{
return false;
}
static unsigned int nObjectives ()
{
return 2;
}
};
typedef moeoRealObjectiveVector < ObjectiveVectorTraits > ObjectiveVector;
class PhyloMOEO : public MOEO< ObjectiveVector, double, double >
{
private:
phylotreeIND *tree;
void copy(const PhyloMOEO &other)
{
if(tree!=NULL)delete tree;
tree=NULL;
if(other.tree!=NULL)tree = new phylotreeIND(other.get_tree());
if(!other.invalidObjectiveVector())this->objectiveVector( other.objectiveVector() );
if(!other.invalidFitness())this->fitness( other.fitness() );
if(!other.invalidDiversity())this->diversity( other.diversity() );
}
public:
PhyloMOEO() : MOEO < ObjectiveVector, double, double >() { tree = NULL; }
PhyloMOEO(const PhyloMOEO &other) : MOEO < ObjectiveVector, double, double >(other) {
tree=NULL;
if(other.tree!=NULL)tree = new phylotreeIND(other.get_tree());
}
void set_tree_template(phylotreeIND &other) { if(tree!=NULL)delete tree; tree = other.clone(); };
void set_random_tree(phylotreeIND &other) { if(tree!=NULL)delete tree; tree = other.randomClone(); };
PhyloMOEO& operator= (const PhyloMOEO& other) {
copy(other);
return *this;
}
void readFrom (std::istream &_is)
{
string s;
_is >> s;
tree->read_newick2( s );
}
phylotreeIND & get_tree() const { return *tree; }
~PhyloMOEO() { if(tree!=NULL)delete tree; }
};
#endif

187
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/PhyloMOEOPartitionStat.h Normal file
View file

@ -0,0 +1,187 @@
/***************************************************************************
* Copyright (C) 2008 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef PHYLOMOEOPARTITONSTAT_H_
#define PHYLOMOEOPARTITONSTAT_H_
#include <moeoObjVecStat.h>
#include <PhyloMOEO.h>
//#include <utils.h>
struct part_info
{
int freq; double prob;
double dp, dl, tdp, tdl;
short int type; // 001=intermediate 1, 010=mp 2, 100=ml 4,
// 110= mp & ml 6, 101 = ml & int 5 011 mp & int 3 , 111 tree 7
bool inverse;
part_info(): freq(0), prob(0.0), dp(0), dl(0), tdp(0), tdl(0), type(0), inverse(false) {};
};
typedef std::map<string,struct part_info> partition_map;
class PhyloMOEOPartitionStat : public eoStat<PhyloMOEO, partition_map>
{
public :
using eoStat<PhyloMOEO, partition_map>::value;
PhyloMOEOPartitionStat(std::string _description = "Partition Stats")
: eoStat<PhyloMOEO, partition_map>(partition_map(), _description) {}
virtual void operator()(const eoPop<PhyloMOEO>& _pop){
doit(_pop); // specializations for scalar and std::vector
}
virtual std::string className(void) const { return "PhyloMOEOPartitionStat"; }
private :
// Specialization for pareto fitness
void doit(const eoPop<PhyloMOEO>& _pop)
{
value().clear();
calculate_frequence_splits( _pop, value() );
}
void calculate_frequence_splits( const eoPop<PhyloMOEO> &pop, partition_map &split_frequences)
{
bool ismp, isml;
isml = isml = false;
double dtp, dtl, maxdl, maxdp, dp, dl;
int n = pop.size();
double best_ml_score, best_mp_score;
string split, split_i;
graph::edge_iterator it, it_e;
//print_media_scores(pop, -1, best_l_idx, best_p_idx);
moeoBestObjVecStat<PhyloMOEO> best_inds;
best_inds(pop);
best_mp_score = (best_inds.value())[0];
best_ml_score = (best_inds.value())[1];
const PhyloMOEO &best_l = best_inds.bestindividuals(1);
const PhyloMOEO &best_p = best_inds.bestindividuals(0);
maxdl = maxdp = 0;
for(int i=0; i<n; i++)
{
dl = dtp = dtl = dp = 0;
isml = ismp = false;
const PhyloMOEO &sol = pop[i];
if(!sol.get_tree().splits_valid())sol.get_tree().calculate_splits_exp();
it = sol.get_tree().TREE.edges_begin();
it_e = sol.get_tree().TREE.edges_end();
// is the mp or ml tree
if( sol.objectiveVector().operator[](1) == best_ml_score) isml=true;
if( sol.objectiveVector().operator[](0) == best_mp_score) ismp= true;
dtl = sol.get_tree().compare_topology_2( best_l.get_tree()) / (2.0*(sol.get_tree().TREE.number_of_edges() - sol.get_tree().number_of_taxons()));
//sol->compare_topology_3( *best_l);
dtp = sol.get_tree().compare_topology_2( best_p.get_tree()) / (2.0*(sol.get_tree().TREE.number_of_edges() - sol.get_tree().number_of_taxons()));
//sol->compare_topology_3( *best_p);
while(it!=it_e)
{
if( sol.get_tree().is_internal( *it))
{
split = sol.get_tree().get_split_key( *it);
split_i = sol.get_tree().get_invert_split_key( *it);
if( split_frequences.find(split) != split_frequences.end() )
{
split_frequences[split].freq++;
split_frequences[split].prob += 1.0/n;
}
else if( split_frequences.find(split_i) != split_frequences.end() )
{
split_frequences[split_i].freq++;
split_frequences[split_i].prob += 1.0/n;
split = split_i;
}
else
{
split_frequences[split].freq = 1;
split_frequences[split].prob = 1.0/n;
}
if( isml )
split_frequences[split].type |= 4; // split belongs ml tree
if( ismp )
split_frequences[split].type |= 2; // split belongs mp tree
if( !(isml || ismp)) split_frequences[split].type |= 1; // split belong intermediate
dp = sol.objectiveVector().operator[](0) - best_mp_score;
dl = sol.objectiveVector().operator[](1) - best_ml_score;
split_frequences[split].dp += dp;
split_frequences[split].dl += dl;
maxdp = (dp > maxdp ? dp : maxdp);
maxdl = (dl > maxdl ? dl : maxdl);
split_frequences[split].tdp += dtp;
split_frequences[split].tdl += dtl;
}
++it;
}
//sol->invalidate_splits();
sol.get_tree().remove_split_memory();
}
partition_map::iterator it1 = split_frequences.begin();
partition_map::iterator it2 = split_frequences.end();
while(it1!=it2)
{
if(maxdp!=0) (*it1).second.dp /= (1.0*(*it1).second.freq*maxdp);
if(maxdl!=0)(*it1).second.dl /= (1.0*(*it1).second.freq*maxdl);
(*it1).second.tdp /= (1.0*(*it1).second.freq);
(*it1).second.tdl /= (1.0*(*it1).second.freq);
++it1;
}
}
};
std::string inverse_split( string split )
{
string s;
string::iterator it1 = split.begin();
string::iterator it2 = split.end();
while(it1!= it2)
{
s += *it1 == '*' ? '.' : '*';
++it1;
}
return s;
}
std::ostream & operator<<(std::ostream & _os, const partition_map & _map) {
partition_map::const_iterator it1 = _map.begin();
partition_map::const_iterator it2 = _map.end();
_os << _map.size() << endl;
_os.setf(ios::fixed);
while(it1!=it2)
{
_os <<
((*it1).second.inverse ? inverse_split((*it1).first) : (*it1).first)
<< '\t' << (*it1).second.freq << '\t' << (*it1).second.prob << '\t'
<< (*it1).second.dp << '\t' << (*it1).second.dl << '\t'
<< (*it1).second.tdp << '\t' << (*it1).second.tdl << '\t' << (*it1).second.type << endl;
++it1;
}
return _os;
}
std::istream & operator>>(std::istream & _is, const partition_map & _map) { std::cout << "not implemented\n"; return _is; }
#endif

35
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/PhyloMOEOProbMatrixContainerUpdater.h Normal file
View file

@ -0,0 +1,35 @@
/***************************************************************************
* Copyright (C) 2008 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef PHYLOMOEOPROBMATRIXCONTAINERUPDATER_H_
#define PHYLOMOEOPROBMATRIXCONTAINERUPDATER_H_
#include <probmatrixcontainer.h>
class PhyloMOEOProbMatrixContainerUpdater: public eoUpdater
{
public:
PhyloMOEOProbMatrixContainerUpdater(ProbMatrixContainer &_pm): pm(_pm) {}
void operator() () { pm.clear(); }
void lastCall() { operator()(); };
private:
ProbMatrixContainer &pm;
};
#endif

116
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/PhyloMOEO_archive.h Normal file
View file

@ -0,0 +1,116 @@
/***************************************************************************
* Copyright (C) 2008 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef PHYLOMOEO_ARCHIVE_H
#define PHYLOMOEO_ARCHIVE_H
#include <PhyloMOEO.h>
typedef moeoArchive<PhyloMOEO> PhyloMOEOArchive;
typedef PhyloMOEOArchive PhyloMOEOPFArchive;
class PhyloMOEOParetoSolutionsArchive:public moeoUnboundedArchive<PhyloMOEO>
{
public:
void save_trees( std::string filename, string title="" )
{
create_file( filename );
if(title.size()>0) os << title << endl;
os << size() << std::endl;
for(int i=0; i<size(); i++)
operator[](i).get_tree().printNewick(os);
os.close();
}
void save_scores( std::string filename, string title = "" )
{
create_file( filename );
if(title.size()>0) os << title << endl;
for(int i=0; i<size(); i++)
os << operator[](i) << std::endl;
os.close();
}
protected:
std::ofstream os;
private:
void create_file( std::string filename)
{
os.open(filename.c_str());
if (!os){
std::string str = "Could not open " + filename;
throw std::runtime_error(str);
}
}
};
class PhyloMOEODummyArchive:public PhyloMOEOParetoSolutionsArchive
{
public:
void update(const eoPop < PhyloMOEO > & _pop)
{
std::copy(_pop.begin(), _pop.end(), back_inserter(*this) );
}
};
class PhyloMOEOFinalSolutionsArchive: public PhyloMOEOParetoSolutionsArchive
{
public:
// overwrite the default operator
void operator()(const eoPop < PhyloMOEO > & _pop) { update(_pop); }
protected:
void update(const eoPop < PhyloMOEO > & _pop)
{
std::copy(_pop.begin(), _pop.end(), back_inserter(*this) );
moeoParetoObjectiveVectorComparator < ObjectiveVector > paretoComparator;
for(int i=0; i< size()-1; i++)
{
//cout << (i+1)/(p.currentSize()*1.0) << "%" << endl;
for(int j=i+1; j< size(); j++)
{
phylotreeIND &sol = operator[](i).get_tree();
phylotreeIND &sol_2 = operator[](j).get_tree();
if( sol.compare_topology_2( sol_2) == 0)
{
if( paretoComparator( operator[](i).objectiveVector(), operator[](j).objectiveVector() ) )
{
operator[](i) = back();
pop_back();
i--; j = size();
}
else
{
operator[](j) = back();
pop_back();
j--;
}
}
}
}
}
};
#endif

55
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/PhyloMOEO_eval.h Normal file
View file

@ -0,0 +1,55 @@
/***************************************************************************
* Copyright (C) 2008 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef PHYLOMOEO_EVAL_H_
#define PHYLOMOEO_EVAL_H_
#include <PhyloMOEO.h>
#include <parsimonycalculator.h>
#include <likelihoodcalculator.h>
class PhyloEval : public moeoEvalFunc < PhyloMOEO >
{
public:
PhyloEval( ParsimonyCalculator &x, LikelihoodCalculator &y): parseval(x), likeval(y) { }
void operator () (PhyloMOEO & _sol)
{
if (_sol.invalidObjectiveVector())
{
ObjectiveVector objVec;
//if(! (_sol.get_tree().splits_valid() ) ) _sol.get_tree().calculate_splits();
parseval.set_tree(_sol.get_tree());
likeval.set_tree(_sol.get_tree());
//objVec[0] = 1;
//objVec[1] = 1;
objVec[0] = parseval.fitch();
objVec[1] = -likeval.calculate_likelihood();
_sol.objectiveVector(objVec);
}
}
private:
LikelihoodCalculator &likeval;
ParsimonyCalculator &parseval;
};
#endif

51
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/PhyloMOEO_init.h Normal file
View file

@ -0,0 +1,51 @@
/***************************************************************************
* Copyright (C) 2008 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef PHYLOMOEO_INIT_H_
#define PHYLOMOEO_INIT_H_
#include <PhyloMOEO.h>
class Phyloraninit : public eoInit <PhyloMOEO>
{
public:
Phyloraninit(phylotreeIND &templatetree): tree(templatetree) { }
void operator()(PhyloMOEO &ind)
{
ind.set_random_tree(tree);
}
private:
phylotreeIND &tree;
};
class Phylonewickinit : public eoInit <PhyloMOEO>
{
public:
Phylonewickinit(std::string newick): tree_string(newick) { }
void operator()(PhyloMOEO &ind)
{
phylotreeIND &tree = ind.get_tree();
tree.read_newick2( tree_string );
}
private:
std::string tree_string;
};
#endif

61
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/PhyloMOEO_operators.h Normal file
View file

@ -0,0 +1,61 @@
/***************************************************************************
* Copyright (C) 2008 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef _PHYLOMOEO_OPERATORS_H_
#define _PHYLOMOEO_OPERATORS_H_
#include <PhyloMOEO.h>
class Phylocross : public eoQuadOp < PhyloMOEO >
{
public:
/**
* the class name (used to display statistics)
*/
std::string className() const { return "Phylocross"; } ;
bool operator()(PhyloMOEO & _tree1, PhyloMOEO & _tree2)
{
phylotreeIND parent_1(_tree1.get_tree());
phylotreeIND parent_2(_tree2.get_tree());
phylotreeIND &son1 = _tree1.get_tree();
phylotreeIND &son2 = _tree2.get_tree();
parent_1.export_subtree( son2 );
parent_2.export_subtree( son1 );
parent_1.invalidate_splits();
parent_2.invalidate_splits();
return true;
}
};
class Phylomutate : public eoMonOp <PhyloMOEO>
{
public:
std::string className() const { return "Phylomut"; }
bool operator() (PhyloMOEO & _tree)
{
_tree.get_tree().mutate(1);
return true;
}
};
#endif

178
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/ProbMatrix.cpp Normal file
View file

@ -0,0 +1,178 @@
/***************************************************************************
* Copyright (C) 2005 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#include "ProbMatrix.h"
#include <cmath>
#include <matrixutils.h>
gsl_cheb_series ** ProbMatrix::cs = NULL;
double f_pij(double x, void *p)
{
struct params *ij = (struct params *)p;
ij->pmatrix->set_branch_length(x);
return ij->pmatrix->p_ij_t(ij->i,ij->j);
}
ProbMatrix::ProbMatrix(SubstModel *m, double bl)
{
model = m;
branchlenght = bl;
calculate_derivate=false;
p = new double[4*4];
p1d = new double[4*4];
p2d = new double[4*4];
tmpexp = new double[4*4];
}
/*
void ProbMatrix::init()
{
struct params par;
double tmpbl = branchlenght;
int k=0;
// first call
gsl_function F;
F.function = f_pij;
F.params = &par;
if(cs == NULL)
{
std::cout << "initializaing polynomials...\n";
cs = new gsl_cheb_series*[16];
calculatepmatrix();
// calculate chevchevy polynomials
for(int i=0; i<4; i++)
for(int j=0;j<4;j++)
{
cs[k] = gsl_cheb_alloc(8);
par.pmatrix = this;
par.i = i; par.j = j;
gsl_cheb_init (cs[k], &F, 0.0, 1.0);
k++;
}
}
//restore the original branch lenght
branchlenght = tmpbl;
calculatepmatrixchev();
}*/
void ProbMatrix::init_derivate()
{
if(calculate_derivate)return;
calculate_p1d_matrix();
calculate_p2d_matrix();
calculate_derivate=true;
}
// set a new branch lenght and recalculate the probability matrix
void ProbMatrix::set_branch_length(double t)
{
// save time if t = older branch lenght
if ( branchlenght != t)
{
branchlenght = t;
calculatepmatrix();
calculate_p1d_matrix();
calculate_p2d_matrix();
}
calculate_derivate=false;
}
// calculate the matrix P = evec*diag( e^eigenvalues )*ievec
void ProbMatrix::calculatepmatrix()
{
double temp;
double **xevec = model->eigensystem->eigenvectors();
double *xeval = model->eigensystem->eigenvalues();
double **xievec = (double **)model->ievec;
//diag( e^eigenvalues )*ievec
for (int i = 0; i < 4; i++)
{
temp = exp( branchlenght * xeval[i] );
for (int j = 0; j < 4; j++)
tmpexp[i*4+j] = xievec[i][j]*temp;
}
// P matrix
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
temp = 0.0;
for (int k = 0; k < 4; k++)
temp += xevec[i][k] * tmpexp[k*4+j];
p[i*4+j] = fabs(temp);
}
}
}
void ProbMatrix::calculatepmatrixchev()
{
int k=0;
for(int i=0; i<4; i++)
for(int j=0; j<4; j++)
{
p[i*4+j] = gsl_cheb_eval(cs[i*4+j], branchlenght);
k++;
}
}
// matrix Q*P
void ProbMatrix::calculate_p1d_matrix()
{
double temp;
double *q=model->rate;
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
temp = 0.0;
for (int k = 0; k < 4; k++)
temp += q[i*4+k] * p[k*4+j];
p1d[i*4+j] = temp;
}
}
}
// calculate Q^2*P
void ProbMatrix::calculate_p2d_matrix()
{
double temp;
double *q=model->rate;
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
temp = 0.0;
for (int k = 0; k < 4; k++)
temp += q[i*4+k]*p1d[i*4+k];
p2d[i*4+j] = temp;
}
}
}

64
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/ProbMatrix.h Normal file
View file

@ -0,0 +1,64 @@
/***************************************************************************
* Copyright (C) 2005 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef PROBMATRIX_H
#define PROBMATRIX_H
#include "SubsModel.h"
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_chebyshev.h>
#include <map>
// store a set of the probability matric according some brach-lenght
class ProbMatrix
{
private:
SubstModel *model;
double branchlenght;
double *tmpexp;
static gsl_cheb_series **cs;
bool calculate_derivate;
void calculatepmatrix();
void calculatepmatrixchev();
void calculate_p1d_matrix();
void calculate_p2d_matrix();
public:
double *p, *p1d, *p2d;
ProbMatrix( SubstModel *m, double bl );
~ProbMatrix() { delete [] tmpexp; delete [] p; delete [] p1d; delete [] p2d; }
inline void init() { calculatepmatrix(); }
//void init();
void init_derivate();
inline double p_ij_t(int i, int j) { return p[i*4+j]; }
inline double p1d_ij_t(int i, int j) { return p1d[i*4+j]; }
inline double p2d_ij_t(int i, int j) { return p2d[i*4+j]; }
void set_branch_length(double l);
};
struct params { ProbMatrix *pmatrix; int i; int j; };
#endif

238
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/RandomNr.cpp Normal file
View file

@ -0,0 +1,238 @@
/*******************************************************************************
RandomNr.cpp
last change: 01/20/1999
version: 0.0.0
design: Eckart Zitzler
Paul E. Sevinc
implementation: Paul E. Sevinc
(c) 1998-1999: Computer Engineering and Networks Laboratory
Swiss Federal Institute of Technology Zurich
description:
See RandomNr.h
*******************************************************************************/
#include "RandomNr.h"
#include <climits>
#include <cstddef>
#include <ctime>
#include <cmath>
#include <iostream>
using namespace std;
static long idum2=123456789;
static long iy=0;
static long iv[NTAB];
static long idum=0;
RandomNr::RandomNr()
{
z = static_cast< long >( time( 0 ) ); // see <ctime>
iy = idum = 0;
sran1(z);
}
RandomNr::RandomNr( long seed )
: z( seed )
{
sran2(seed);
}
void
RandomNr::setSeed( long seed )
{
z = seed;
sran2(seed);
}
double
RandomNr::uniform01()
{
// see Reiser, Martin / Wirth, Niklaus.
// - Programming in Oberon: Steps beyond Pascal and Modula.
const int a = 16807,
m = LONG_MAX, // see <climits> - replace by
// numeric_limits< long >::max()
// when <limits> is available
q = m / a,
r = m % a;
long gamma;
gamma = a * ( z % q ) - r * ( z / q );
z = ( gamma > 0 ) ? gamma : ( gamma + m );
return z * ( 1.0 / m );
// NOTE: we're in trouble if at some point z becomes -LONG_MAX, 0, or LONG_MAX...
}
void
RandomNr::sran1(unsigned int seed) {
int j;
long k;
idum = seed;
if (idum == 0) idum=1;
if (idum < 0) idum = -idum;
for (j=NTAB+7;j>=0;j--) {
k=(idum)/IQ;
idum=IA*(idum-k*IQ)-IR*k;
if (idum < 0) idum += IM;
if (j < NTAB) iv[j] = idum;
}
iy=iv[0];
}
double
RandomNr::ran1() {
int j;
long k;
double temp;
k=(idum)/IQ;
idum=IA*(idum-k*IQ)-IR*k;
if (idum < 0) idum += IM;
j=iy/NDIV;
iy=iv[j];
iv[j] = idum;
if ((temp=AM*iy) > RNMX) return RNMX;
else return temp;
}
void
RandomNr::sran2(unsigned int seed) {
int j;
long k;
idum = static_cast<long>(seed);
if (idum == 0) idum=1;
if (idum < 0) idum = -idum;
idum2=(idum);
for (j=NTAB+7;j>=0;j--) {
k=(idum)/IQ1;
idum=IA1*(idum-k*IQ1)-k*IR1;
if (idum < 0) idum += IM1;
if (j < NTAB) iv[j] = idum;
}
iy=iv[0];
}
double
RandomNr::ran2() {
int j;
long k;
float temp;
k=(idum)/IQ1;
idum=IA1*(idum-k*IQ1)-k*IR1;
if (idum < 0) idum += IM1;
k=idum2/IQ2;
idum2=IA2*(idum2-k*IQ2)-k*IR2;
if (idum2 < 0) idum2 += IM2;
j=iy/NDIV;
iy=iv[j]-idum2;
iv[j] = idum;
if (iy < 1) iy += IMM1;
if ((temp=AM*iy) > RNMX) {
return RNMX;
}
else{
return temp;
}
}
bool
RandomNr::flipP( double p )
{
return ran2() < p ? true : false;
}
unsigned int
RandomNr::uniform0Max( unsigned int max ) // exclusive
{
return static_cast< unsigned short >( max * ran2() );
}
int
RandomNr::uniformMinMax( int min, // inclusive
int max ) // exclusive
{
return min >= 0
? static_cast< int >( min + ( max - min ) * ran2() )
: static_cast< int >( min - 1 + ( max - min ) * ran2() );
}
int
RandomNr::gaussMinMax( int min, // inclusive
int max ) // exclusive
{
double tmp = UnitGaussian();
if ((max*tmp > max) || ( max*tmp < min))
{
return max;
}
return min >= 0
? static_cast< int >( max * tmp )
: static_cast< int >( max * tmp );
}
double
RandomNr::doubleuniformMinMax( double min, // inclusive
double max ) // exclusive
{
return static_cast< double >( min + ( max - min ) * ran2() );
}
double RandomNr::UnitGaussian(){
static int cached=0;
static double cachevalue;
if(cached == 1){
cached = 0;
return cachevalue;
}
double rsquare, factor, var1, var2;
do{
var1 = 2.0 * ran2() - 1.0;
var2 = 2.0 * ran2() - 1.0;
rsquare = var1*var1 + var2*var2;
} while(rsquare >= 1.0 || rsquare == 0.0);
double val = -2.0 * log(rsquare) / rsquare;
if(val > 0.0) factor = sqrt(val);
else factor = 0.0; // should not happen, but might due to roundoff
cachevalue = var1 * factor;
cached = 1;
return (var2 * factor);
}
void RandomNr::resumo(){
cout << "semilla =" << z << endl;
}

116
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/RandomNr.h Normal file
View file

@ -0,0 +1,116 @@
/*******************************************************************************
RandomNr.h
last change: 01/20/1999
version: 0.0.0
design: Eckart Zitzler
Paul E. Sevinc
implementation: Paul E. Sevinc
(c) 1998-1999: Computer Engineering and Networks Laboratory
Swiss Federal Institute of Technology Zurich
description:
RandomNr is a helper class for pseudo-random number
generation that doesn't need to be subclassed unless
one wants to replace the algorithm used or add
another method.
Usually, one only instance per optimization problem
is necessary, and classes or functions that make use
of random numbers keep a reference to that instance.
*******************************************************************************/
#ifndef RANDOM_NR_H
#define RANDOM_NR_H
#include <cstddef>
//#include "TIKEAFExceptions.h"
#define IM1 2147483563L
#define IM2 2147483399L
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014L
#define IA2 40692L
#define IQ1 53668L
#define IQ2 52774L
#define IR1 12211L
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
#define IA 16807L
#define IM 2147483647L
//#define AM (1.0/IM)
#define IQ 127773L
#define IR 2836L
//#define NTAB 32
//#define NDIV (1+(IM-1)/NTAB)
//#define EPS 1.2e-7
//#define RNMX (1.0-EPS)
//typedef unsigned int size_t;
/*static long iy = 0;
static long iv[NTAB];
static long idum = 0;*/
class RandomNr
{
protected:
long z;
public:
// use the current time as seed
RandomNr();
// use the argument as seed
RandomNr( long );
void setSeed( long );
virtual ~RandomNr() {};
// return the next uniformly distributed
// random real in ]0; 1[
virtual double uniform01();
// return the next uniformly distributed
// random integer in [Min; Max[
int uniformMinMax( int, int );
// throw ( LimitsException );
double doubleuniformMinMax(double, double);
int gaussMinMax(int, int ); // exclusive
double UnitGaussian();
// return the next uniformly distributed
// random integer (size_t) in [0; Max[
unsigned int uniform0Max( unsigned int );
// throw ( LimitsException );
// flip an unfair coin (head is up with
// probability P), return true if head
// is up and false otherwise
bool flipP( double );
//throw ( ProbabilityException );
void sran1(unsigned int seed);
double ran1();
void sran2(unsigned int seed);
double ran2();
void resumo();
long get_seed() { return z; }
};
#endif

276
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/Sequences.cpp Normal file
View file

@ -0,0 +1,276 @@
/***************************************************************************
* Copyright (C) 2005 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#include <Sequences.h>
#include <iostream>
#include <fstream>
#include <map>
#include "assert.h"
#include <ExceptionManager.h>
void Sequences::list_taxons() const {
std::cout << "lista de taxons" << std::endl;
for(int i=0; i< seqnames.size(); i++)
std::cout << seqnames[i] << std::endl;
}
using namespace std;
Sequences::Sequences(const char *fn)
{
numseqs = 0;
seqlen = 0;
seqs = NULL;
fstream fin;
string tmpseq;
// ambiguos characters
meaning[0][0] = meaning[0][1] = 1; meaning[0][2] = meaning[0][3] = meaning[0][4] = 0; // M
meaning[1][0] = meaning[1][2] = 1; meaning[1][1] = meaning[1][3] = meaning[1][4] = 0; // R
meaning[2][0] = meaning[2][3] = 1; meaning[2][1] = meaning[2][2] = meaning[2][4] = 0; // W
meaning[3][1] = meaning[3][2] = 1; meaning[3][0] = meaning[3][3] = meaning[3][4] = 0; // S
meaning[4][1] = meaning[4][3] = 1; meaning[4][0] = meaning[4][2] = meaning[4][4] = 0; // Y
meaning[5][2] = meaning[5][3] = 1; meaning[5][0] = meaning[5][1] = meaning[5][4] = 0; // K
meaning[6][0] = meaning[6][1] = meaning[6][2] = 1 ; meaning[6][3] = meaning[6][4] = 0; // V
meaning[7][0] = meaning[7][1] = meaning[7][3] = 1 ; meaning[7][2] = meaning[7][4] = 0; // H
meaning[8][0] = meaning[8][2] = meaning[8][3] = 1 ; meaning[8][1] = meaning[8][4] = 0; // D
meaning[9][1] = meaning[9][2] = meaning[9][3] = 1 ; meaning[9][0] = meaning[9][4] = 0; // B
meaning[10][1] = meaning[10][2] = meaning[10][3] = meaning[10][0] = 0; meaning[10][4] = 1; // gap
meaning[11][1] = meaning[11][2] = meaning[11][3] = meaning[11][0] = 1; meaning[11][4] = 0; // undefined
// open the file for reading
try{
fin.open(fn, ios::in);
if(!fin.is_open())
{
cout << "\n" << fn << endl;
throw ExceptionManager(10);
}
// read numseq and seqlen
if(!(fin >> numseqs >> seqlen))
{
throw ExceptionManager(11);
}
seqs = new unsigned char[numseqs*seqlen];
seqnames.resize(numseqs);
for(int i=0; i < numseqs; i++)
{
fin >> seqnames[i];
fin >> tmpseq;
for(int j=0; j<seqlen; j++)
{
switch(toupper(tmpseq[j]))
{
case 'A': seqs[i*seqlen+j] = 0; break;
case 'C': seqs[i*seqlen+j] = 1; break;
case 'G': seqs[i*seqlen+j] = 2; break;
case 'T': seqs[i*seqlen+j] = 3; break;
case 'M': seqs[i*seqlen+j] = 4; break;
case 'R': seqs[i*seqlen+j] = 5; break;
case 'W': seqs[i*seqlen+j] = 6; break;
case 'S': seqs[i*seqlen+j] = 7; break;
case 'Y': seqs[i*seqlen+j] = 8; break;
case 'K': seqs[i*seqlen+j] = 9; break;
case 'V': seqs[i*seqlen+j] = 10; break;
case 'H': seqs[i*seqlen+j] = 11; break;
case 'D': seqs[i*seqlen+j] = 12; break;
case 'B': seqs[i*seqlen+j] = 13; break;
case '-': seqs[i*seqlen+j] = 14; break; // gap character
case '?':
case 'X':
case 'N' :
case 'O' : seqs[i*seqlen+j] = 15; break; // undefined character
default :
{
cout << "invalid character:" << tmpseq[j] << " seq:" << seqnames[i] << " column " << j << endl;
throw ExceptionManager(11);
assert(1==0); break; // invalid character
}
}
}
}
}
catch( ExceptionManager e)
{
e.Report();
}
fin.close();
}
int Sequences::search_taxon_name(std::string &s) const
{
int i;
for(i=numseqs-1; i>=0 && seqnames[i].substr(0,10)!=s; i--);
return i;
}
// return true is the site is informative
bool Sequences::is_informative(unsigned int col) const
{
unsigned char temp[16];
int sum = 0;
memset( temp, 0, 16*sizeof(unsigned char));
for(int j=0; j< numseqs; j++)
{
// check if the site is repeated
unsigned char l = (*this)[j][col]; //SeqData->pattern_pos(i,j);
if( temp[l] == 0)
{
temp[l] = 1; sum++;
}
}
// informative sites have at least two differents chars
return(sum > 1);
}
void Sequences::calculate_patterns()
{
// map pattern and index in the vector pattern
num_patterns=num_inf_sites=0;
map<string, int> pattern_map;
struct PatternInfo p;
string s;
for(int i=0; i< seqlen; i++)
{
s = position(i);
//if(s.find('?',0)== string::npos)cout << i << endl;
// new pattern
if( pattern_map.find(s) == pattern_map.end() )
{
pattern_map[s] = num_patterns;
p.idx = i;
p.count = 1;
patterns.push_back(p);
// calculate informative site
if( is_informative(i) )
{
inf_sites.push_back(num_patterns);
num_inf_sites++;
}
num_patterns++;
}
// old pattern
else
patterns[ pattern_map[s] ].count++;
}
//cout << "Padrones reconhecidos:" << num_patterns << endl;
//cout << "Sites informativos :" << num_inf_sites << endl;
}
void Sequences::calculate_frequences()
{
long num_undefined = 0; // number of undefined positions
freqs[0] = freqs[1] = freqs[2] = freqs[3] = 0.25;
double tmp[4];
double divisor;
// i don't know why most program do this
for(int m=0; m<8; m++)
{
tmp[0] = tmp[1] = tmp[2] = tmp[3] = 0;
for(int i=0; i< num_patterns; i++)
for(int j=0; j< numseqs; j++)
{
unsigned char l = pattern_pos(i,j);
if ( is_ambiguous(l) )
{
divisor = 0;
for(int k=0; k<4; k++)
divisor += ambiguos_meaning(l)[k] * freqs[k];
for(int k=0; k<4; k++)
tmp[k] +=
(ambiguos_meaning(l)[k]*freqs[k]/divisor)*
patterns[i].count;
}
else if( is_undefined(l) || is_gap(l) )
for(int k=0; k<4; k++) tmp[k] += freqs[k]*patterns[i].count;
else
tmp[l] += patterns[i].count;
}
for(int k=0; k<4; k++) freqs[k] = tmp[k] / (tmp[0]+tmp[1]+tmp[2]+tmp[3]);
}
// final frequences
//for(int i=0; i<4; i++){
// cout << "frequencia de " << i << " -->:" << freqs[i] << endl;
//}
// cout << "suma de frequencias:" << freqs[0] + freqs[1] + freqs[2] + freqs[3] << endl;
}
// get the
string Sequences::position(int pos)
{
string column;
unsigned char c;
char nucleotide;
for(int i=0; i < numseqs; i++)
{
c = seq_pos(i, pos);
switch(c)
{
case 0: nucleotide = 'A'; break;
case 1: nucleotide = 'C'; break;
case 2: nucleotide = 'G'; break;
case 3: nucleotide = 'T'; break;
case 4: nucleotide = 'M'; break;
case 5: nucleotide = 'R'; break;
case 6: nucleotide = 'W'; break;
case 7: nucleotide = 'S'; break;
case 8: nucleotide = 'Y'; break;
case 9: nucleotide = 'K'; break;
case 10: nucleotide = 'V'; break;
case 11: nucleotide = 'H'; break;
case 12: nucleotide = 'D'; break;
case 13: nucleotide = 'B'; break;
case 14: nucleotide = '-'; break;
case 15: nucleotide = '?'; break;
default: assert(1==0); break;
}
column+=nucleotide;
}
return column;
}
void Sequences::save_seq_data(ostream &of=cout)
{
of << "\nSequence Datafile Statistics\n";
of << "----------------------------------------------\n";
of << "Number of Sequences : " << this->numseqs << endl;
of << "Sequence Length : " << this->seqlen << endl;
of << "Number of Patterns : " << this->num_patterns << endl;
of << "Frequency A : " << this->freqs[0] << endl;
of << "Frequency C : " << this->freqs[1] << endl;
of << "Frequency G : " << this->freqs[2] << endl;
of << "Frequency T : " << this->freqs[3] << endl;
}

89
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/Sequences.h Normal file
View file

@ -0,0 +1,89 @@
/***************************************************************************
* Copyright (C) 2005 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef SEQUENCES_H
#define SEQUENCES_H
#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <cstring>
// class to store the sequences in a file
struct PatternInfo
{
int idx;
int count;
};
class Sequences
{
private:
//enum nucleotides = { A, C, G, T};
unsigned int num_patterns;
unsigned int num_inf_sites;
unsigned int numseqs;
unsigned int seqlen;
unsigned char *seqs;
std::vector<std::string> seqnames;
std::vector<struct PatternInfo> patterns;
std::vector<unsigned int> inf_sites;
unsigned char meaning[12][5];
double freqs[4];
bool is_informative(unsigned int col) const;
public:
inline unsigned char* operator[](int n) const { return &(seqs[n*seqlen]); }
Sequences(const char *filename);
~Sequences() {
delete [] seqs;
}
inline unsigned char seq_pos( int i, int j) const { return (*this)[i][j]; }
inline unsigned char pattern_pos(int i, int j) const
{ return seq_pos( j, patterns[i].idx); }
inline unsigned char infsite_pos(int i, int j) const
{ return pattern_pos( inf_sites[i], j); }
std::string position(int i);
inline int num_seqs() const { return numseqs; }
inline int seq_len() const { return seqlen; }
inline int pattern_count() const { return num_patterns; }
inline int pattern_count(int i) const { return patterns[i].count; }
inline int infsite_count() const { return num_inf_sites; }
inline int infsite_count(int i) const { return patterns[inf_sites[i]].count; }
inline double *frequences() { return freqs; }
inline double frequence(int i) const { return freqs[i]; }
inline std::string seqname(int i) const { return seqnames.at(i); }
inline bool is_defined(short int l) const { return l<4; }
inline bool is_ambiguous(short int l) const { return (l>3 && l<14); }
inline bool is_gap(short int l) const { return (l == 14); }
inline bool is_undefined(short int l) const { return (l ==15); }
inline std::vector<struct PatternInfo>& get_patterns() { return patterns; }
inline unsigned char* ambiguos_meaning(short int l) const { return
(unsigned char*) meaning[l-4]; }
void list_taxons() const;
int search_taxon_name(std::string &s) const;
void calculate_patterns();
void calculate_frequences();
void save_seq_data(std::ostream &of);
};
#endif

214
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/SubsModel.cpp Normal file
View file

@ -0,0 +1,214 @@
/***************************************************************************
* Copyright (C) 2005 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#include "SubsModel.h"
//#include <iostream>
using namespace std;
SubstModel::SubstModel(Sequences &p, int m)
{
model = m;
patterns = &p;
eigensystem = NULL;
patterns->calculate_frequences();
frequences = patterns->frequences();
kappa = 4;
//kappa = 32.451;
a=b=c=d=e=1;
ievec = new double*[4];
for(int i=0; i<4; i++)ievec[i] = new double[4];
}
void SubstModel::init()
{
init_rate_matrix();
switch(model)
{
case JC69:
jc69();break;
case F81:
f81();break;
case K2P:
k2p();break;
case HKY85:
hky85();break;
case GTR:
gtr();break;
default:
{
model = 1;
f81(); break;
}
}
// calculate eigenvalues and eigenvectors
if(eigensystem!=NULL)delete eigensystem;
eigensystem = new EigenSolver(rate, 4);
eigensystem->solve();
// calculat inverse eigenvectors
luinverse(eigensystem->eigenvectors(), ievec, NUM_AA);
}
void SubstModel::init_rate_matrix()
{
rate[0] = rate[1] = rate[2] = rate[3] =
rate[4] = rate[5] = rate[6] = rate[7] =
rate[8] = rate[9] = rate[10] = rate[11] =
rate[12] = rate[13] = rate[14] = rate[15] = 1.0;
}
void SubstModel::construct_rate_matrix()
{
// multiply the off-diagonal elements by the frequences
mult_frequences();
// calculate digonal elements
set_diagonal();
// normalize the rate matrix
normalize();
}
// multiple the off-diagonal elementes by the frequences
void SubstModel::mult_frequences()
{
for (int i = 0; i < 4; i++) {
for (int j = i+1; j < 4; j++) {
set_rate(i,j, get_rate(i,j) * frequences[j]);
set_rate(j,i, get_rate(j,i) * frequences[i]);
}
}
}
void SubstModel::set_diagonal()
{
// set the diagonal elements
for(int i=0; i<4; i++)
{
double sum = 0;
for(int j=0; j<4; j++) if(i!=j)sum+=get_rate(i,j);
set_rate(i,i, -sum);
}
}
// normaliza the rate matrix
void SubstModel::normalize()
{
double tmp = 0.0;
for (int i = 0; i < 4; i++)
{
tmp += -get_rate(i,i)*frequences[i];
}
//std::cout << "Rate matrix:" << std::endl;
//std::cout << "----------------------------------" << std::endl;
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
set_rate(i,j, get_rate(i,j)/tmp);
//std::cout << get_rate(i,j) << " | ";
}
//std::cout << std::endl;
}
}
void SubstModel::print_rate_matrix()
{
std::cout << "Rate matrix:" << std::endl;
std::cout << "----------------------------------" << std::endl;
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
std::cout << get_rate(i,j) << " | ";
}
std::cout << std::endl;
}
}
// useful for models that don't complain about nucleotide frequences
void SubstModel::set_equal_frequences()
{
frequences[0] = frequences[1] = frequences[2] = frequences[3] = 0.25;
}
// initialize the variables and matrizes
void SubstModel::f81()
{
// set the matrix
// parameter mu = 1.0
// nothing to do, only construct the matrix
construct_rate_matrix();
}
void SubstModel::jc69()
{
// set equal frequences
set_equal_frequences();
construct_rate_matrix();
}
void SubstModel::hky85()
{
//cout << "criando hky85 com kappa:" << kappa << endl;
// transition/traversion ratio
set_rate(0,2, kappa);
set_rate(1,3, kappa);
set_rate(2,0, kappa);
set_rate(3,1, kappa);
// multfrequences
construct_rate_matrix();
}
void SubstModel::k2p()
{
// set equal frequences
set_equal_frequences();
// the same of HKY85 except by the equal frequences
hky85();
}
void SubstModel::gtr()
{
a=1.23767538; b=3.58902963;
c=2.16811705; d=0.73102339; e=6.91039679;
set_rate(0,1, a);
set_rate(1,0, a);
set_rate(0,2, b);
set_rate(2,0, b);
set_rate(0,3, c);
set_rate(3,0, c);
set_rate(1,2, d);
set_rate(2,1, d);
set_rate(1,3, e);
set_rate(3,1, e);
// multfrequences
construct_rate_matrix();
}

74
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/SubsModel.h Normal file
View file

@ -0,0 +1,74 @@
/***************************************************************************
* Copyright (C) 2005 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef SUBSMODEL_H
#define SUBSMODEL_H
// substitution model class
// only nucleotides
#include "eigensolver.h"
#include "Patterns.h"
#include <iostream>
class SubstModel
{
private:
// instantaneus rate matrix 4x4
double rate[16]; // rate matrix
double **ievec; // inverse eigenvectors
double *frequences;
double a,b,c,d,e; // model parameterss
double kappa;
int model;
Sequences *patterns;
EigenSolver *eigensystem;
// define the substituion model
void construct_rate_matrix();
void setsubstmodel( int );
void calculatepmatrix(double);
void mult_frequences();
void set_diagonal();
void normalize();
void set_equal_frequences();
void init_rate_matrix();
void f81();
void k2p();
void jc69();
void hky85();
void gtr();
inline double get_rate(int i, int j) { return rate[4*i+j]; }
inline void set_rate( int i, int j, double val) { rate[4*i +j] = val; }
public:
enum models { JC69, F81, K2P, HKY85, GTR };
friend class ProbMatrix;
~SubstModel() { delete eigensystem; for(int i=0; i<4; i++) delete[] ievec[i]; delete [] ievec; }
SubstModel(Sequences &p, int m = 0);
void init();
void set_kappa(double t) { kappa = t; init(); }
double get_kappa() { return kappa; }
void print_rate_matrix();
};
#endif

77
contribution/branches/PhyloMOEA/PhyloMOEA-serial/PhyloMOEA/eigensolver.cpp Normal file
View file

@ -0,0 +1,77 @@
/***************************************************************************
* Copyright (C) 2005 by Waldo Cancino *
* wcancino@icmc.usp.br *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* 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, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#include "eigensolver.h"
EigenSolver::EigenSolver(double *matrix, int dim)
{
// allocate copy matrix
copymatrix = new double*[dim];
evec = new double*[dim];
eval = new double[dim];
dimension = dim;
for(int i=0; i<dim; i++)
{
copymatrix[i] = new double[dim];
evec[i] = new double[dim];
for(int j=0; j<dim; j++)
copymatrix[i][j] = matrix[i*dim+j];
}
}
EigenSolver::EigenSolver(double **matrix, int dim)
{
// allocate copy matrix
copymatrix = new double*[dim];
evec = new double*[dim];
eval = new double[dim];
dimension = dim;
for(int i=0; i<dim; i++)
{
copymatrix[i] = new double[dim];
evec[i] = new double[dim];
memcpy( copymatrix, matrix, dim*sizeof(double) );
}
}
EigenSolver::~EigenSolver()
{
// deallocate matrices
for(int i=0; i<dimension; i++)
{
delete [] evec[i];
delete [] copymatrix[i];
}
delete [] eval;
delete [] evec;
delete [] copymatrix;
}
void EigenSolver::solve()
{
int ordr[dimension];
double evali[dimension];
// resolve eigensystem
elmhes( copymatrix, ordr, dimension);
eltran( copymatrix, evec, ordr, dimension);
hqr2( dimension, 1, dimension, copymatrix, evec, eval, evali);
}

Some files were not shown because too many files have changed in this diff Show more