#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <vector>
#include <string>
#include <fstream>
#include <cstdlib>
#include <sstream>
#include <iostream>
#include <unistd.h>
#include <limits.h>
#include <assert.h>
#include <set>
#include <algorithm>
#include <limits>
#include <time.h>
#include <map>

/*
  A graph package for representing protein-protein interactions.

  The input file has the following format
  <number of nodes>
  comment
  
  node name, which is the gene name in the original Ito data
  node information, which is the (Bait/Prey) description in Ito data
  node number, which is a number assigned to the node uniquely 
  neighboring node numbers, space separated

  <repeated>
*/

using namespace std;
typedef unsigned int uint;
typedef vector<uint> SubGraph;

struct SubGraphs
{
    int n; //number of subgraphs
    vector<SubGraph> subgraphs;
    vector<int> L; //internal links
    vector<double> P0; //null probability
};

bool verbose = false;
bool summarize = false;
bool prune = false;

//store the pairwise_alignment of two subgraphs
struct pairwise_alignment
{
    SubGraph a;
    SubGraph b;
    int score; //score of this alignment
};

/* We represent the graph as a vector of vectors.
   graph[i] is the adjacency list for node i.
   All our graphs are undirected so if j is in
   graph[i], i is in graph[j]

   These could probably all be encapsulated inside a class..
   */

class Graph
{
    typedef unsigned int uint;
    vector< vector<bool> > matrix;
    vector< vector<uint> > adjlist;
    vector<string> node_name;
    vector<string> node_info;
    string file_info;
    int size;
    int num_edges;
    
public:
    Graph(uint n, string info): matrix(n), adjlist(n),
	node_name(n), node_info(n), file_info(info), size(n), num_edges(0)
    {
	//initialize n x n matrix
	for(int i = 0; i < n; i++)
	    matrix[i] = vector<bool>(n,false);
    }

    ~Graph() {}

    //set the name and info for node n
    void setNode(uint n, string name, string info)
    {
	assert(n < size);
	node_name[n] = name;
	node_info[n] = info;
    }

    //return reference to node n's name
    string& getName(uint n)
    {
	assert(n < size);
	return node_name[n];
    }

    //return reference to node n's info
    string& getInfo(uint n)
    {
	assert( n < size);
	return node_info[n];
    }

    //return the degree of node n
    int degree(uint n) const
    {
	return adjlist[n].size();
    }
    
    //insert directed edges i->j into graph
    //does nothing if edge already exists
    void addEdge(uint i, uint j)
    {
	assert(i < size);
	assert(j < size);

	if(!matrix[i][j])
	{ //edge not already present
	    matrix[i][j] = true;
	    //keep adjlist sorted
	    vector<uint>::iterator itr;
	    for(itr = adjlist[i].begin(); itr != adjlist[i].end(); itr++)
	    {
		if(*itr > j)
		    break;
	    }
	    adjlist[i].insert(itr, j);
	    num_edges++;
	}
    }

    //returns vector of neighbor nodes to n
    const vector<uint>& neighbors(uint n) const
    {
	assert(n < size);
	return adjlist[n];
    }

    //return true if there is an edge i->j
    bool isEdge(uint i, uint j) const
    {
	assert(i < size);
	assert(j < size);
	return matrix[i][j];
    }

    //return number of nodes
    int numNodes() const
    {
	return size;
    }

    //return number of (undirected)
    int numEdges() const
    {
	return num_edges;
    }
    
    //returns true if graph is undirected
    bool isUndirected() const
    {
	for(int i = 0; i < size; i++)
	    for(int j = 0; j < i; j++)
	    {
		if(matrix[i][j] != matrix[j][i])
		    return false;
	    }
	return true;
    }


    // compute the probability of each subgraph against the null ensemble
    // which is the unbiased sum of all graphs with the same connectivity
    // this is approximated by w(i,j) = deg(i)*deg(j)/(number of edges)
    // It is unclear to me whether we use the connectivity information
    // from the whole graph, or just the subgraph, so I'm going to use the
    // whole graph
    double null_probability(SubGraph& g)
    {
	double result = 1.0;

	//for every pair of nodes in the subgraph
	for(int i = 0; i < g.size(); i++)
	{
	    for(int j = 0; j < i; j++)
	    {
		int ni = g[i];
		int nj = g[j];
		int cij = isEdge(ni,nj);
		double wij = ((double)degree(ni)*degree(nj))/(numEdges()*2.0);

		if(cij)
		    result *= wij;
		else
		    result *= (1.0 - wij);
	    }
	}
	return result;
    }

    int subgraph_cnt; //number of total subgraphs

    //given a subgraph sub, will find all connected subgraphs containing
    //sub of one size larger not using nodes from marked.  Then calls itself
    //recursively on the result.  If the size of sub equals k, we add it to subgraphs.
    //A node is marked if it is in sub or if we've already found all connected subgraphs
    //which contain that node AND sub
    void compute_subgraphs(set<SubGraph >& subgraphs, int k, int l, vector<uint>& sub, vector<bool> marked,
	    map<double, SubGraph>& threshold_map, int t)
    {

	if(l == k) //terminate, not interested in larger graphs
	{
	   // printf("%d\r",cnt);
	    vector<uint> sorted_sub = sub; //don't mutate sub!
	    sort(sorted_sub.begin(),sorted_sub.end());

	    if(prune) //prune more
	    {
		//only accept subgraphs where all the nodes
		//have more than one internal link
		for(int i = 0; i < sub.size(); i++)
		{
		    int linkcnt = 0;
		    for(int j = 0; j < sub.size(); j++)
		    {
			if(i != j && isEdge(sub[i],sub[j]))
			    linkcnt++;
		    }
		    if(linkcnt <= 1)
			return;
		}
	    }
	    subgraph_cnt++;
	    
	    if(t > 0) //thresholding on
	    {
		double P0 = null_probability(sorted_sub); //compute probability
		if(threshold_map.size() >= t)
		{
		    //see if it's better than what we already have
		    if(P0 < (*threshold_map.rbegin()).first)
		    {
			subgraphs.erase((*threshold_map.rbegin()).second);
			threshold_map.erase((*threshold_map.rbegin()).first);
			subgraphs.insert(sorted_sub);
			threshold_map[P0] = sorted_sub;
			//note that if many graphs have the same P0 we don't
			//do exactly the right thing here, but that shouldn't cause too much trouble
		    }
		    //otherwise too likely in a random graph..
		}
		else
		{
		    subgraphs.insert(sorted_sub);
		    threshold_map[P0] = sorted_sub;
		}
	    }
	    else //just insert
	    {
		subgraphs.insert(sorted_sub);
	    }
	}
	else //find all subgraphs containing sub of size l+1
	{
	    //for each node of sub
	    for(int i = 0; i < l; i++)
	    {
		uint n = sub[i];
		//for each neighbor
		vector<uint> neigh = neighbors(n);
		for(int j = 0; j < neigh.size(); j++)
		{
	        //if it isn't marked, mark and add it
		    int m = neigh[j];
		    if(!marked[m])
		    {
			marked[m] = true;
			sub[l] = m;
			compute_subgraphs(subgraphs, k, l+1, sub, marked,threshold_map,t);
		    }
		}
	    }
	}
    }


    // find all the connected subgraphs of graph of size k and put them
    // into subgraphs.  A subgraph is just a sorted vector of node numbers
    // We use a recursive algorithm that given all the subgraphs of size n,
    // finds all the subgraphs of size n+1
    // If t > 0, return only the t many subgraphs with the highest probability
    int subgraphs(SubGraphs& subs, int k, int t)
    {
	set<SubGraph> setsubs; //we use a set to make detecting duplicates easy
	int n = size;
	vector<bool> marked(n,false); //keep track of nodes we've seen
	vector<uint> sub(k);
	subs.subgraphs.clear();

	//keep top t subgraphs sorted by probability
	map<double, SubGraph> threshold_map;
	
	if(k <= 0)
	    return 0;

	subgraph_cnt = 0;
	//base case.. single node as a subgraph of size 1
	for(int i = 0; i < n; i++)
	{
	    sub[0] = i;
	    marked[i] = true;
	    compute_subgraphs(setsubs, k, 1, sub, marked, threshold_map, t);
	}

	//copy setsubs into subs
	//We want the final version in vector form so we
	//can refer to specific subgraphs via their index

	subs.subgraphs.resize(setsubs.size());
	subs.n = subs.subgraphs.size();
	int i;
	set<SubGraph>::iterator itr;
	for(i = 0, itr = setsubs.begin(); itr != setsubs.end(); itr++, i++)
	{
	    subs.subgraphs[i] = *itr;
	}

	subs.L.resize(subs.n);
	//compute the number of internal links in each subgraph and store into L
	//for each subgraph
	for(int g = 0; g < subs.n; g++)
	{
	    int cnt = 0;
	    //for every pair of nodes in the subgraph
	    for(int i = 0; i < subs.subgraphs[g].size(); i++)
	    {
		for(int j = 0; j < i; j++)
		{
		    int ni = subs.subgraphs[g][i];
		    int nj = subs.subgraphs[g][j];		    
		    if(isEdge(ni,nj))
			cnt++;
		}
	    }
	    subs.L[g] = cnt;
	}

	subs.P0.resize(subs.n);

	for(int g = 0; g < subs.n; g++)
	{

	    double result = null_probability(subs.subgraphs[g]);
	    subs.P0[g] = result;
	}
    
	return subgraph_cnt;
    }

};


//return the score for the alignment between a and b of graph
//where we align a[i] with b[i]
//This is equation 3 from Lassig.  Each edge contained in a not
//in b and vice versa increments the score by 1.
//Since our graphs are all undirected, we count each edge only once
int score(const Graph& graph, const SubGraph& a, const SubGraph& b)
{
    //we assume the size of subgraphs is small, so we just look
    //at all k^2 possible edges rather than just looking at existing
    //edges in a and b and matching them
    int k = a.size();
    int ret = 0;
    
    for(int i = 0; i < k; i++)
    {
	for(int j = i; j < k; j++)
	{
	    if(graph.isEdge(a[i],a[j]) != graph.isEdge(b[i],b[j]))
	    {
		ret += 1;
	    }
	}
//	if(graph.degree(a[i]) != graph.degree(b[i]))
//	    ret += abs(graph.degree(a[i]) - graph.degree(b[i]));
    }
    
    return ret;
}

//Given a graph and a collection of subgraphs of that graph,
//compute the minimal mismatch score for each pair of subgraphs
//This requires us to compare all factorially many alignments of
//each pair of subgraphs.
//mismatches is returned as a lower triangular matrix
//We currently aren't storing the actual alignment that results
//in the minimal score, but that could be fixed easily enough.
void compute_pairwise_mismatches(const Graph& graph, const SubGraphs& subs, vector<vector<pairwise_alignment> >& mismatches)
{
    int p = subs.subgraphs.size(); //number of subgraphs
    //initialize mismatches
    mismatches.clear();
    mismatches.resize(p);
    for(int i = 1; i < p; i++)
	mismatches[i].resize(i);

    //consider each pair of subgraphs
    for(int i = 0; i < p; i++)
    {
	for(int j = 0; j < i; j++)
	{
	    //try all possible alignments and choose the minimum alignment
	    int min = INT_MAX;
	    SubGraph a = subs.subgraphs[i];
	    SubGraph b = subs.subgraphs[j];
	    SubGraph best_b;
	    
	    //keep permuting b
	    do
	    {
		int s = score(graph, a, b);
		if(s < min)
		{
		    best_b = b;
		    min = s; //found a better alignment
		}
	    }
	    while(next_permutation(b.begin(), b.end()));
	    mismatches[i][j].score = min;
	    mismatches[i][j].a = a;
	    mismatches[i][j].b = best_b;
	}
    }
}


//generate a random bitvector into s
void random_alignment(vector<bool>& s)
{
    int n = s.size();
    for(int i = 0; i < n; i++)
    {
	s[i] = rand()&1; //even or odd
    }
}


//Use equations 9/10 to score a given multiple alignment
//For now omitting mysterious logZ factor..
double score_alignment(const Graph& graph,
	               SubGraphs& subs,
		       vector<vector<pairwise_alignment> >& alignments,
		       vector<bool> s, //alignment to score
		       double sigma, //sigma parameter -> increases number of internal links
		       double mu //mu -> decreases fuzziness
		       )
{
    double result = 0;
    int p = 0; //number of subgraphs in this alignment
    //compute
    // sigma*sum(L) - mu/(2p) * sum(minimal pairwise alignments) - logZ

    //first sum(L)
    for(int i = 0; i < s.size(); i++)
    {
	if(s[i])
	{
	    p++;
	    result += subs.L[i];
	}
    }
    //multiply by sigma
    result *= sigma;

    //now compute sum over pairwise alignments
    double align_sum = 0.0;

    for(int i = 0; i < s.size(); i++)
    {
	for(int j = 0; j < i; j++)
	{
	    //only consider pairwise alignments that between
	    //members of the multiple alignment
	    if(s[i] && s[j])
	    {
		align_sum += alignments[i][j].score;
	    }
	}
    }

    //multiply be mu/2p
    align_sum *= (mu/(2*p));

    //subtract
    result -= align_sum;

    //TODO: subtract logZ

    return result;
}


//use simulated annealing to find a maximal multiple alignment
//score over all the subgraphs
double simulated_annealing(const Graph& graph,
	                 SubGraphs& subs,
			 vector<vector<pairwise_alignment> >& alignments,
			 double T0, //initial temp
			 double P, //rate of cooling
			 int iterations, //number of sweep to perform
			 double sigma, //scoring paramater
			 double mu, //scoring parameter
			 vector<bool>& s //bitvector of aligned subgraphs
			 )
{
    int cnt = 0; //count number of iterations
    int p = alignments.size(); //number of subgraphs
    int E, Ei, oldEi;
    double Ti = T0; //initial temp
    
    s = vector<bool>(p,false); //if s[i] is true, then the ith subgraph is in the alignment
    
    //random starting point
    random_alignment(s);
    
    oldEi = Ei = E = score_alignment(graph, subs, alignments, s, sigma, mu);	
   
    while(cnt < iterations)
    {
	//our moveset is just to flip a single bit
	int k = rand()%p;

	s[k] = !s[k];

	//evaluate the new alignment
	oldEi = Ei;
	Ei = score_alignment(graph, subs, alignments, s, sigma, mu);

	if(Ei > E)
	{
	    //this is a good move, keep it
	    E = Ei; 
	}
	else //bad move, keep with some probability
	{
	    double e = -(E-Ei)/Ti;
	    e = exp(e);
	    double r = rand();
	    //accept with probability r
	    if(r < e*RAND_MAX) 
	    {
		E = Ei;
	    }
	    else //reject this move
	    {
		//undo
		s[k] = !s[k];
		Ei = oldEi;
	    }

	    Ti = Ti*P;
	}
	cnt++;
    }

    return score_alignment(graph, subs, alignments, s, sigma, mu);
}

//now it would be nice to actually output the motif, which means
//we have to find the actual multiple alignment (what we've found so
//far is actually just the members of the multiple alignment)

//the way we construct this is we just select the alignment of
//a subgraph which is the most popular minimum alignment with all
//the other members of the multiple alignment

void construct_multiple_alignment(const Graph& graph,
	                 SubGraphs& subs,
			 vector<vector<pairwise_alignment> >& alignments,
			 vector<bool>& s,
			 vector<SubGraph>& malign)
{
    malign.clear();

    //first get a list of all the members of the multiple alignment
    vector<uint> members;
    for(int i = 0; i < s.size(); i++)
    {
	if(s[i])
	    members.push_back(i);
    }
    
    malign.reserve(members.size());

    for(int i = 0; i < members.size(); i++)
    {
	//to figure out the best alignment for the ith member
	//pick the alignment that is the most popular
	map<SubGraph, int> vote;
	int bindex = members[i];

	//compute the pairwise alignment to each already chosen
	//alignment in the malign matrix
	for(int j = 0; j < i; j++)
	{
	    SubGraph a = malign[j];
	    
	    //try all possible alignments and choose the minimum alignment
	    int min = INT_MAX;
	    SubGraph b = subs.subgraphs[bindex];
	    SubGraph best_b;

	    //keep permuting b, note that is starts out sorted so we don't miss any permutations
	    do
	    {
		int s = score(graph, a, b);
		if(s < min)
		{
		    best_b = b;
		    min = s; //found a better alignment
		}
	    }
	    while(next_permutation(b.begin(), b.end()));

	    //subgraph a votes for best_b alignment of b
	    vote[best_b]++;
	}

	if(vote.empty()) //first one
	    malign.push_back(subs.subgraphs[bindex]);
	else
	{
	    //find the most popular
	    map<SubGraph, int>::iterator maxitr = vote.end();
	    int maxvote = 0;
	    for(map<SubGraph, int>::iterator itr = vote.begin(); itr != vote.end(); itr++)
	    {
		if((*itr).second > maxvote)
		{
		    maxvote = (*itr).second;
		    maxitr = itr;
		}
	    }
	    assert(maxvote > 0);
	    malign.push_back((*maxitr).first);
	}
    }       
   
}


//create a motif from a multiple alignment
//This is just a graph with edges whose weights are
//between 0 and 1
void create_motif(const Graph& graph, const vector<SubGraph>& malign, vector<vector<double> >& motif)
{
    motif.clear();

    if(malign.size() <= 0)
	return;

    int n = malign[0].size(); //size of the subgraphs

    motif.resize(n);
    for(int i = 0; i < n; i++)
	motif[i] = vector<double>(n, 0.0);

    //first just compute the sum of edges for each subgraph
    for(int g = 0; g < malign.size(); g++)
    {
	for(int i = 0; i < n; i++)
	{
	    for(int j = 0; j < n; j++)
	    {
		if(graph.isEdge(malign[g][i], malign[g][j]))
		    motif[i][j] += 1.0;
	    }
	}
    }

    //then divide by the total number of subgraphs
    for(int i = 0; i < n; i++)
    {
	for(int j = 0; j < n; j++)
	{
	    motif[i][j] /= malign.size();
	}
    }
    //and that's it :-)
}


void myerror(char *msg)
{
    fprintf(stderr, "ERROR: %s\n", msg);
    exit(-1);
}


// print out usage message and exit
void usage(char * progname)
{
    fprintf(stderr,"Usage: %s [-v] [-s] [-t <threshold>] [-mu <mu>] [-sigma <sigma>] -f file -k <subgraph size>\n",progname);
    exit(-1);
}



int main(int argc, char *argv[])
{
    unsigned int i, j, n;
    ifstream gfile;
    string tmpstr, name, file_info;
    int k = 0;
    double mu = 5.0;
    double sigma = 10.0;
    int sim_iter = 5000; //number of simulated annealing iterations
    int threshold = -1; //use the threshold most significant subgraphs
    enum times {start, build_graph, find_subgraphs, pairwise_alignments, sim_anneal, multiple_align, build_motif, end};
    clock_t time_arr[end+1];
    
    /* process cmdline arguments */

    for(i = 1; i < argc; i++)
    {
	if(strcmp(argv[i],"-v") == 0)
	{
	    verbose = true;
	    summarize = true;
	}
	else if(strcmp(argv[i],"-f") == 0)
	{
	    i++;
	    if(i == argc)
		usage(argv[0]);
	    gfile.open(argv[i]);
	    if(!gfile)
		myerror("Could not open file.");	    
	}
	else if(strcmp(argv[i],"-k") == 0)
	{
	    i++;
	    if(i == argc)
		usage(argv[0]);
	    k = atoi(argv[i]);		
	}
	else if(strcmp(argv[i],"-s") == 0)
	{
	    summarize = true;
	}
	else if(strcmp(argv[i],"-t") == 0)
	{
	    i++;
	    threshold = atoi(argv[i]);
	    if(threshold <= 0)
		myerror("Threshold must be positive.");
	}
	else if(strcmp(argv[i],"-mu") == 0)
	{
	    i++;
	    mu = atof(argv[i]);
	    if(mu == 0)
		myerror("mu can't be zero. Or at least I don't think it can be..");
	}
	else if(strcmp(argv[i],"-sigma") == 0)
	{
	    i++;
	    sigma = atof(argv[i]);
	    if(sigma == 0)
		myerror("sigma can't be zero. Or at least I don't think it can be..");
	}
	else if(strcmp(argv[i],"-simiter") == 0)
	{
	    i++;
	    sim_iter = atoi(argv[i]);
	    if(sim_iter <= 0)
		myerror("simiter can't be zero");
	}
	else if(strcmp(argv[i],"-prune") == 0)
	{
	    prune = true;
	}
    }

    if(!gfile.is_open())
	usage(argv[0]);

    if(k == 0)
	usage(argv[0]);


    time_arr[start] = clock();
    
    /* load the graph */
    if(!(gfile >> n))
	myerror("Couldn't read number of nodes.");

    getline(gfile, tmpstr); //absorbe newline
    getline(gfile, file_info); //file info
    getline(gfile, tmpstr); //absorb newline

    n++; //apparently we're indexing by 1..
    
    Graph graph(n, file_info);

    /* read in the nodes and edges.
       node name, which is the gene name in the original Ito data
       node information, which is the (Bait/Prey) description in Ito data
       node number, which is a number assigned to the node uniquely 
       neighboring node numbers, space separated
    */
    while(getline(gfile, name))  //node name
    {
	string info;
	getline(gfile, info);

	gfile >> i; //node number	
	gfile.ignore(INT_MAX,'\n'); //absorbe newline
	if(i >= n) myerror("Node out of range.");

	graph.setNode(i,name,info);

	//process edges
	getline(gfile, tmpstr);
	istringstream ist(tmpstr);
	while(ist >> j)
	{
	    graph.addEdge(i,j);
	    graph.addEdge(j,i); //this forces the graph to be undirected..
	}
	
	gfile.ignore(INT_MAX,'\n'); //absorbe newline
    }

    time_arr[build_graph] = clock();
    if(summarize)
    {
	cout << "Number of nodes: " << graph.numNodes() << endl;
    }

    //as a sanity check, count number of unnamed graph nodes
    int cnt = 0;
    for(int i = 0; i < graph.numNodes(); i++)
    {
	if(graph.getName(i) == "")
	    cnt++;
    }

    if(cnt)
	cout << cnt << " Unnamed Nodes\n";

    fflush(stdout);
    SubGraphs subs;    
    int subs_cnt = graph.subgraphs(subs, k, threshold);

    time_arr[find_subgraphs] = clock();

    if(summarize)
    {
	cout << "Number of subgraphs: " << subs_cnt << endl;
    }

    if(threshold > 0 && threshold < subs.subgraphs.size()) //prune the number of subgraphs
    {
	SubGraphs newsubs;
	newsubs.n = threshold;
	
	//select the threshold many subgraphs with the lowest null probability
	//(those that are unlikely to occur by chance will probably be more meaningful)
	for(int i = 0; i < threshold; i++)
	{
	    double min = 1.0;
	    int min_loc = -1;
	    //find the least likely subgraph
	   
	    for(int g = 0; g < subs.n; g++)
	    {
		if(subs.P0[g] < min)
		{
		    min_loc = g;
		    min = subs.P0[g];
		}
	    }
	    assert(min_loc >= 0);

	    newsubs.subgraphs.push_back(subs.subgraphs[min_loc ]);
	    newsubs.P0.push_back(subs.P0[min_loc ]);
	    newsubs.L.push_back(subs.L[min_loc ]);
	    //remove this subgraph from consideration
	    subs.P0[min_loc ] = 2.0;
	}

	assert(newsubs.subgraphs.size() == newsubs.n);
	subs = newsubs;
    }
    
    if(verbose)
    {
	cout << "Subgraphs\n---------" << endl;
	for(int i = 0; i < subs.subgraphs.size(); i++)
	{
	    cout << i << ":\t";
	    for(int j = 0; j < subs.subgraphs[i].size(); j++)
	    {
		cout << subs.subgraphs[i][j] << "(" << graph.getName(subs.subgraphs[i][j]) << ") ";
	    }
	    cout << "\t" << subs.P0[i];
	    cout << "\t" << subs.L[i];
	    cout << endl;
	}
    }
    

    fflush(stdout);
    //now find the pairwise minimum mismatch of all pairs of subgraphs
    vector< vector<pairwise_alignment> > mismatches;
    compute_pairwise_mismatches(graph, subs, mismatches);

    time_arr[pairwise_alignments] = clock();
    
    if(summarize) //create and print histogram
    {
	vector<int> cnt;
	for(int i = 0; i < mismatches.size(); i++)
	{
	    for(int j = 0; j < i; j++)
	    {
		int s = mismatches[i][j].score;
		if(s >= cnt.size())
		    cnt.resize(s+1);
		cnt[s]++;
	    }
	}

	cout << "Minimum mismatch histogram" << endl;
	for(int i = 0; i < cnt.size(); i++)
	{
	    cout << i << ":\t" << cnt[i] << endl;
	}	
    }
    
    if(verbose && 0)
    {
	cout << "\nMinimal Mismatches\n------------------" << endl;

	for(int i = 0; i < mismatches.size(); i++)
	{
	    for(int j = 0; j < i; j++)
	    {
		cout << mismatches[i][j].score << " ";
	    }
	    cout << endl;
	}
    }


    fflush(stdout);
    
    vector<bool> s; //bitvector of subgraphs in multiple alignment
    double score = simulated_annealing(graph, subs, mismatches,
			 20.0, //initial temp
			 .97, //rate of cooling
			 sim_iter, //number of sweep to perform
			 sigma, //sigma scoring paramater -> more internal links
			 mu, //mu scoring parameter -> less fuzziness
			 s //bitvector of aligned subgraphs
	    );

    time_arr[sim_anneal] = clock();

    time_arr[end] = clock();
    if(verbose || summarize)
    {
	cout << "Score of found maximum multiple alignment: " << score << endl;
	cout << "Members of multiple alignment:\n";
	int cnt = 0;
	for(int i = 0; i < s.size(); i++)
	{
	    if(s[i])
	    {
		cnt++;
		if(verbose)
		    cout << i << " ";
	    }
	}
	cout << endl;
	cout << "Size of multiple alignment: " << cnt << endl;
    }

    //construct actual motif
    vector<SubGraph> malign;
    construct_multiple_alignment(graph, subs, mismatches, s, malign);
    
    time_arr[multiple_align] = clock();

    if(verbose)
    {
	cout << "Multiple alignment:\n";
	for(int i = 0; i < malign.size(); i++)
	{
	    for(int j = 0; j < malign[i].size(); j++)
	    {
		cout << malign[i][j] << " " << graph.getName(malign[i][j]) << "\t";
	    }
	    cout << endl;
	}
    }
    
    //from the multiple alignment, get the motif
    vector<vector<double> > motif;
    create_motif(graph, malign, motif);

    if(summarize || verbose)
    {
	//print out motif
	cout << "Motif:" << endl;
	for(int i = 0; i < motif.size(); i++)
	{
	    for(int j = 0; j < motif[i].size(); j++)
	    {
		printf("%1.4f\t",motif[i][j]);
	    }
	    cout << endl;
	}
    }
    
    time_arr[build_motif] = clock();
    
    //output timing information
    if(verbose)
    {
	cout << "--------"<< endl;
	cout << "Building Graph........\t" << (time_arr[build_graph]-time_arr[start])/(double)CLOCKS_PER_SEC << endl;
	cout << "Finding Subgraphs.....\t" << (time_arr[find_subgraphs]-time_arr[build_graph])/(double)CLOCKS_PER_SEC << endl;
	cout << "Pairwise Alignments...\t" << (time_arr[pairwise_alignments]-time_arr[find_subgraphs])/(double)CLOCKS_PER_SEC << endl;
	cout << "Simulated Annealing...\t" << (time_arr[sim_anneal]-time_arr[pairwise_alignments])/(double)CLOCKS_PER_SEC << endl;
	cout << "Multiple Alignment....\t" << (time_arr[multiple_align]-time_arr[sim_anneal])/(double)CLOCKS_PER_SEC << endl;
	cout << "Build Motif...........\t" << (time_arr[build_motif]-time_arr[multiple_align])/(double)CLOCKS_PER_SEC << endl;
	cout << "Total.................\t" << (time_arr[end]-time_arr[start])/(double)CLOCKS_PER_SEC << endl;
    }
    
    return 0;
}