#ifndef ITERATOR_H
#define ITERATOR_H
#include <stdlib.h>
#include <string>
#include <vector>
#include "matrix.h"
#include "probability.h"

#define LET 4
#define GLET 8
#define G0 3
#define G1 4
#define G2 5
#define G3 6
#define GN 7

#define IN(a,b,c) (a+b*GLET+c*GLET*GLET)
#define IN2(a,b,c) (a+b*LET+c*LET*LET)
#define MI(a,b) (a+b*GLET)
#define FL(a,b,c,d) (a+b*GLET+c*GLET*GLET+d*GLET*GLET*GLET)

class Stack
{
	public:
	Stack() {next=0;}
	Stack *next;
	int one;
	int two;
	int thr;
};

class TreeP
{
	public:
		TreeP() {prob=1;count=0;}
		int count;
		double prob;
		TreeP operator += (TreeP s) {
			prob *= s.prob;
			count += s.count;
			return *this;
		}
};

class Iterator
{
	public:
	Iterator();
	Iterator(char *t, Iterator *p, double scale);
	~Iterator();
	Iterator *right;
	Iterator *left;
	Iterator *next;
	Iterator *parent;
	Iterator *tail;
	Stack *record;
	double len;
	double scale;
	double actual_time;
	double tree_len();
	double longest_branch();
	double distance(std::string s1, std::string s2);
	double find_node(std::string s);
	double prob, rprob, nprob;
	short int *gm;
	double *gv;
	double *gl;
	int *gn;
	std::string species;
	int val;
	int ref;
	int sis;
	double mod;
	int is_leaf() { return (!right && !left); }
	Iterator *leaf_prev();
	int count_subs(int record);
	int count_nodes();
	int count_empty();
	int count_gaps(int in);
	double cons_score(int l);
// 	TreeP tree_score(int black, int center);
	TreeP tree_score(int black, int center, double (*f)(Iterator*, Stack*) );
	void set_base_state(int a, int b, int c);
	void set_species(std::string a) {species = a;}
	void set_actual_time(double time) { actual_time = time; }
	void reset();
	void flash_record(int l);
	void record_from_flash();

	int flag(std::string s);
	void reset_flags();
	void linked_list(Iterator **c);
	Iterator *flag_tree(int f, double l, pfunc p);
	int tflag;
	int flag_ref(std::string s, int unflag=1);
	int flag_sis(std::string s);
	
	void rebuild(double s, pfunc p);
	int subs_in(std::string spec);
	Iterator *clone_tree();

	int set_nuc(std::string species, std::string nuc);

	double substitute(int n1, int n2, int let);
	double record_likelihood( rpfunc f );
	double tree_likelihood( rpfunc f );
	void record_state(int a, int b, int c, double p, double n);
	void make_record();
	void print_state(int depth=0);
	void print_root();
	void print_nucs();
	void print_tree();
	void print_rec_nucs();
	void get_nucs(std::string s);
	double golden_max(std::vector<double> t, std::vector<double> s);
	double rate_score();
	void rate_branch_score(double &c, std::vector<double> &mr1, std::vector<double> &mr2);
	std::string key();
	std::vector<std::string> keystring();
	std::string keycache;
	double ins(int i);
	double del(int d, int isize=1);
	Matrix<double> *m1, *m2, *m3, *m24, *m33, *m33c, *m31;
	int b1, b2, b3;
	int g1, g2, g3;
	int p1, p2, p3;
	int ins_state;
	void reflag();
	Row<double> insertion, deletion;
	void build_matrices(double time, pfunc p);
};

void load_matrices(std::string, std::string, std::string, std::string);

extern Row<double> f0;
extern Row<double> f1;
extern Row<double> f2;
extern Row<double> f3;
#endif    // ITERATOR_H