/*
 * This program takes two input sequences, one from human and one from mouse,
 * and calculates a conservation score for them.
 */

#include <stdio.h>
#include <strings.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include "matrix.h"
#include "iterator.h"
#include "probability.h"
#include "const.h"
#include "distrib.h"
#include <iostream>

extern int runmode;
int rebuilding=0;
extern Matrix<double> *m24ihash[];
extern Row<double> f3c;

extern int mfcount;
extern int m24count;
extern int memuse;
extern int footprint;
extern int computing;
extern int gapless;
extern int cposition;
extern int recompute;
extern double interval;
using namespace std;
void clearhash();
int gsize=GLET;
#define IN(a,b,c) (a+b*GLET+c*GLET*GLET)
double gval_ml(Iterator *i, int i1);
double compute(Iterator *i);

void shift2(double &a, double &b, const double c);
void shift3(double &a, double &b, double &c, const double d);

#define G_LEFT 0
#define G_MID 1
#define G_RIGHT 2

extern int use_bayes;
int gmode = G_LEFT;

double gval_one(Iterator *i, int i1, rpfunc f);

double gs_one(int i1, Iterator *i, rpfunc f)
{
	if (!i) return 1;

	double r = 1e20;

	// you're a leaf
	if (i->is_leaf()) {
		if (gmode == G_LEFT) {
			i->record_state(i->b1,0,0, 0,0);
			r = subst1( i1, i->b1, r);
		} else if (gmode == G_RIGHT) {
			i->record_state(i->b3,0,0, 0,0);
			r = subst1( i1,i->b3, r);
		} else {
			i->record_state(i->b2,0,0, 0,0);
			r = subst1( i1,i->b2, r);
		}
		return r;
	}

	// you're not a leaf - loop through possible states and find one that minimizes tree
	int j;
	for (j=0;j<gsize;j++) {
		// skip non-legal states
		if (gmode != G_MID && i->p2 >= G1 && j >= G1 && j != i->p2) continue;
		double s = 0;
		s = gval_one(i, j, f);
		s += subst1( i1, j, r);
		if (s <= r) {
			if (r == s) {
				// collision!
				int temp = i->gn[i1];
				
				i->gn[i1] = j;
				i->record_state(j,0,0,s,0);
				i->make_record();
				double lik = i->record_likelihood(f);
				
				if (lik > i->gl[i1]) {
					// take the more likely scenario.
					i->gl[i1] = lik;
				} else {
					i->gn[i1] = temp;
					i->record_state(temp,0,0, s,0);
					i->make_record();
				}
			} else {
				i->gn[i1] = j;
				i->record_state(j,0,0,s,0);
				i->make_record();
				i->gl[i1] = i->record_likelihood(f);
			}
			r = s;
		}
	}
	return r;
}

double gval_one(Iterator *i, int i1, rpfunc f)
{
	if (i->gm[i1]) return i->gv[i1];

	double a,b;
	
	a = gs_one(i1,i->right,f);
	b = gs_one(i1,i->left,f);

	i->gv[i1] = a+b;
	i->gm[i1] = 1;
	return i->gv[i1];
}

double compute_one(Iterator *i, rpfunc f) {
	double val = high();
	int low = -1;
	Stack *b = new Stack;
	int i1,i2,i3, j;
	for (j=0;j<gsize;j++) {
		double s = gval_one(i,j,f);
		if (s <= val) {
			if (s == val) {
				// always prefer insertions in the root, so we can reparent.
				if  (i->record->one < G1 || j >= G1) {
					Stack temp = *(i->record);
					i->record_state(j,0,0,s,0);
					i->make_record();
					double l = i->record_likelihood(f);
					if (l > i->gl[temp.one]) {
						i->gl[j] = l;
						low = j;
					} else {
						i->record_state(temp.one,0,0,s,0);
						low = temp.one;
					}
				}
			} else {
				i->record_state(j,0,0,s,0);
				i->make_record();
				i->gl[j] = i->record_likelihood(f);
				low = j;
			}
			val = s;
		}
	}
	i->record_state(low,0,0,0,0);
	i->make_record();
 	
	return val;
}

void label_tree(Iterator *i)
{
	i->reset();
	gmode = G_MID;
	compute_one(i,single_prob);
	i->flash_record(G_MID);
}

void print_sis(Iterator *i) {
	if (i->sis) {
		if (i->record->one >= G1)
		cout << '-';
		else
		cout << let[i->record->one];
		return;
	}
	if (i->right) print_sis(i->right);
	if (i->left) print_sis(i->left);
}

void run_states_anc(Iterator *i)
{
	// get most parsimonious tree - how?
	int j,k,l,m,n;
	label_tree(i);

	print_sis(i);
}