gvsu/cs677/pa4/maximum-parsimony.cc

#include <iostream>
#include <string>
#include <sys/time.h>
#include <mpi.h>
#include <stdlib.h>         /* exit() */
#include <limits.h>         /* INT_MAX */
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include "Sequence.h"
#include "PTree.h"
using namespace std;

static vector<Sequence *> sequences;
static int num_seqs;
static int my_rank;
static int world_size;          /* the number of processes */
static int lowest_cost = INT_MAX;
static PTree * lowest_cost_tree = NULL;

void permute(int * array, int index);
void usage();
string getsequence(int fd);
void evalPermutation(int * array);

void usage()
{
    cerr << "Usage: maximum-parsimony <sequence-file>" << endl;
    exit(42);
}

string getsequence(int fd)
{
    string s;
    char chr;
    bool saw_name = false;
    for (;;)
    {
        int bytes_read = read(fd, &chr, 1);
        if (bytes_read < 1)
        {
            if (s == "")
                break;
            else
            {
                cerr << "Unexpected EOF in the middle of a line!" << endl;
                exit(2);
            }
        }
        else if (chr == '\n')
        {
            break;
        }
        else if (chr == ' ' || chr == '\t')
        {
            saw_name = true;
        }
        else
        {
            if (saw_name)
                s += chr;
        }
    }
    return s;
}

int main(int argc, char * argv[])
{
    char * filename = NULL;

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
    MPI_Comm_size(MPI_COMM_WORLD, &world_size);

    for (int i = 1; i < argc; i++)
    {
        filename = argv[i];
    }

    if (filename == NULL)
        usage();

    int fd = open(filename, O_RDONLY);
    if (fd < 0)
    {
        cerr << "Couldn't open '" << filename << "'!" << endl;
        exit(1);
    }

    for (;;)
    {
        string seqstr = getsequence(fd);
        if (seqstr == "")
            break;
        Sequence * seq = new Sequence(seqstr.c_str());
        sequences.push_back(seq);
    }

    close(fd);

    num_seqs = sequences.size();
    int indices[num_seqs];
    for (int i = 0; i < num_seqs; i++)
        indices[i] = 0;
    permute(&indices[0], 0);

    /* after all processes have evaluated their trees, figure out
     * which process discovered the lowest cost tree */
    int processes_lowest_costs[world_size];
    MPI_Gather(&lowest_cost, 1, MPI_INT,
               &processes_lowest_costs[0], world_size, MPI_INT,
               0, MPI_COMM_WORLD);
    int min_process = 0;
    if (my_rank == 0)
    {
        cout << "Lowest costs determined by processes:" << endl;
        for (int i = 0; i < world_size; i++)
        {
            cout << processes_lowest_costs[i] << "  ";
            if (processes_lowest_costs[i] < processes_lowest_costs[min_process])
                min_process = i;
        }
        cout << endl;
    }
    MPI_Bcast(&min_process, 1, MPI_INT, 0, MPI_COMM_WORLD);
    if (my_rank == min_process)
    {
        cout << "Process " << my_rank << " lowest-cost tree:" << endl;
        cout << (*lowest_cost_tree) << endl;
    }
    delete lowest_cost_tree;

    MPI_Finalize();

    return 0;
}

void evalPermutation(int * array)
{
    static int invocation = 0;
    bool save_tree = false;
    PTree * tree;
    if (invocation % world_size == my_rank)
    {
        /* I am responsible for evaluating this permutation */
        vector<Sequence *> theseLeaves(num_seqs);
        for (int i = 0; i < num_seqs; i++)
            theseLeaves[i] = sequences[array[i] - 1];
        tree = new PTree(theseLeaves);
        int cost = tree->getCount();
        if (cost < lowest_cost)
        {
            lowest_cost = cost;
            if (lowest_cost_tree != NULL)
                delete lowest_cost_tree;
            lowest_cost_tree = tree;
            save_tree = true;
        }
    }
    if (!save_tree)
        delete tree;
    invocation++;
}

void permute(int * array, int index)
{
    static int level = -1;
    level++;
    array[index] = level;

    if (level == num_seqs)
    {
        evalPermutation(array);
    }
    else
    {
        for (int i = 0; i < num_seqs; i++)
        {
            if (array[i] == 0)
            {
                permute(array, i);
            }
        }
    }

    level--;
    array[index] = 0;
}