
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "nrutil.h"
#include <iostream.h>
#include <fstream.h>
#include <iomanip.h>
#include <time.h>
#include <string.h>
#include <ctype.h>
#include <malloc.h>

//#ifndef __ANCESTRAL_H
//#define __ANCESTRAL_H

//#define MAXSTR 500
//#define MAXNODE 6000
//#define MAXLENGTH 1000
#define NP 20
#define E 2.718281828
//for readali
#define MAXSTR   10001
#define NR_END 1
#define FTOL 1.0e-6

struct NODE {
	char nodeName[50];
	int nodeNum;
	int nodeLevel;
	int childL;
	int childR;
	bool rootFlag;
	int childE;
	int parent;
	double dist;
	int seqIndex;
	char *seq_aa;//[MAXLENGTH];
	int *seq_num;//[MAXLENGTH];
	int C[NP+1];
	double L[NP+1];
	double amp;
	double **prob;//[MAXLENGTH][NP+1];
	double tmpL[NP+1][NP+1];//used for goldstein method
};

struct topology{
	int nodeLevel;
	int childL;
	int childR;
	int childE;
	int parent;
	bool rootFlag;
	double dist;
	bool gapFlag;
	double L[NP+1];
	double amp;
};

void usage();
void readTree(char treeFile[]);
void init_S();
void searchGap();
void cal_PI();
void cal_Q(double pi[]);
void cal_Eign();
void cal_dist();
double dist(int node1, int node2);
char *itoa(int b);
double functionValue(int **Nab, double d);
double dpij(int s1,int s2,double d, double alfa);
double pij(int s1,int s2,double d, double alfa);
void optimalAlfa(int refineFlag,int condProb);
void noOptimalAlfa(int condProb);
double calTree(int siteNo,int root, double alfa);
void optTree(int siteNo, double alfa);
void optTree_Goldstein(int siteNo, double alfa);
void assignTree();
int rebuildTree(int siteNo);
void print_tree(int nodeNo);
topology **topmatrix(long nrl, long nrh, long ncl, long nch);
void setLevel(int nodeNum);
double functionValue1(int **Nab, double d);
void setRLevel(int siteNo,int nodeNum);
void calRL(int siteNo,int nodeNo,double alfa);
void weighbor(int N, int sequenceLength, NODE *node, double **distMatrix);
void changeRoot(int siteNo, int nodeNo, int root, int child);


void balanc(double **a, int n);
void elmhes(double **a, int n);
void hqr(double **a, int n, double wr[], double wi[]);
void svdcmp(double **a, int m, int n, double w[], double **v);
void ludcmp(double **a, int n, int *indx, double *d);
void lubksb(double **a, int n, int *indx, double b[]);

//for readali
static char *mymalloc(int size);
char *ccvector(long nl, long nh);
char **cmatrix(long nrl, long nrh, long ncl, long nch);
//void nrerror(char error_text[]);
char *strsave(char *str);
char *strnsave(char *str, int l);
static char **incbuf(int n, char **was);
static int *incibuf(int n, int *was);

void readali(char *filename);


//#endif
