-
-/*********************************************************************************************************************************/
-
-int PCOACommand::read_phylip(istream& f, int square_m, vector<string>& name_list, vector<vector<double> >& d){
- try {
- // int count1=0;
- // int count2=0;
-
- int rank;
- f >> rank;
-
- name_list.resize(rank);
- d.resize(rank);
- if(square_m == 1){
- for(int i=0;i<rank;i++)
- d[i].resize(rank);
- for(int i=0;i<rank;i++) {
- f >> name_list[i];
- // cout << i << "\t" << name_list[i] << endl;
- for(int j=0;j<rank;j++) {
- if (m->control_pressed) { return 0; }
-
- f >> d[i][j];
- if (d[i][j] == -0.0000)
- d[i][j] = 0.0000;
- }
- }
- }
- else if(square_m == 2){
- for(int i=0;i<rank;i++){
- d[i].resize(rank);
- }
- d[0][0] = 0.0000;
- f >> name_list[0];
- for(int i=1;i<rank;i++){
- f >> name_list[i];
- d[i][i]=0.0000;
- for(int j=0;j<i;j++){
- if (m->control_pressed) { return 0; }
- f >> d[i][j];
- if (d[i][j] == -0.0000)
- d[i][j] = 0.0000;
- d[j][i]=d[i][j];
- }
- }
- }
-
- return 0;
- }
- catch(exception& e) {
- m->errorOut(e, "PCOACommand", "read_phylip");
- exit(1);
- }
-
-}
-
-/*********************************************************************************************************************************/
-
-void PCOACommand::read(string fname, vector<string>& names, vector<vector<double> >& D){
- try {
- ifstream f;
- m->openInputFile(fname, f);
-
- //check whether matrix is square
- char d;
- int q = 1;
- int numSeqs;
- string name;
-
- f >> numSeqs >> name;
-
- while((d=f.get()) != EOF){
-
- //is d a number meaning its square
- if(isalnum(d)){
- q = 1;
- break;
- }
-
- //is d a line return meaning its lower triangle
- if(d == '\n'){
- q = 2;
- break;
- }
- }
- f.close();
-
- //reopen to get back to beginning
- m->openInputFile(fname, f);
- read_phylip(f, q, names, D);
- }
- catch(exception& e) {
- m->errorOut(e, "PCOACommand", "read");
- exit(1);
- }
-}
-
-/*********************************************************************************************************************************/
-
-void PCOACommand::recenter(double offset, vector<vector<double> > D, vector<vector<double> >& G){
- try {
- int rank = D.size();
-
- vector<vector<double> > A(rank);
- vector<vector<double> > C(rank);
- for(int i=0;i<rank;i++){
- A[i].resize(rank);
- C[i].resize(rank);
- }
-
- double scale = -1.0000 / (double) rank;
-
- for(int i=0;i<rank;i++){
- A[i][i] = 0.0000;
- C[i][i] = 1.0000 + scale;
- for(int j=i+1;j<rank;j++){
- A[i][j] = A[j][i] = -0.5 * D[i][j] * D[i][j] + offset;
- C[i][j] = C[j][i] = scale;
- }
- }
-
- A = linearCalc.matrix_mult(C,A);
- G = linearCalc.matrix_mult(A,C);
- }
- catch(exception& e) {
- m->errorOut(e, "PCOACommand", "recenter");
- exit(1);
- }
-
-}
-