Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

source: branches/nameserver/GDE/PHYLIP/promlk.c

Visit:

Last change on this file was 2175, checked in by westram, 21 years ago
upgrade to PHYLIP 3.6
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 83.7 KB

Line
1
2	#include "phylip.h"
3	#include "seq.h"
4
5	/* version 3.6. (c) Copyright 1986-2000 by the University of Washington
6	and by Joseph Felsenstein. Written by Joseph Felsenstein and Lucas Mix.
7	Permission is granted to copy and use this program provided no fee is
8	charged for it and provided that this copyright notice is not removed. */
9
10	#define epsilon 0.0001 /* used in makenewv, getthree, update */
11	#define over 60
12
13	typedef long vall[maxcategs];
14	typedef double contribarr[maxcategs];
15
16	#ifndef OLDC
17	/* function prototypes */
18	void init_protmats(void);
19	void getoptions(void);
20	void makeprotfreqs(void);
21	void allocrest(void);
22	void doinit(void);
23	void inputoptions(void);
24	void input_protdata(long);
25	void makeweights(void);
26	void prot_makevalues(long, pointarray, long, long, sequence, steptr);
27	void getinput(void);
28
29	void prot_inittable(void);
30	void alloc_pmatrix(long);
31	void make_pmatrix(double , double , double **, long, double, double,
32	double , double *);
33	void prot_nuview(node *);
34	void getthree(node *p, double thigh, double tlow);
35	void makenewv(node *);
36	void update(node *);
37	void smooth(node *);
38	void promlk_add(node , node , node *, boolean);
39	void promlk_re_move(node , node , boolean);
40
41	double prot_evaluate(node *);
42	void tryadd(node , node , node *);
43	void addpreorder(node , node , node *, boolean, boolean);
44	void restoradd(node , node , node *, double);
45	void tryrearr(node , boolean );
46	void repreorder(node , boolean );
47	void rearrange(node **);
48	void nodeinit(node *);
49	void initrav(node *);
50	void travinit(node *);
51
52	void travsp(node *);
53	void treevaluate(void);
54	void promlk_coordinates(node , long );
55	void promlk_drawline(long, double);
56	void promlk_printree(void);
57	void describe(node *);
58	void prot_reconstr(node *, long);
59	void rectrav(node *, long, long);
60	void summarize(void);
61	void promlk_treeout(node *);
62	void initpromlnode(node , node , node , long, long, long , long *,
63	initops, pointarray, pointarray, Char , Char , FILE *);
64	void tymetrav(node , double );
65
66	void free_all_protx(long, pointarray);
67	void maketree(void);
68	void clean_up(void);
69	/* function prototypes */
70	#endif
71
72
73	extern sequence y;
74	double **tbl;
75
76	Char infilename[100], outfilename[100], intreename[100], outtreename[100],
77	catfilename[100], weightfilename[100];
78	double *rrate;
79	long sites, weightsum, categs, datasets, ith, njumble, jumb, numtrees;
80	/* sites = number of sites in actual sequences
81	numtrees = number of user-defined trees */
82	long inseed, inseed0, mx, mx0, mx1;
83	boolean global, jumble, lngths, trout, usertree, weights, rctgry, ctgry,
84	auto_, progress, mulsets, firstset, hypstate, reconsider,
85	smoothit, polishing, justwts, gama, invar, usepam, usejtt;
86	tree curtree, bestree, bestree2;
87	node qwhere, grbg;
88	double sumrates, cv, alpha, lambda, lambda1, invarfrac;
89	long *enterorder;
90	steptr aliasweight;
91	double *rate;
92	longer seed;
93	double *probcat;
94	contribarr *contribution;
95	char aachar[26]="ARNDCQEGHILKMFPSTWYVBZX?*-";
96	char *progname;
97	long rcategs, nonodes2;
98
99
100	/* Local variables for maketree, propagated globally for C version: */
101	long k, maxwhich, col;
102	double like, bestyet, tdelta, lnlike, slope, curv, maxlogl;
103	boolean lastsp, smoothed, succeeded;
104	double *l0gl;
105	double x[3], lnl[3];
106	double expon1i[maxcategs], expon1v[maxcategs],
107	expon2i[maxcategs], expon2v[maxcategs];
108	node *there;
109	double **l0gf;
110	Char ch, ch2;
111	vall *mp;
112
113
114	/* Variables introduced to allow for protein probability calculations */
115	long max_num_sibs; /* maximum number of siblings used in a */
116	/* nuview calculation. determines size */
117	/* final size of pmatrices */
118	double eigmat; / eig matrix variable */
119	double *probmat; / prob matrix variable */
120	double ***dpmatrix; / derivative of pmatrix */
121	double ***ddpmatrix; / derivative of xpmatrix */
122	double ****pmatrices; / matrix of probabilities of protien */
123	/* conversion. The 5 subscripts refer */
124	/* to sibs, rcategs, categs, final and */
125	/* initial states, respectively. */
126	double freqaa[20]; /* amino acid frequencies */
127
128	static double pameigmat[20] =
129	{
130	-0.022091252,-0.019297602, 0.000004760,-0.017477817,
131	-0.016575549,-0.015504543,-0.002112213,-0.002685727,
132	-0.002976402,-0.013440755,-0.012926992,-0.004293227,
133	-0.005356688,-0.011064786,-0.010480731,-0.008760449,
134	-0.007142318,-0.007381851,-0.007806557,-0.008127024
135	};
136
137	static double pamprobmat[20][20] =
138	{
139	{-0.01522976,-0.00746819,-0.13934468, 0.11755315,-0.00212101,
140	0.01558456,-0.07408235,-0.00322387, 0.01375826, 0.00448826,
141	0.00154174, 0.02013313,-0.00159183,-0.00069275,-0.00399898,
142	0.08414055,-0.01188178,-0.00029870, 0.00220371, 0.00042546},
143	{-0.07765582,-0.00712634,-0.03683209,-0.08065755,-0.00462872,
144	-0.03791039, 0.10642147,-0.00912185, 0.01436308,-0.00133243,
145	0.00166346, 0.00624657,-0.00003363,-0.00128729,-0.00690319,
146	0.17442028,-0.05373336,-0.00078751,-0.00038151, 0.01382718},
147	{-0.08810973,-0.04081786,-0.04066004,-0.04736004,-0.03275406,
148	-0.03761164,-0.05047487,-0.09086213,-0.03269598,-0.03558015,
149	-0.08407966,-0.07970977,-0.01504743,-0.04011920,-0.05182232,
150	-0.07026991,-0.05846931,-0.01016998,-0.03047472,-0.06280511},
151	{ 0.02513756,-0.00578333, 0.09865453, 0.01322314,-0.00310665,
152	0.05880899,-0.09252443,-0.02986539,-0.03127460, 0.01007539,
153	-0.00360119,-0.01995024, 0.00094940,-0.00145868,-0.01388816,
154	0.11358341,-0.12127513,-0.00054696,-0.00055627, 0.00417284},
155	{ 0.16517316,-0.00254742,-0.03318745,-0.01984173, 0.00031890,
156	-0.02817810, 0.02661678,-0.01761215, 0.01665112, 0.10513343,
157	-0.00545026, 0.01827470,-0.00207616,-0.00763758,-0.01322808,
158	-0.02202576,-0.07434204, 0.00020593, 0.00119979,-0.10827873},
159	{ 0.16088826, 0.00056313,-0.02579303,-0.00319655, 0.00037228,
160	-0.03193150, 0.01655305,-0.03028640, 0.01367746,-0.11248153,
161	0.00778371, 0.02675579, 0.00243718, 0.00895470,-0.01729803,
162	-0.02686964,-0.08262584, 0.00011794,-0.00225134, 0.09415650},
163	{-0.01739295, 0.00572017,-0.00712592,-0.01100922,-0.00870113,
164	-0.00663461,-0.01153857,-0.02248432,-0.00382264,-0.00358612,
165	-0.00139345,-0.00971460,-0.00133312, 0.01927783,-0.01053838,
166	-0.00911362,-0.01010908, 0.09417598, 0.01763850,-0.00955454},
167	{ 0.01728888, 0.01344211, 0.01200836, 0.01857259,-0.17088517,
168	0.01457592, 0.01997839, 0.02844884, 0.00839403, 0.00196862,
169	0.01391984, 0.03270465, 0.00347173,-0.01940984, 0.01233979,
170	0.00542887, 0.01008836, 0.00126491,-0.02863042, 0.00449764},
171	{-0.02881366,-0.02184155,-0.01566086,-0.02593764,-0.04050907,
172	-0.01539603,-0.02576729,-0.05089606,-0.00597430, 0.02181643,
173	0.09835597,-0.04040940, 0.00873512, 0.12139434,-0.02427882,
174	-0.02945238,-0.01566867,-0.01606503, 0.09475319, 0.02238670},
175	{ 0.04080274,-0.02869626,-0.05191093,-0.08435843, 0.00021141,
176	0.13043842, 0.00871530, 0.00496058,-0.02797641,-0.00636933,
177	0.02243277, 0.03640362,-0.05735517, 0.00196918,-0.02218934,
178	-0.00608972, 0.02872922, 0.00047619, 0.00151285, 0.00883489},
179	{-0.02623824, 0.00331152, 0.03640692, 0.04260231,-0.00038223,
180	-0.07480340,-0.01022492,-0.00426473, 0.01448116, 0.01456847,
181	0.05786680, 0.03368691,-0.10126924,-0.00147454, 0.01275395,
182	0.00017574,-0.01585206,-0.00015767,-0.00231848, 0.02310137},
183	{-0.00846258,-0.01508106,-0.01967505,-0.02772004, 0.01248253,
184	-0.01331243,-0.02569382,-0.04461524,-0.02207075, 0.04663443,
185	0.19347923,-0.02745691, 0.02288515,-0.04883849,-0.01084597,
186	-0.01947187,-0.00081675, 0.00516540,-0.07815919, 0.08035585},
187	{-0.06553111, 0.09756831, 0.00524326,-0.00885098, 0.00756653,
188	0.02783099,-0.00427042,-0.16680359, 0.03951331,-0.00490540,
189	0.01719610, 0.15018204, 0.00882722,-0.00423197,-0.01919217,
190	-0.02963619,-0.01831342,-0.00524338, 0.00011379,-0.02566864},
191	{-0.07494341,-0.11348850, 0.00241343,-0.00803016, 0.00492438,
192	0.00711909,-0.00829147, 0.05793337, 0.02734209, 0.02059759,
193	-0.02770280, 0.14128338, 0.01532479, 0.00364307, 0.05968116,
194	-0.06497960,-0.08113941, 0.00319445,-0.00104222, 0.03553497},
195	{ 0.05948223,-0.08959930, 0.03269977,-0.03272374,-0.00365667,
196	-0.03423294,-0.06418925,-0.05902138, 0.05746317,-0.02580596,
197	0.01259572, 0.05848832, 0.00672666, 0.00233355,-0.05145149,
198	0.07348503, 0.11427955, 0.00142592,-0.01030651,-0.04862799},
199	{-0.01606880, 0.05200845,-0.01212967,-0.06824429,-0.00234304,
200	0.01094203,-0.07375538, 0.08808629, 0.12394822, 0.02231351,
201	-0.03608265,-0.06978045,-0.00618360, 0.00274747,-0.01921876,
202	-0.01541969,-0.02223856,-0.00107603,-0.01251777, 0.05412534},
203	{ 0.01688843, 0.05784728,-0.02256966,-0.07072251,-0.00422551,
204	-0.06261233,-0.08502830, 0.08925346,-0.08529597, 0.01519343,
205	-0.05008258, 0.10931873, 0.00521033, 0.02593305,-0.00717855,
206	0.02291527, 0.02527388,-0.00266188,-0.00871160, 0.02708135},
207	{-0.04233344, 0.00076379, 0.01571257, 0.04003092, 0.00901468,
208	0.00670577, 0.03459487, 0.12420216,-0.00067366,-0.01515094,
209	0.05306642, 0.04338407, 0.00511287, 0.01036639,-0.17867462,
210	-0.02289440,-0.03213205, 0.00017924,-0.01187362,-0.03933874},
211	{ 0.01284817,-0.01685622, 0.00724363, 0.01687952,-0.00882070,
212	-0.00555957, 0.01676246,-0.05560456,-0.00966893, 0.06197684,
213	-0.09058758, 0.00880607, 0.00108629,-0.08308956,-0.08056832,
214	-0.00413297, 0.02973107, 0.00092948, 0.07010111, 0.13007418},
215	{ 0.00700223,-0.01347574, 0.00691332, 0.03122905, 0.00310308,
216	0.00946862, 0.03455040,-0.06712536,-0.00304506, 0.04267941,
217	-0.10422292,-0.01127831,-0.00549798, 0.11680505,-0.03352701,
218	-0.00084536, 0.01631369, 0.00095063,-0.09570217, 0.06480321}
219	};
220
221	/* this jtt matrix decomposition due to Elisabeth Tillier */
222	static double jtteigmat[] =
223	{0.0, -0.007031123, -0.006484345, -0.006086499, -0.005514432,
224	-0.00772664, -0.008643413, -0.010620756, -0.009965552, -0.011671808,
225	-0.012222418,-0.004589201, -0.013103714, -0.014048038, -0.003170582,
226	-0.00347935, -0.015311677, -0.016021194, -0.017991454, -0.018911888};
227
228	static double jttprobmat[20][20] =
229	{{0.076999996, 0.051000003, 0.043000004, 0.051999998, 0.019999996, 0.041,
230	0.061999994, 0.073999997, 0.022999999, 0.052000004, 0.090999997, 0.058999988,
231	0.024000007, 0.04, 0.050999992, 0.069, 0.059000006, 0.014000008, 0.032000004,
232	0.066000005},
233	{0.015604455, -0.068062363, 0.020106264, 0.070723273, 0.011702977, 0.009674053,
234	0.074000798, -0.169750458, 0.005560808, -0.008208636, -0.012305869,
235	-0.063730179, -0.005674643, -0.02116828, 0.104586169, 0.016480839, 0.016765139,
236	0.005936994, 0.006046367, -0.0082877},
237	{-0.049778281, -0.007118197, 0.003801272, 0.070749616, 0.047506147,
238	0.006447017, 0.090522425, -0.053620432, -0.008508175, 0.037170603,
239	0.051805545, 0.015413608, 0.019939916, -0.008431976, -0.143511376,
240	-0.052486072, -0.032116542, -0.000860626, -0.02535993, 0.03843545},
241	{-0.028906423, 0.092952047, -0.009615343, -0.067870117, 0.031970392,
242	0.048338335, -0.054396304, -0.135916654, 0.017780083, 0.000129242,
243	0.031267424, 0.116333586, 0.007499746, -0.032153596, 0.033517051,
244	-0.013719269, -0.00347293, -0.003291821, -0.02158326, -0.008862168},
245	{0.037181176, -0.023106564, -0.004482225, -0.029899635, 0.118139633,
246	-0.032298569, -0.04683198, 0.05566988, -0.012622847, 0.002023096,
247	-0.043921088, -0.04792557, -0.003452711, -0.037744513, 0.020822974,
248	0.036580187, 0.02331425, -0.004807711, -0.017504496, 0.01086673},
249	{0.044754061, -0.002503471, 0.019452517, -0.015611487, -0.02152807,
250	-0.013131425, -0.03465365, -0.047928912, 0.020608851, 0.067843095,
251	-0.122130014, 0.002521499, 0.013021646, -0.082891087, -0.061590119,
252	0.016270856, 0.051468938, 0.002079063, 0.081019713, 0.082927944},
253	{0.058917882, 0.007320741, 0.025278141, 0.000357541, -0.002831285,
254	-0.032453034, -0.010177288, -0.069447924, -0.034467324, 0.011422358,
255	-0.128478324, 0.04309667, -0.015319944, 0.113302422, -0.035052393,
256	0.046885372, 0.06185183, 0.00175743, -0.06224497, 0.020282093},
257	{-0.014562092, 0.022522921, -0.007094389, 0.03480089, -0.000326144,
258	-0.124039037, 0.020577906, -0.005056454, -0.081841576, -0.004381786,
259	0.030826152, 0.091261631, 0.008878828, -0.02829487, 0.042718836,
260	-0.011180886, -0.012719227, -0.000753926, 0.048062375, -0.009399129},
261	{0.033789571, -0.013512235, 0.088010984, 0.017580292, -0.006608005,
262	-0.037836971, -0.061344686, -0.034268357, 0.018190209, -0.068484614,
263	0.120024744, -0.00319321, -0.001349477, -0.03000546, -0.073063759,
264	0.081912399, 0.0635245, 0.000197, -0.002481798, -0.09108114},
265	{-0.113947615, 0.019230545, 0.088819683, 0.064832765, 0.001801467,
266	-0.063829682, -0.072001633, 0.018429333, 0.057465965, 0.043901014,
267	-0.048050874, -0.001705918, 0.022637173, 0.017404665, 0.043877902,
268	-0.017089594, -0.058489485, 0.000127498, -0.029357194, 0.025943972},
269	{0.01512923, 0.023603725, 0.006681954, 0.012360216, -0.000181447,
270	-0.023011838, -0.008960024, -0.008533239, 0.012569835, 0.03216118,
271	0.061986403, -0.001919083, -0.1400832, -0.010669741, -0.003919454,
272	-0.003707024, -0.026806029, -0.000611603, -0.001402648, 0.065312824},
273	{-0.036405351, 0.020816769, 0.011408213, 0.019787053, 0.038897829,
274	0.017641789, 0.020858533, -0.006067252, 0.028617353, -0.064259496,
275	-0.081676567, 0.024421823, -0.028751676, 0.07095096, -0.024199434,
276	-0.007513119, -0.028108766, -0.01198095, 0.111761119, -0.076198809},
277	{0.060831772, 0.144097327, -0.069151377, 0.023754576, -0.003322955,
278	-0.071618574, 0.03353154, -0.02795295, 0.039519769, -0.023453968,
279	-0.000630308, -0.098024591, 0.017672997, 0.003813378, -0.009266499,
280	-0.011192111, 0.016013873, -0.002072968, -0.010022044, -0.012526904},
281	{-0.050776604, 0.092833081, 0.044069596, 0.050523021, -0.002628417,
282	0.076542572, -0.06388631, -0.00854892, -0.084725311, 0.017401063,
283	-0.006262541, -0.094457679, -0.002818678, -0.0044122, -0.002883973,
284	0.028729685, -0.004961596, -0.001498627, 0.017994575, -0.000232779},
285	{-0.01894566, -0.007760205, -0.015160993, -0.027254587, 0.009800903,
286	-0.013443561, -0.032896517, -0.022734138, -0.001983861, 0.00256111,
287	0.024823166, -0.021256768, 0.001980052, 0.028136263, -0.012364384,
288	-0.013782446, -0.013061091, 0.111173981, 0.021702122, 0.00046654},
289	{-0.009444193, -0.042106824, -0.02535015, -0.055125574, 0.006369612,
290	-0.02945416, -0.069922064, -0.067221068, -0.003004999, 0.053624311,
291	0.128862984, -0.057245803, 0.025550508, 0.087741073, -0.001119043,
292	-0.012036202, -0.000913488, -0.034864475, 0.050124813, 0.055534723},
293	{0.145782464, -0.024348311, -0.031216873, 0.106174443, 0.00202862,
294	0.02653866, -0.113657267, -0.00755018, 0.000307232, -0.051241158,
295	0.001310685, 0.035275877, 0.013308898, 0.002957626, -0.002925034,
296	-0.065362319, -0.071844582, 0.000475894, -0.000112419, 0.034097762},
297	{0.079840455, 0.018769331, 0.078685899, -0.084329807, -0.00277264,
298	-0.010099754, 0.059700608, -0.019209715, -0.010442992, -0.042100476,
299	-0.006020556, -0.023061786, 0.017246106, -0.001572858, -0.006703785,
300	0.056301316, -0.156787357, -0.000303638, 0.001498195, 0.051363455},
301	{0.049628261, 0.016475144, 0.094141653, -0.04444633, 0.005206131,
302	-0.001827555, 0.02195624, 0.013066683, -0.010415582, -0.022338403,
303	0.007837197, -0.023397671, -0.002507095, 0.005177694, 0.017109561,
304	-0.202340113, 0.069681441, 0.000120736, 0.002201146, 0.004670849},
305	{0.089153689, 0.000233354, 0.010826822, -0.004273519, 0.001440618,
306	0.000436077, 0.001182351, -0.002255508, -0.000700465, 0.150589876,
307	-0.003911914, -0.00050154, -0.004564983, 0.00012701, -0.001486973,
308	-0.018902754, -0.054748555, 0.000217377, -0.000319302, -0.162541651}};
309
310
311	void init_protmats()
312	{
313	long l, m;
314
315	eigmat = (double ) Malloc (20 sizeof(double));
316	for (l = 0; l <= 19; l++)
317	if (usejtt)
318	eigmat[l] = jtteigmat[l];
319	else
320	eigmat[l] = pameigmat[l];
321	probmat = (double *) Malloc (20 sizeof(double *));
322	for (l = 0; l < 20; l++)
323	probmat[l] = (double ) Malloc (20 sizeof(double));
324	for (l = 0; l <= 19; l++)
325	for (m= 0; m <= 19; m++)
326	if (usejtt)
327	probmat[l][m] = jttprobmat[l][m];
328	else
329	probmat[l][m] = pamprobmat[l][m];
330	} /* init_protmats */
331
332
333	void getoptions()
334	{
335	/* interactively set options */
336	long i, loopcount, loopcount2;
337	Char ch;
338	boolean done;
339	boolean didchangecat, didchangercat;
340	double probsum;
341
342	fprintf(outfile, "\nAmino acid sequence\n");
343	fprintf(outfile, " Maximum Likelihood method with molecular ");
344	fprintf(outfile, "clock, version %s\n\n", VERSION);
345
346	putchar('\n');
347	auto_ = false;
348	ctgry = false;
349	didchangecat = false;
350	rctgry = false;
351	didchangercat = false;
352	categs = 1;
353	rcategs = 1;
354	gama = false;
355	global = false;
356	hypstate = false;
357	invar = false;
358	jumble = false;
359	njumble = 1;
360	lambda = 1.0;
361	lambda1 = 0.0;
362	lngths = false;
363	trout = true;
364	usepam = false;
365	usejtt = true;
366	usertree = false;
367	weights = false;
368	printdata = false;
369	progress = true;
370	treeprint = true;
371	interleaved = true;
372	init_protmats();
373	loopcount = 0;
374	do {
375	cleerhome();
376	printf("\nAmino acid sequence\n");
377	printf(" Maximum Likelihood method with molecular clock, version %s\n\n",
378	VERSION);
379	printf("Settings for this run:\n");
380	printf(" U Search for best tree?");
381	if (usertree)
382	printf(" No, use user trees in input file\n");
383	else
384	printf(" Yes\n");
385	printf(" P JTT or PAM amino acid change model? %s\n",
386	(usejtt ? "Jones-Taylor-Thornton model" : "Dayhoff PAM model"));
387	if (usertree) {
388	printf(" L Use lengths from user tree?");
389	if (lngths)
390	printf(" Yes\n");
391	else
392	printf(" No\n");
393	}
394	printf(" C One category of substitution rates?");
395	if (!ctgry)
396	printf(" Yes\n");
397	else
398	printf(" %ld categories\n", categs);
399	printf(" R Rate variation among sites?");
400	if (!rctgry)
401	printf(" constant rate of change\n");
402	else {
403	if (gama)
404	printf(" Gamma distributed rates\n");
405	else {
406	if (invar)
407	printf(" Gamma+Invariant sites\n");
408	else
409	printf(" user-defined HMM of rates\n");
410	}
411	printf(" A Rates at adjacent sites correlated?");
412	if (!auto_)
413	printf(" No, they are independent\n");
414	else
415	printf(" Yes, mean block length =%6.1f\n", 1.0 / lambda);
416	}
417	if (!usertree) {
418	printf(" G Global rearrangements?");
419	if (global)
420	printf(" Yes\n");
421	else
422	printf(" No\n");
423	}
424	printf(" W Sites weighted? %s\n",
425	(weights ? "Yes" : "No"));
426	if (!usertree) {
427	printf(" J Randomize input order of sequences?");
428	if (jumble)
429	printf(" Yes (seed = %8ld, %3ld times)\n", inseed0, njumble);
430	else
431	printf(" No. Use input order\n");
432	}
433	printf(" M Analyze multiple data sets?");
434	if (mulsets)
435	printf(" Yes, %2ld %s\n", datasets,
436	(justwts ? "sets of weights" : "data sets"));
437	else
438	printf(" No\n");
439	printf(" I Input sequences interleaved?");
440	if (interleaved)
441	printf(" Yes\n");
442	else
443	printf(" No, sequential\n");
444	printf(" 0 Terminal type (IBM PC, ANSI, none)?");
445	if (ibmpc)
446	printf(" IBM PC\n");
447	if (ansi)
448	printf(" ANSI\n");
449	if (!(ibmpc \|\| ansi))
450	printf(" (none)\n");
451	printf(" 1 Print out the data at start of run");
452	if (printdata)
453	printf(" Yes\n");
454	else
455	printf(" No\n");
456	printf(" 2 Print indications of progress of run");
457	if (progress)
458	printf(" Yes\n");
459	else
460	printf(" No\n");
461	printf(" 3 Print out tree");
462	if (treeprint)
463	printf(" Yes\n");
464	else
465	printf(" No\n");
466	printf(" 4 Write out trees onto tree file?");
467	if (trout)
468	printf(" Yes\n");
469	else
470	printf(" No\n");
471	printf(" 5 Reconstruct hypothetical sequences? %s\n",
472	(hypstate ? "Yes" : "No"));
473	printf("\nAre these settings correct? ");
474	printf("(type Y or the letter for one to change)\n");
475	#ifdef WIN32
476	phyFillScreenColor();
477	#endif
478	scanf("%c%*[^\n]", &ch);
479	getchar();
480	if (ch == '\n')
481	ch = ' ';
482	uppercase(&ch);
483	done = (ch == 'Y');
484	if (!done) {
485	uppercase(&ch);
486	if (strchr("UPCRJAFWGLTMI012345", ch) != NULL){
487	switch (ch) {
488
489	case 'C':
490	ctgry = !ctgry;
491	if (ctgry) {
492	printf("\nSitewise user-assigned categories:\n\n");
493	initcatn(&categs);
494	if (rate){
495	free(rate);
496	}
497	rate = (double ) Malloc( categs sizeof(double));
498	didchangecat = true;
499	initcategs(categs, rate);
500	}
501	break;
502
503	case 'P':
504	if (usejtt) {
505	usejtt = false;
506	usepam = true;
507	} else {
508	usepam = false;
509	usejtt = true;
510	}
511	break;
512
513	case 'R':
514	if (!rctgry) {
515	rctgry = true;
516	gama = true;
517	} else {
518	if (gama) {
519	gama = false;
520	invar = true;
521	} else {
522	if (invar)
523	invar = false;
524	else
525	rctgry = false;
526	}
527	}
528	break;
529
530	case 'A':
531	auto_ = !auto_;
532	if (auto_) {
533	initlambda(&lambda);
534	lambda1 = 1.0 - lambda;
535	}
536	break;
537
538	case 'G':
539	global = !global;
540	break;
541
542	case 'W':
543	weights = !weights;
544	break;
545
546	case 'J':
547	jumble = !jumble;
548	if (jumble)
549	initjumble(&inseed, &inseed0, seed, &njumble);
550	else njumble = 1;
551	break;
552
553	case 'L':
554	lngths = !lngths;
555	break;
556
557	case 'U':
558	usertree = !usertree;
559	break;
560
561	case 'M':
562	mulsets = !mulsets;
563	if (mulsets) {
564	printf("Multiple data sets or multiple weights?");
565	loopcount2 = 0;
566	do {
567	printf(" (type D or W)\n");
568	#ifdef WIN32
569	phyFillScreenColor();
570	#endif
571	scanf("%c%*[^\n]", &ch2);
572	getchar();
573	if (ch2 == '\n')
574	ch2 = ' ';
575	uppercase(&ch2);
576	countup(&loopcount2, 10);
577	} while ((ch2 != 'W') && (ch2 != 'D'));
578	justwts = (ch2 == 'W');
579	if (justwts)
580	justweights(&datasets);
581	else
582	initdatasets(&datasets);
583	if (!jumble) {
584	jumble = true;
585	initjumble(&inseed, &inseed0, seed, &njumble);
586	}
587	}
588	break;
589
590	case 'I':
591	interleaved = !interleaved;
592	break;
593
594	case '0':
595	initterminal(&ibmpc, &ansi);
596	break;
597
598	case '1':
599	printdata = !printdata;
600	break;
601
602	case '2':
603	progress = !progress;
604	break;
605
606	case '3':
607	treeprint = !treeprint;
608	break;
609
610	case '4':
611	trout = !trout;
612	break;
613
614	case '5':
615	hypstate = !hypstate;
616	break;
617	}
618	} else
619	printf("Not a possible option!\n");
620	}
621	countup(&loopcount, 100);
622	} while (!done);
623	if (gama \|\| invar) {
624	loopcount = 0;
625	do {
626	printf(
627	"\nCoefficient of variation of substitution rate among sites (must be positive)\n");
628	printf(
629	" In gamma distribution parameters, this is 1/(square root of alpha)\n");
630	#ifdef WIN32
631	phyFillScreenColor();
632	#endif
633	scanf("%lf%*[^\n]", &cv);
634	getchar();
635	countup(&loopcount, 10);
636	} while (cv <= 0.0);
637	alpha = 1.0 / (cv * cv);
638	}
639	if (!rctgry)
640	auto_ = false;
641	if (rctgry) {
642	printf("\nRates in HMM");
643	if (invar)
644	printf(" (including one for invariant sites)");
645	printf(":\n");
646	initcatn(&rcategs);
647	if (probcat){
648	free(probcat);
649	free(rrate);
650	}
651	probcat = (double ) Malloc(rcategs sizeof(double));
652	rrate = (double ) Malloc(rcategs sizeof(double));
653	didchangercat = true;
654	if (gama)
655	initgammacat(rcategs, alpha, rrate, probcat);
656	else {
657	if (invar) {
658	loopcount = 0;
659	do {
660	printf("Fraction of invariant sites?\n");
661	scanf("%lf%*[^\n]", &invarfrac);
662	getchar();
663	countup (&loopcount, 10);
664	} while ((invarfrac <= 0.0) \|\| (invarfrac >= 1.0));
665	initgammacat(rcategs-1, alpha, rrate, probcat);
666	for (i = 0; i < rcategs-1; i++)
667	probcat[i] = probcat[i]*(1.0-invarfrac);
668	probcat[rcategs-1] = invarfrac;
669	rrate[rcategs-1] = 0.0;
670	} else {
671	initcategs(rcategs, rrate);
672	initprobcat(rcategs, &probsum, probcat);
673	}
674	}
675	}
676	if (!didchangercat){
677	rrate = Malloc( rcategs*sizeof(double));
678	probcat = Malloc( rcategs*sizeof(double));
679	rrate[0] = 1.0;
680	probcat[0] = 1.0;
681	}
682	if (!didchangecat){
683	rate = Malloc( categs*sizeof(double));
684	rate[0] = 1.0;
685	}
686	} /* getoptions */
687
688
689	void makeprotfreqs()
690	{
691	/* calculate amino acid frequencies based on eigmat */
692	long i, mineig;
693
694	mineig = 0;
695	for (i = 0; i <= 19; i++)
696	if (fabs(eigmat[i]) < fabs(eigmat[mineig]))
697	mineig = i;
698	memcpy(freqaa, probmat[mineig], 20 * sizeof(double));
699	for (i = 0; i <= 19; i++)
700	freqaa[i] = fabs(freqaa[i]);
701	} /* makeprotfreqs */
702
703
704	void allocrest()
705	{
706	long i;
707
708	y = (Char *)Malloc(sppsizeof(Char *));
709	nayme = (naym )Malloc(sppsizeof(naym));
710	for (i = 0; i < spp; i++)
711	y[i] = (char )Malloc(sites sizeof(char));
712	enterorder = (long )Malloc(sppsizeof(long));
713	weight = (long )Malloc(sitessizeof(long));
714	category = (long )Malloc(sitessizeof(long));
715	alias = (long )Malloc(sitessizeof(long));
716	aliasweight = (long )Malloc(sitessizeof(long));
717	ally = (long )Malloc(sitessizeof(long));
718	location = (long )Malloc(sitessizeof(long));
719	} /* allocrest */
720
721
722	void doinit()
723	{
724	/* initializes variables */
725
726	inputnumbers(&spp, &sites, &nonodes, 1);
727	getoptions();
728	makeprotfreqs();
729	if (printdata)
730	fprintf(outfile, "%2ld species, %3ld sites\n", spp, sites);
731	alloctree(&curtree.nodep, nonodes, usertree);
732	allocrest();
733	if (usertree)
734	return;
735	alloctree(&bestree.nodep, nonodes, 0);
736	if (njumble <= 1)
737	return;
738	alloctree(&bestree2.nodep, nonodes, 0);
739	} /* doinit */
740
741
742	void inputoptions()
743	{
744	long i;
745
746	if (!firstset)
747	samenumsp(&sites, ith);
748	if (firstset) {
749	for (i = 0; i < sites; i++)
750	category[i] = 1;
751	for (i = 0; i < sites; i++)
752	weight[i] = 1;
753	}
754	if (justwts \|\| weights)
755	inputweights(sites, weight, &weights);
756	weightsum = 0;
757	for (i = 0; i < sites; i++)
758	weightsum += weight[i];
759	if ((ctgry && categs > 1) && (firstset \|\| !justwts)) {
760	inputcategs(0, sites, category, categs, "ProMLK");
761	if (printdata)
762	printcategs(outfile, sites, category, "Site categories");
763	}
764	if (weights && printdata)
765	printweights(outfile, 0, sites, weight, "Sites");
766	fprintf(outfile, "%s model of amino acid change\n\n",
767	(usejtt ? "Jones-Taylor-Thornton" : "Dayhoff PAM"));
768	} /* inputoptions */
769
770
771	void input_protdata(long chars)
772	{
773	/* input the names and sequences for each species */
774	/* used by proml */
775	long i, j, k, l, basesread, basesnew;
776	Char charstate;
777	boolean allread, done;
778
779	if (printdata)
780	headings(chars, "Sequences", "---------");
781	basesread = 0;
782	basesnew = 0;
783	allread = false;
784	while (!(allread)) {
785	allread = true;
786	if (eoln(infile))
787	scan_eoln(infile);
788	i = 1;
789	while (i <= spp) {
790	if ((interleaved && basesread == 0) \|\| !interleaved)
791	initname(i - 1);
792	j = (interleaved) ? basesread : 0;
793	done = false;
794	while (!done && !eoff(infile)) {
795	if (interleaved)
796	done = true;
797	while (j < chars && !(eoln(infile) \|\| eoff(infile))) {
798	charstate = gettc(infile);
799	if (charstate == '\n')
800	charstate = ' ';
801	if (charstate == ' ' \|\| (charstate >= '0' && charstate <= '9'))
802	continue;
803	uppercase(&charstate);
804	if ((strchr("ABCDEFGHIKLMNPQRSTVWXYZ*?-", charstate)) == NULL){
805	printf("ERROR: bad amino acid: %c at position %ld of species %ld\n",
806	charstate, j, i);
807	if (charstate == '.') {
808	printf(" Periods (.) may not be used as gap characters.\n");
809	printf(" The correct gap character is (-)\n");
810	}
811	exxit(-1);
812	}
813	j++;
814	y[i - 1][j - 1] = charstate;
815	}
816	if (interleaved)
817	continue;
818	if (j < chars)
819	scan_eoln(infile);
820	else if (j == chars)
821	done = true;
822	}
823	if (interleaved && i == 1)
824	basesnew = j;
825
826	scan_eoln(infile);
827
828	if ((interleaved && j != basesnew) \|\|
829	(!interleaved && j != chars)) {
830	printf("ERROR: SEQUENCES OUT OF ALIGNMENT AT POSITION %ld.\n", j);
831	exxit(-1);
832	}
833	i++;
834	}
835
836	if (interleaved) {
837	basesread = basesnew;
838	allread = (basesread == chars);
839	} else
840	allread = (i > spp);
841	}
842	if (!printdata)
843	return;
844	for (i = 1; i <= ((chars - 1) / 60 + 1); i++) {
845	for (j = 1; j <= spp; j++) {
846	for (k = 0; k < nmlngth; k++)
847	putc(nayme[j - 1][k], outfile);
848	fprintf(outfile, " ");
849	l = i * 60;
850	if (l > chars)
851	l = chars;
852	for (k = (i - 1) * 60 + 1; k <= l; k++) {
853	if (j > 1 && y[j - 1][k - 1] == y[0][k - 1])
854	charstate = '.';
855	else
856	charstate = y[j - 1][k - 1];
857	putc(charstate, outfile);
858	if (k % 10 == 0 && k % 60 != 0)
859	putc(' ', outfile);
860	}
861	putc('\n', outfile);
862	}
863	putc('\n', outfile);
864	}
865	putc('\n', outfile);
866	} /* input_protdata */
867
868
869	void makeweights()
870	{
871	/* make up weights vector to avoid duplicate computations */
872	long i;
873
874	for (i = 1; i <= sites; i++) {
875	alias[i - 1] = i;
876	ally[i - 1] = 0;
877	aliasweight[i - 1] = weight[i - 1];
878	location[i - 1] = 0;
879	}
880	sitesort2(sites, aliasweight);
881	sitecombine2(sites, aliasweight);
882	sitescrunch2(sites, 1, 2, aliasweight);
883	for (i = 1; i <= sites; i++) {
884	if (aliasweight[i - 1] > 0)
885	endsite = i;
886	}
887	for (i = 1; i <= endsite; i++) {
888	ally[alias[i - 1] - 1] = alias[i - 1];
889	location[alias[i - 1] - 1] = i;
890	}
891	mp = (vall ) Malloc(sitessizeof(vall));
892	contribution = (contribarr ) Malloc( endsitesizeof(contribarr));
893	} /* makeweights */
894
895
896	void prot_makevalues(long categs, pointarray treenode, long endsite,
897	long spp, sequence y, steptr alias)
898	{
899	/* set up fractional likelihoods at tips */
900	/* a version of makevalues2 found in seq.c */
901	/* used by proml */
902	long i, j, k, l;
903	long b;
904
905	for (k = 0; k < endsite; k++) {
906	j = alias[k];
907	for (i = 0; i < spp; i++) {
908	for (l = 0; l < categs; l++) {
909	memset(treenode[i]->protx[k][l], 0, sizeof(double)*20);
910	switch (y[i][j - 1]) {
911
912	case 'A':
913	treenode[i]->protx[k][l][0] = 1.0;
914	break;
915
916	case 'R':
917	treenode[i]->protx[k][l][(long)arginine - (long)alanine] = 1.0;
918	break;
919
920	case 'N':
921	treenode[i]->protx[k][l][(long)asparagine - (long)alanine] = 1.0;
922	break;
923
924	case 'D':
925	treenode[i]->protx[k][l][(long)aspartic - (long)alanine] = 1.0;
926	break;
927
928	case 'C':
929	treenode[i]->protx[k][l][(long)cysteine - (long)alanine] = 1.0;
930	break;
931
932	case 'Q':
933	treenode[i]->protx[k][l][(long)glutamine - (long)alanine] = 1.0;
934	break;
935
936	case 'E':
937	treenode[i]->protx[k][l][(long)glutamic - (long)alanine] = 1.0;
938	break;
939
940	case 'G':
941	treenode[i]->protx[k][l][(long)glycine - (long)alanine] = 1.0;
942	break;
943
944	case 'H':
945	treenode[i]->protx[k][l][(long)histidine - (long)alanine] = 1.0;
946	break;
947
948	case 'I':
949	treenode[i]->protx[k][l][(long)isoleucine - (long)alanine] = 1.0;
950	break;
951
952	case 'L':
953	treenode[i]->protx[k][l][(long)leucine - (long)alanine] = 1.0;
954	break;
955
956	case 'K':
957	treenode[i]->protx[k][l][(long)lysine - (long)alanine] = 1.0;
958	break;
959
960	case 'M':
961	treenode[i]->protx[k][l][(long)methionine - (long)alanine] = 1.0;
962	break;
963
964	case 'F':
965	treenode[i]->protx[k][l][(long)phenylalanine - (long)alanine] = 1.0;
966	break;
967
968	case 'P':
969	treenode[i]->protx[k][l][(long)proline - (long)alanine] = 1.0;
970	break;
971
972	case 'S':
973	treenode[i]->protx[k][l][(long)serine - (long)alanine] = 1.0;
974	break;
975
976	case 'T':
977	treenode[i]->protx[k][l][(long)threonine - (long)alanine] = 1.0;
978	break;
979
980	case 'W':
981	treenode[i]->protx[k][l][(long)tryptophan - (long)alanine] = 1.0;
982	break;
983
984	case 'Y':
985	treenode[i]->protx[k][l][(long)tyrosine - (long)alanine] = 1.0;
986	break;
987
988	case 'V':
989	treenode[i]->protx[k][l][(long)valine - (long)alanine] = 1.0;
990	break;
991
992	case 'B':
993	treenode[i]->protx[k][l][(long)asparagine - (long)alanine] = 1.0;
994	treenode[i]->protx[k][l][(long)aspartic - (long)alanine] = 1.0;
995	break;
996
997	case 'Z':
998	treenode[i]->protx[k][l][(long)glutamine - (long)alanine] = 1.0;
999	treenode[i]->protx[k][l][(long)glutamic - (long)alanine] = 1.0;
1000	break;
1001
1002	case 'X': /* unknown aa */
1003	for (b = 0; b <= 19; b++)
1004	treenode[i]->protx[k][l][b] = 1.0;
1005	break;
1006
1007	case '?': /* unknown aa */
1008	for (b = 0; b <= 19; b++)
1009	treenode[i]->protx[k][l][b] = 1.0;
1010	break;
1011
1012	case '': / stop codon symbol */
1013	for (b = 0; b <= 19; b++)
1014	treenode[i]->protx[k][l][b] = 1.0;
1015	break;
1016
1017	case '-': /* deletion event-absent data or aa */
1018	for (b = 0; b <= 19; b++)
1019	treenode[i]->protx[k][l][b] = 1.0;
1020	break;
1021	}
1022	}
1023	}
1024	}
1025	} /* prot_makevalues */
1026
1027
1028	void getinput()
1029	{
1030	long grcategs;
1031
1032	/* reads the input data */
1033	if (!justwts \|\| firstset)
1034	inputoptions();
1035	if (!justwts \|\| firstset)
1036	input_protdata(sites);
1037	makeweights();
1038	setuptree2(curtree);
1039	if (!usertree) {
1040	setuptree2(bestree);
1041	if (njumble > 1)
1042	setuptree2(bestree2);
1043	}
1044	grcategs = (categs > rcategs) ? categs : rcategs;
1045	prot_allocx(nonodes, grcategs, curtree.nodep, usertree);
1046	if (!usertree) {
1047	prot_allocx(nonodes, grcategs, bestree.nodep, 0);
1048	if (njumble > 1)
1049	prot_allocx(nonodes, grcategs, bestree2.nodep, 0);
1050	}
1051	prot_makevalues(rcategs, curtree.nodep, endsite, spp, y, alias);
1052	} /* getinput */
1053
1054
1055	void prot_inittable()
1056	{
1057	/* Define a lookup table. Precompute values and print them out in tables */
1058	/* Allocate memory for the pmatrices, dpmatices and ddpmatrices */
1059	long i, j, k, l;
1060	double sumrates;
1061
1062	/* Allocate memory for pmatrices, the array of pointers to pmatrices */
1063
1064	pmatrices = (double ****) Malloc (spp sizeof(double ****));
1065
1066	/* Allocate memory for the first 2 pmatrices, the matrix of conversion */
1067	/* probabilities, but only once per run (aka not on the second jumble. */
1068
1069	alloc_pmatrix(0);
1070	alloc_pmatrix(1);
1071
1072	/* Allocate memory for one dpmatrix, the first derivative matrix */
1073
1074	dpmatrix = (double ***) Malloc( rcategs sizeof(double ***));
1075	for (j = 0; j < rcategs; j++) {
1076	dpmatrix[j] = (double **) Malloc( categs sizeof(double **));
1077	for (k = 0; k < categs; k++) {
1078	dpmatrix[j][k] = (double *) Malloc( 20 sizeof(double *));
1079	for (l = 0; l < 20; l++)
1080	dpmatrix[j][k][l] = (double ) Malloc( 20 sizeof(double));
1081	}
1082	}
1083
1084	/* Allocate memory for one ddpmatrix, the second derivative matrix */
1085	ddpmatrix = (double ***) Malloc( rcategs sizeof(double ***));
1086	for (j = 0; j < rcategs; j++) {
1087	ddpmatrix[j] = (double **) Malloc( categs sizeof(double **));
1088	for (k = 0; k < categs; k++) {
1089	ddpmatrix[j][k] = (double *) Malloc( 20 sizeof(double *));
1090	for (l = 0; l < 20; l++)
1091	ddpmatrix[j][k][l] = (double ) Malloc( 20 sizeof(double));
1092	}
1093	}
1094
1095	/* Allocate memory and assign values to tbl, the matrix of possible rates*/
1096
1097	tbl = (double *) Malloc( rcategs sizeof(double *));
1098	for (j = 0; j < rcategs; j++)
1099	tbl[j] = (double ) Malloc( categs sizeof(double));
1100
1101	for (j = 0; j < rcategs; j++)
1102	for (k = 0; k < categs; k++)
1103	tbl[j][k] = rrate[j]*rate[k];
1104
1105	sumrates = 0.0;
1106	for (i = 0; i < endsite; i++) {
1107	for (j = 0; j < rcategs; j++)
1108	sumrates += aliasweight[i] * probcat[j]
1109	* tbl[j][category[alias[i] - 1] - 1];
1110	}
1111	sumrates /= (double)sites;
1112	for (j = 0; j < rcategs; j++)
1113	for (k = 0; k < categs; k++) {
1114	tbl[j][k] /= sumrates;
1115	}
1116	if (gama \|\| invar) {
1117	fprintf(outfile, "\nDiscrete approximation to gamma distributed rates\n");
1118	fprintf(outfile,
1119	" Coefficient of variation of rates = %f (alpha = %f)\n", cv, alpha);
1120	}
1121	if (rcategs > 1) {
1122	fprintf(outfile, "\nState in HMM Rate of change Probability\n\n");
1123	for (i = 0; i < rcategs; i++)
1124	if (probcat[i] < 0.0001)
1125	fprintf(outfile, "%9ld%16.3f%20.6f\n", i+1, rrate[i], probcat[i]);
1126	else if (probcat[i] < 0.001)
1127	fprintf(outfile, "%9ld%16.3f%19.5f\n", i+1, rrate[i], probcat[i]);
1128	else if (probcat[i] < 0.01)
1129	fprintf(outfile, "%9ld%16.3f%18.4f\n", i+1, rrate[i], probcat[i]);
1130	else
1131	fprintf(outfile, "%9ld%16.3f%17.3f\n", i+1, rrate[i], probcat[i]);
1132	putc('\n', outfile);
1133	if (auto_) {
1134	fprintf(outfile,
1135	"Expected length of a patch of sites having the same rate = %8.3f\n",
1136	1/lambda);
1137	putc('\n', outfile);
1138	}
1139	}
1140	if (categs > 1) {
1141	fprintf(outfile, "\nSite category Rate of change\n\n");
1142	for (k = 0; k < categs; k++)
1143	fprintf(outfile, "%9ld%16.3f\n", k+1, rate[k]);
1144	fprintf(outfile, "\n\n");
1145	}
1146	} /* prot_inittable */
1147
1148
1149	void alloc_pmatrix(long sib)
1150	{
1151	/* Allocate memory for a new pmatrix. Called iff num_sibs>max_num_sibs */
1152	long j, k, l;
1153	double ****temp_matrix;
1154
1155	temp_matrix = (double ***) Malloc (rcategs sizeof(double ***));
1156	for (j = 0; j < rcategs; j++) {
1157	temp_matrix[j] = (double **) Malloc(categs sizeof(double **));
1158	for (k = 0; k < categs; k++) {
1159	temp_matrix[j][k] = (double *) Malloc(20 sizeof (double *));
1160	for (l = 0; l < 20; l++)
1161	temp_matrix[j][k][l] = (double ) Malloc(20 sizeof(double));
1162	}
1163	}
1164	pmatrices[sib] = temp_matrix;
1165	max_num_sibs++;
1166	} /* alloc_pmatrix */
1167
1168
1169	void make_pmatrix(double matrix, double dmat, double **ddmat,
1170	long derivative, double lz, double rat,
1171	double eigmat, double *probmat)
1172	{
1173	/* Computes the R matrix such that matrix[m][l] is the joint probability */
1174	/* of m and l. */
1175	/* Computes a P matrix such that matrix[m][l] is the conditional */
1176	/* probability of m given l. This is accomplished by dividing all terms */
1177	/* in the R matrix by freqaa[m], the frequency of l. */
1178
1179	long k, l, m; /* (l) original character state */
1180	/* (m) final character state */
1181	/* (k) lambda counter */
1182	double p0, p1, p2, q;
1183	double elambdat[20], delambdat[20], ddelambdat[20];
1184	/* exponential term for matrix */
1185	/* and both derivative matrices */
1186
1187	for (k = 0; k <= 19; k++) {
1188	elambdat[k] = exp(lz * rat * eigmat[k]);
1189	if(derivative != 0) {
1190	delambdat[k] = (elambdat[k] * rat * eigmat[k]);
1191	ddelambdat[k] = (delambdat[k] * rat * eigmat[k]);
1192	}
1193	}
1194	for (m = 0; m <= 19; m++) {
1195	for (l = 0; l <= 19; l++) {
1196	p0 = 0.0;
1197	p1 = 0.0;
1198	p2 = 0.0;
1199	for (k = 0; k <= 19; k++) {
1200	q = probmat[k][m] * probmat[k][l];
1201	p0 += (q * elambdat[k]);
1202	if(derivative !=0) {
1203	p1 += (q * delambdat[k]);
1204	p2 += (q * ddelambdat[k]);
1205	}
1206	}
1207	matrix[m][l] = p0 / freqaa[m];
1208	if(derivative != 0) {
1209	dmat[m][l] = p1 / freqaa[m];
1210	ddmat[m][l] = p2 / freqaa[m];
1211	}
1212	}
1213	}
1214	} /* make_pmatrix */
1215
1216
1217	void prot_nuview(node *p)
1218	{
1219	long b, i, j, k, l, m, num_sibs, sib_index;
1220	node sib_ptr, sib_back_ptr;
1221	psitelike prot_xx, x2;
1222	double lw, prod7;
1223	double **pmat;
1224
1225	/* Figure out how many siblings the current node has */
1226	/* and be sure that pmatrices is large enough */
1227	num_sibs = count_sibs(p);
1228	for (i = 0; i < num_sibs; i++)
1229	if (pmatrices[i] == NULL)
1230	alloc_pmatrix(i);
1231
1232	/* Recursive calls, should be called for all children */
1233	sib_ptr = p;
1234	for (i=0 ; i < num_sibs; i++) {
1235	sib_ptr = sib_ptr->next;
1236	sib_back_ptr = sib_ptr->back;
1237
1238	if (!(sib_back_ptr == NULL))
1239	if (!sib_back_ptr->tip && !sib_back_ptr->initialized)
1240	prot_nuview(sib_back_ptr);
1241	}
1242
1243	/* Make pmatrices for all possible combinations of category, rcateg */
1244	/* and sib */
1245	sib_ptr = p; /* return to p */
1246	for (sib_index=0; sib_index < num_sibs; sib_index++) {
1247	sib_ptr = sib_ptr->next;
1248	sib_back_ptr = sib_ptr->back;
1249
1250	if (sib_back_ptr != NULL)
1251	lw = fabs(p->tyme - sib_back_ptr->tyme);
1252	else
1253	lw = 0.0;
1254
1255	for (j = 0; j < rcategs; j++)
1256	for (k = 0; k < categs; k++)
1257	make_pmatrix(pmatrices[sib_index][j][k], NULL, NULL, 0, lw,
1258	tbl[j][k], eigmat, probmat);
1259	}
1260
1261	for (i = 0; i < endsite; i++) {
1262	k = category[alias[i]-1] - 1;
1263	for (j = 0; j < rcategs; j++) {
1264
1265	/* initialize to 1 all values of prot_xx */
1266	for (m = 0; m <= 19; m++)
1267	prot_xx[m] = 1;
1268
1269	sib_ptr = p; /* return to p */
1270	/* loop through all sibs and calculate likelihoods for all possible*/
1271	/* amino acid combinations */
1272	for (sib_index=0; sib_index < num_sibs; sib_index++) {
1273	sib_ptr = sib_ptr->next;
1274	sib_back_ptr = sib_ptr->back;
1275
1276	if (sib_back_ptr != NULL)
1277	memcpy(x2, sib_back_ptr->protx[i][j], sizeof(psitelike));
1278	else
1279	for (b = 0; b <= 19; b++)
1280	x2[b] = 1.0;
1281	pmat = pmatrices[sib_index][j][k];
1282	for (m = 0; m <= 19; m++) {
1283	prod7 = 0;
1284	for (l = 0; l <= 19; l++)
1285	prod7 += (pmat[m][l] * x2[l]);
1286	prot_xx[m] *= prod7;
1287	}
1288	}
1289	/* And the final point of this whole function: */
1290	memcpy(p->protx[i][j], prot_xx, sizeof(psitelike));
1291	}
1292	}
1293
1294	p->initialized = true;
1295	} /* prot_nuview */
1296
1297
1298	void getthree(node *p, double thigh, double tlow)
1299	{
1300	/* compute likelihood at a new triple of points */
1301	int i;
1302	double tt = p->tyme;
1303	double td = fabs(tdelta);
1304
1305	x[0] = tt - td;
1306	x[1] = tt;
1307	x[2] = tt + td;
1308
1309	if ( x[0] < tlow + epsilon ) {
1310	x[0] = tlow + epsilon;
1311	x[1] = ( x[0] + x[2] ) / 2;
1312	}
1313
1314	if ( x[2] > thigh - epsilon ) {
1315	x[2] = thigh - epsilon;
1316	x[1] = ( x[0] + x[2] ) / 2;
1317	}
1318
1319	for ( i = 0 ; i < 3 ; i++ ) {
1320	p->tyme = x[i];
1321	prot_nuview(p);
1322	lnl[i] = prot_evaluate(p);
1323	}
1324	} /* getthree */
1325
1326	void makenewv(node *p)
1327	{
1328	/* improve a node time */
1329	long it, imin, imax, i, num_sibs;
1330	double tt, tfactor, tlow, thigh, oldlike, ymin, ymax, s32, s21, yold;
1331	boolean done, already;
1332	node s, sdown, sib_ptr, sib_back_ptr;
1333
1334	s = curtree.nodep[p->index - 1];
1335	sdown = s->back;
1336	if (s == curtree.root)
1337	tlow = -10.0;
1338	else
1339	tlow = sdown->tyme;
1340
1341	sib_ptr = s;
1342	num_sibs = count_sibs(p);
1343
1344	thigh = s->next->back->tyme;
1345	for (i=0 ; i < num_sibs; i++) {
1346	sib_ptr = sib_ptr->next;
1347	sib_back_ptr = sib_ptr->back;
1348	if (sib_back_ptr->tyme < thigh)
1349	thigh = sib_back_ptr->tyme;
1350	}
1351	done = (thigh - tlow < 4.0*epsilon);
1352	it = 1;
1353	if (s != curtree.root)
1354	tdelta = (thigh - tlow) / 10.0;
1355	else
1356	tdelta = (thigh - s->tyme) / 5.0;
1357	tfactor = 1.0;
1358	if (!done)
1359	getthree(s, thigh, tlow);
1360	while (it < iterations && !done) {
1361	tt = s->tyme;
1362	ymax = lnl[0];
1363	imax = 1;
1364	for (i = 2; i <= 3; i++) {
1365	if (lnl[i - 1] > ymax) {
1366	ymax = lnl[i - 1];
1367	imax = i;
1368	}
1369	}
1370	if (imax != 2) {
1371	ymax = x[1];
1372	x[1] = x[imax - 1];
1373	x[imax - 1] = ymax;
1374	ymax = lnl[1];
1375	lnl[1] = lnl[imax - 1];
1376	lnl[imax - 1] = ymax;
1377	}
1378	tt = x[1];
1379	oldlike = lnl[1];
1380	yold = tt;
1381	s32 = (lnl[2] - lnl[1]) / (x[2] - x[1]);
1382	s21 = (lnl[1] - lnl[0]) / (x[1] - x[0]);
1383	if (fabs(x[2] - x[0]) > epsilon)
1384	curv = (s32 - s21) / ((x[2] - x[0]) / 2);
1385	else
1386	curv = 0.0;
1387	slope = (s32 + s21) / 2 - curv * (x[2] - 2 * x[1] + x[0]) / 4;
1388	if (curv >= 0.0) {
1389	if (slope < 0)
1390	tdelta = -fabs(tdelta);
1391	else
1392	tdelta = fabs(tdelta);
1393	} else
1394	tdelta = -(tfactor * slope / curv);
1395	if (tt + tdelta <= tlow + epsilon)
1396	tdelta = tlow + epsilon - tt;
1397	if (tt + tdelta >= thigh - epsilon)
1398	tdelta = thigh - epsilon - tt;
1399	tt += tdelta;
1400	done = (fabs(yold - tt) < epsilon \|\| fabs(tdelta) < epsilon);
1401	s->tyme = tt;
1402	prot_nuview(s);
1403	lnlike = prot_evaluate(s);
1404	ymin = lnl[0];
1405	imin = 1;
1406	for (i = 2; i <= 3; i++) {
1407	if (lnl[i - 1] < ymin) {
1408	ymin = lnl[i - 1];
1409	imin = i;
1410	}
1411	}
1412	already = (tt == x[0]) \|\| (tt == x[1]) \|\| (tt == x[2]);
1413	if (!already && ymin < lnlike) {
1414	x[imin - 1] = tt;
1415	lnl[imin - 1] = lnlike;
1416	}
1417	if (already \|\| lnlike < oldlike) {
1418	tt = x[1];
1419	tfactor /= 2;
1420	tdelta /= 2;
1421	curtree.likelihood = oldlike;
1422	lnlike = oldlike;
1423	} else
1424	tfactor = 1.0;
1425
1426	if (!done) {
1427	sib_ptr = p;
1428	num_sibs = count_sibs(p);
1429	p->tyme = tt;
1430	for (i=0 ; i < num_sibs; i++) {
1431	sib_ptr = sib_ptr->next;
1432	sib_ptr->tyme = tt;
1433	}
1434
1435	sib_ptr = p;
1436	prot_nuview(p);
1437	for (i=0 ; i < num_sibs; i++) {
1438	sib_ptr = sib_ptr->next;
1439	prot_nuview(sib_ptr);
1440	}
1441	}
1442
1443	it++;
1444	}
1445	sib_ptr = p;
1446	for (i=0 ; i < num_sibs; i++) {
1447	sib_ptr = sib_ptr->next;
1448	inittrav (sib_ptr);
1449	}
1450	smoothed = smoothed && done;
1451	} /* makenewv */
1452
1453
1454	void update(node *p)
1455	{
1456	node sib_ptr, sib_back_ptr;
1457	long i, num_sibs;
1458
1459	/* improve time and recompute views at a node */
1460	if (p == NULL)
1461	return;
1462	if (p->back != NULL) {
1463	if (!p->back->tip && !p->back->initialized)
1464	prot_nuview(p->back);
1465	}
1466
1467	sib_ptr = p;
1468	num_sibs = count_sibs(p);
1469	for (i=0 ; i < num_sibs; i++) {
1470	sib_ptr = sib_ptr->next;
1471	sib_back_ptr = sib_ptr->back;
1472	if (sib_back_ptr != NULL) {
1473	if (!sib_back_ptr->tip && !sib_back_ptr->initialized)
1474	prot_nuview(sib_back_ptr);
1475	}
1476	}
1477
1478	if ((!usertree) \|\| (usertree && !lngths) \|\| p->iter) {
1479	makenewv(p);
1480	return;
1481	}
1482	prot_nuview(p);
1483
1484	sib_ptr = p;
1485	num_sibs = count_sibs(p);
1486	for (i=0 ; i < num_sibs; i++) {
1487	sib_ptr = sib_ptr->next;
1488	prot_nuview(sib_ptr);
1489	}
1490	} /* update */
1491
1492
1493	void smooth(node *p)
1494	{
1495	node *sib_ptr;
1496	long i, num_sibs;
1497
1498	if (p == NULL)
1499	return;
1500	if (p->tip)
1501	return;
1502
1503	update(p);
1504
1505	smoothed = false;
1506	sib_ptr = p;
1507	num_sibs = count_sibs(p);
1508	for (i=0; i < num_sibs; i++) {
1509	sib_ptr = sib_ptr->next;
1510	if (polishing \|\| (smoothit && !smoothed)) {
1511	smooth(sib_ptr->back);
1512	p->initialized = false;
1513	sib_ptr->initialized = false;
1514	}
1515	update(p);
1516	}
1517	} /* smooth */
1518
1519
1520	void promlk_add(node below, node newtip, node *newfork, boolean tempadd)
1521	{
1522	/* inserts the nodes newfork and its descendant, newtip, into the tree. */
1523	long i;
1524	boolean done;
1525	node *p;
1526
1527	below = curtree.nodep[below->index - 1];
1528	newfork = curtree.nodep[newfork->index-1];
1529	newtip = curtree.nodep[newtip->index-1];
1530	if (below->back != NULL)
1531	below->back->back = newfork;
1532	newfork->back = below->back;
1533	below->back = newfork->next->next;
1534	newfork->next->next->back = below;
1535	newfork->next->back = newtip;
1536	newtip->back = newfork->next;
1537	if (newtip->tyme < below->tyme)
1538	p = newtip;
1539	else p = below;
1540	newfork->tyme = p->tyme;
1541	if (curtree.root == below)
1542	curtree.root = newfork;
1543	if (newfork->back != NULL) {
1544	if (p->tyme > newfork->back->tyme)
1545	newfork->tyme = (p->tyme + newfork->back->tyme) / 2.0;
1546	else newfork->tyme = p->tyme - epsilon;
1547	newfork->next->tyme = newfork->tyme;
1548	newfork->next->next->tyme = newfork->tyme;
1549	do {
1550	p = curtree.nodep[p->back->index - 1];
1551	done = (p == curtree.root);
1552	if (!done)
1553	done = (curtree.nodep[p->back->index - 1]->tyme < p->tyme - epsilon);
1554	if (!done) {
1555	curtree.nodep[p->back->index - 1]->tyme = p->tyme - epsilon;
1556	curtree.nodep[p->back->index - 1]->next->tyme = p->tyme - epsilon;
1557	curtree.nodep[p->back->index - 1]->next->next->tyme = p->tyme - epsilon;
1558	}
1559	} while (!done);
1560	} else {
1561	newfork->tyme = newfork->tyme - 2*epsilon;
1562	newfork->next->tyme = newfork->tyme;
1563	newfork->next->next->tyme = newfork->tyme;
1564	}
1565	inittrav(newtip);
1566	inittrav(newtip->back);
1567	smoothed = false;
1568	i = 1;
1569	while (i < smoothings && !smoothed) {
1570	smoothed = true;
1571	smooth(newfork);
1572	smooth(newfork->back);
1573	i++;
1574	}
1575	} /* promlk_add */
1576
1577
1578	void promlk_re_move(node item, node fork, boolean tempadd)
1579	{
1580	/* removes nodes item and its ancestor, fork, from the tree.
1581	the new descendant of fork's ancestor is made to be
1582	fork's second descendant (other than item). Also
1583	returns pointers to the deleted nodes, item and fork */
1584	node p, q;
1585	long i;
1586
1587	if ((*item)->back == NULL) {
1588	*fork = NULL;
1589	return;
1590	}
1591	item = curtree.nodep[(item)->index-1];
1592	fork = curtree.nodep[(item)->back->index - 1];
1593	if (curtree.root == *fork) {
1594	if (item == (fork)->next->back)
1595	curtree.root = (*fork)->next->next->back;
1596	else
1597	curtree.root = (*fork)->next->back;
1598	}
1599	p = (*item)->back->next->back;
1600	q = (*item)->back->next->next->back;
1601	/* debug replace by hookup calls? Does that have NULL protection? */
1602	if (p != NULL)
1603	p->back = q;
1604	if (q != NULL)
1605	q->back = p;
1606	(*fork)->back = NULL;
1607	p = (*fork)->next;
1608	while (p != *fork) {
1609	p->back = NULL;
1610	p = p->next;
1611	}
1612	(*item)->back = NULL;
1613	inittrav(p);
1614	inittrav(q);
1615	if (tempadd)
1616	return;
1617	i = 1;
1618	while (i <= smoothings) {
1619	smooth(q);
1620	if (smoothit)
1621	smooth(q->back);
1622	i++;
1623	}
1624	} /* promlk_re_move */
1625
1626
1627	double prot_evaluate(node *p)
1628	{
1629	contribarr tterm;
1630	static contribarr like, nulike, clai;
1631	double sum, sum2, sumc=0, y, prod4, prodl, frexm, sumterm, lterm;
1632	double **pmat;
1633	long i, j, k, l, m, lai;
1634	node q, r;
1635	psitelike x1, x2;
1636
1637	sum = 0.0;
1638
1639	if (p == curtree.root && (count_sibs(p) == 2)) {
1640	r = p->next->back;
1641	q = p->next->next->back;
1642	y = r->tyme + q->tyme - 2 * p->tyme;
1643	if (!r->tip && !r->initialized) prot_nuview (r);
1644	if (!q->tip && !q->initialized) prot_nuview (q);
1645	} else if (p == curtree.root) {
1646	/* the next two lines copy tyme and x to p->next. Normally they are
1647	not initialized for an internal node. */
1648	/* assumes bifurcation */
1649	p->next->tyme = p->tyme;
1650	prot_nuview(p->next);
1651	r = p->next;
1652	q = p->next->back;
1653	y = fabs(p->next->tyme - q->tyme);
1654	} else {
1655	r = p;
1656	q = p->back;
1657	if (!r->tip && !r->initialized) prot_nuview (r);
1658	if (!q->tip && !q->initialized) prot_nuview (q);
1659	y = fabs(r->tyme - q->tyme);
1660	}
1661
1662	for (j = 0; j < rcategs; j++)
1663	for (k = 0; k < categs; k++)
1664	make_pmatrix(pmatrices[0][j][k],NULL,NULL,0,y,tbl[j][k],eigmat,probmat);
1665	for (i = 0; i < endsite; i++) {
1666	k = category[alias[i]-1] - 1;
1667	for (j = 0; j < rcategs; j++) {
1668	memcpy(x1, r->protx[i][j], sizeof(psitelike));
1669	memcpy(x2, q->protx[i][j], sizeof(psitelike));
1670	prod4 = 0.0;
1671	pmat = pmatrices[0][j][k];
1672	for (m = 0; m <= 19; m++) {
1673	prodl = 0.0;
1674	for (l = 0; l <= 19; l++)
1675	prodl += (pmat[m][l] * x2[l]);
1676	frexm = x1[m] * freqaa[m];
1677	prod4 += (prodl * frexm);
1678	}
1679	tterm[j] = prod4;
1680	}
1681	sumterm = 0.0;
1682	for (j = 0; j < rcategs; j++)
1683	sumterm += probcat[j] * tterm[j];
1684	if (sumterm < 0.0)
1685	sumterm = 0.00000001; /* ??? */
1686	lterm = log(sumterm);
1687	for (j = 0; j < rcategs; j++)
1688	clai[j] = tterm[j] / sumterm;
1689	memcpy(contribution[i], clai, rcategs * sizeof(double));
1690	if (!auto_ && usertree)
1691	l0gf[which - 1][i] = lterm;
1692	sum += aliasweight[i] * lterm;
1693	}
1694	if (auto_) {
1695	for (j = 0; j < rcategs; j++)
1696	like[j] = 1.0;
1697	for (i = 0; i < sites; i++) {
1698	if ((ally[i] > 0) && (location[ally[i]-1] > 0)) {
1699	sumc = 0.0;
1700	for (k = 0; k < rcategs; k++)
1701	sumc += probcat[k] * like[k];
1702	sumc *= lambda;
1703	lai = location[ally[i] - 1];
1704	memcpy(clai, contribution[lai - 1], rcategs*sizeof(double));
1705	for (j = 0; j < rcategs; j++)
1706	nulike[j] = ((1.0 - lambda) * like[j] + sumc) * clai[j];
1707	} else {
1708	for (j = 0; j < rcategs; j++)
1709	nulike[j] = ((1.0 - lambda) * like[j] + sumc);
1710	}
1711	memcpy(like, nulike, rcategs * sizeof(double));
1712	}
1713	sum2 = 0.0;
1714	for (i = 0; i < rcategs; i++)
1715	sum2 += probcat[i] * like[i];
1716	sum += log(sum2);
1717	}
1718	curtree.likelihood = sum;
1719	if (auto_ \|\| !usertree)
1720	return sum;
1721	l0gl[which - 1] = sum;
1722	if (which == 1) {
1723	maxwhich = 1;
1724	maxlogl = sum;
1725	return sum;
1726	}
1727	if (sum > maxlogl) {
1728	maxwhich = which;
1729	maxlogl = sum;
1730	}
1731	return sum;
1732	} /* prot_evaluate */
1733
1734
1735	void tryadd(node p, node item, node *nufork)
1736	{ /* temporarily adds one fork and one tip to the tree.
1737	if the location where they are added yields greater
1738	likelihood than other locations tested up to that
1739	time, then keeps that location as there */
1740
1741	long grcategs;
1742	grcategs = (categs > rcategs) ? categs : rcategs;
1743
1744	promlk_add(p, item, nufork, true);
1745	like = prot_evaluate(p);
1746	if (lastsp) {
1747	if (like >= bestyet \|\| bestyet == UNDEFINED)
1748	prot_copy_(&curtree, &bestree, nonodes, grcategs);
1749	}
1750	if (like > bestyet \|\| bestyet == UNDEFINED) {
1751	bestyet = like;
1752	there = p;
1753	}
1754	promlk_re_move(item, nufork, true);
1755	} /* tryadd */
1756
1757
1758	void addpreorder(node p, node item_, node *nufork_, boolean contin,
1759	boolean continagain)
1760	{
1761	/* traverses a binary tree, calling function tryadd
1762	at a node before calling tryadd at its descendants */
1763	node item, nufork;
1764
1765	item = item_;
1766	nufork = nufork_;
1767	if (p == NULL)
1768	return;
1769	tryadd(p, &item, &nufork);
1770	contin = continagain;
1771	if ((!p->tip) && contin) {
1772	addpreorder(p->next->back, item, nufork, contin, continagain);
1773	addpreorder(p->next->next->back, item, nufork, contin, continagain);
1774	}
1775	} /* addpreorder */
1776
1777
1778	void restoradd(node below, node newtip, node *newfork, double prevtyme)
1779	{
1780	/* restore "new" tip and fork to place "below". restore tymes */
1781	/* assumes bifurcation */
1782	hookup(newfork, below->back);
1783	hookup(newfork->next, below);
1784	hookup(newtip, newfork->next->next);
1785	curtree.nodep[newfork->index-1] = newfork;
1786	newfork->tyme = prevtyme;
1787	/* assumes bifurcations */
1788	newfork->next->tyme = prevtyme;
1789	newfork->next->next->tyme = prevtyme;
1790	} /* restoradd */
1791
1792
1793	void tryrearr(node p, boolean success)
1794	{
1795	/* evaluates one rearrangement of the tree.
1796	if the new tree has greater likelihood than the old
1797	one sets success = TRUE and keeps the new tree.
1798	otherwise, restores the old tree */
1799	node frombelow, whereto, *forknode;
1800	double oldlike, prevtyme;
1801	boolean wasonleft;
1802
1803	if (p == curtree.root)
1804	return;
1805	forknode = curtree.nodep[p->back->index - 1];
1806	if (forknode == curtree.root)
1807	return;
1808	oldlike = bestyet;
1809	prevtyme = forknode->tyme;
1810	/* the following statement presumes bifurcating tree */
1811	if (forknode->next->back == p) {
1812	frombelow = forknode->next->next->back;
1813	wasonleft = true;
1814	}
1815	else {
1816	frombelow = forknode->next->back;
1817	wasonleft = false;
1818	}
1819	whereto = curtree.nodep[forknode->back->index - 1];
1820	promlk_re_move(&p, &forknode, true);
1821	promlk_add(whereto, p, forknode, true);
1822	like = prot_evaluate(p);
1823	if (like <= oldlike && oldlike != UNDEFINED) {
1824	promlk_re_move(&p, &forknode, true);
1825	restoradd(frombelow, p, forknode, prevtyme);
1826	if (wasonleft && (forknode->next->next->back == p)) {
1827	hookup (forknode->next->back, forknode->next->next);
1828	hookup (forknode->next, p);
1829	}
1830	curtree.likelihood = oldlike;
1831	inittrav(forknode);
1832	inittrav(forknode->next);
1833	inittrav(forknode->next->next);
1834	} else {
1835	(*success) = true;
1836	bestyet = like;
1837	}
1838	} /* tryrearr */
1839
1840
1841	void repreorder(node p, boolean success)
1842	{
1843	/* traverses a binary tree, calling function tryrearr
1844	at a node before calling tryrearr at its descendants */
1845	if (p == NULL)
1846	return;
1847	tryrearr(p, success);
1848	if (p->tip)
1849	return;
1850	if (!(*success))
1851	repreorder(p->next->back, success);
1852	if (!(*success))
1853	repreorder(p->next->next->back, success);
1854	} /* repreorder */
1855
1856
1857	void rearrange(node **r)
1858	{
1859	/* traverses the tree (preorder), finding any local
1860	rearrangement which increases the likelihood.
1861	if traversal succeeds in increasing the tree's
1862	likelihood, function rearrange runs traversal again */
1863	boolean success;
1864	success = true;
1865	while (success) {
1866	success = false;
1867	repreorder(*r, &success);
1868	}
1869	} /* rearrange */
1870
1871
1872	void nodeinit(node *p)
1873	{
1874	/* set up times at one node */
1875	node sib_ptr, sib_back_ptr;
1876	long i, num_sibs;
1877	double lowertyme;
1878
1879	sib_ptr = p;
1880	num_sibs = count_sibs(p);
1881
1882	/* lowertyme = lowest of children's times */
1883	lowertyme = p->next->back->tyme;
1884	for (i=0 ; i < num_sibs; i++) {
1885	sib_ptr = sib_ptr->next;
1886	sib_back_ptr = sib_ptr->back;
1887	if (sib_back_ptr->tyme < lowertyme)
1888	lowertyme = sib_back_ptr->tyme;
1889	}
1890
1891	p->tyme = lowertyme - 0.1;
1892
1893	sib_ptr = p;
1894	for (i=0 ; i < num_sibs; i++) {
1895	sib_ptr = sib_ptr->next;
1896	sib_back_ptr = sib_ptr->back;
1897
1898	sib_ptr->tyme = p->tyme;
1899	sib_back_ptr->v = sib_back_ptr->tyme - p->tyme;
1900	sib_ptr->v = sib_back_ptr->v;
1901	}
1902	} /* nodeinit */
1903
1904
1905	void initrav(node *p)
1906	{
1907
1908	long i, num_sibs;
1909	node sib_ptr, sib_back_ptr;
1910
1911	/* traverse to set up times throughout tree */
1912	if (p->tip)
1913	return;
1914
1915	sib_ptr = p;
1916	num_sibs = count_sibs(p);
1917	for (i=0 ; i < num_sibs; i++) {
1918	sib_ptr = sib_ptr->next;
1919	sib_back_ptr = sib_ptr->back;
1920	initrav(sib_back_ptr);
1921	}
1922
1923	nodeinit(p);
1924	} /* initrav */
1925
1926
1927	void travinit(node *p)
1928	{
1929	long i, num_sibs;
1930	node sib_ptr, sib_back_ptr;
1931
1932	/* traverse to set up initial values */
1933	if (p == NULL)
1934	return;
1935	if (p->tip)
1936	return;
1937	if (p->initialized)
1938	return;
1939
1940
1941	sib_ptr = p;
1942	num_sibs = count_sibs(p);
1943	for (i=0 ; i < num_sibs; i++) {
1944	sib_ptr = sib_ptr->next;
1945	sib_back_ptr = sib_ptr->back;
1946	travinit(sib_back_ptr);
1947	}
1948
1949	prot_nuview(p);
1950	p->initialized = true;
1951	} /* travinit */
1952
1953
1954	void travsp(node *p)
1955	{
1956	long i, num_sibs;
1957	node sib_ptr, sib_back_ptr;
1958
1959	/* traverse to find tips */
1960	if (p == curtree.root)
1961	travinit(p);
1962	if (p->tip)
1963	travinit(p->back);
1964	else {
1965	sib_ptr = p;
1966	num_sibs = count_sibs(p);
1967	for (i=0 ; i < num_sibs; i++) {
1968	sib_ptr = sib_ptr->next;
1969	sib_back_ptr = sib_ptr->back;
1970	travsp(sib_back_ptr);
1971	}
1972	}
1973	} /* travsp */
1974
1975
1976	void treevaluate()
1977	{
1978	/* evaluate likelihood of tree, after iterating branch lengths */
1979	long i, j, num_sibs;
1980	node sib_ptr, sib_back_ptr;
1981	double dummy;
1982
1983	polishing = true;
1984	smoothit = true;
1985	for (i = 0; i < spp; i++)
1986	curtree.nodep[i]->initialized = false;
1987	for (i = spp; i < nonodes; i++) {
1988	sib_ptr = curtree.nodep[i];
1989	sib_ptr->initialized = false;
1990	num_sibs = count_sibs(sib_ptr);
1991	for (j=0 ; j < num_sibs; j++) {
1992	sib_ptr = sib_ptr->next;
1993	sib_back_ptr = sib_ptr->back;
1994	sib_ptr->initialized = false;
1995	}
1996
1997	}
1998	if (!lngths)
1999	initrav(curtree.root);
2000	travsp(curtree.root);
2001	for (i = 1; i <= smoothings * 4; i++)
2002	smooth(curtree.root);
2003	dummy = prot_evaluate(curtree.root);
2004	} /* treevaluate */
2005
2006
2007	void promlk_coordinates(node p, long tipy)
2008	{
2009	/* establishes coordinates of nodes */
2010	node q, first, last, pp1 =NULL, *pp2 =NULL;
2011	long num_sibs, p1, p2, i;
2012
2013	if (p->tip) {
2014	p->xcoord = 0;
2015	p->ycoord = (*tipy);
2016	p->ymin = (*tipy);
2017	p->ymax = (*tipy);
2018	(*tipy) += down;
2019	return;
2020	}
2021	q = p->next;
2022	do {
2023	promlk_coordinates(q->back, tipy);
2024	q = q->next;
2025	} while (p != q);
2026	num_sibs = count_sibs(p);
2027	p1 = (long)((num_sibs+1)/2.0);
2028	p2 = (long)((num_sibs+2)/2.0);
2029	i = 1;
2030	q = p->next;
2031	first = q->back;
2032	do {
2033	if (i == p1) pp1 = q->back;
2034	if (i == p2) pp2 = q->back;
2035	last = q->back;
2036	q = q->next;
2037	i++;
2038	} while (q != p);
2039	p->xcoord = (long)(0.5 - over * p->tyme);
2040	p->ycoord = (pp1->ycoord + pp2->ycoord) / 2;
2041	p->ymin = first->ymin;
2042	p->ymax = last->ymax;
2043	} /* promlk_coordinates */
2044
2045
2046	void promlk_drawline(long i, double scale)
2047	{
2048	/* draws one row of the tree diagram by moving up tree */
2049	node p, q, r, first =NULL, *last =NULL;
2050	long n, j;
2051	boolean extra, done;
2052
2053	p = curtree.root;
2054	q = curtree.root;
2055	extra = false;
2056	if ((long)(p->ycoord) == i) {
2057	if (p->index - spp >= 10)
2058	fprintf(outfile, "-%2ld", p->index - spp);
2059	else
2060	fprintf(outfile, "--%ld", p->index - spp);
2061	extra = true;
2062	} else
2063	fprintf(outfile, " ");
2064	do {
2065	if (!p->tip) {
2066	r = p->next;
2067	done = false;
2068	do {
2069	if (i >= r->back->ymin && i <= r->back->ymax) {
2070	q = r->back;
2071	done = true;
2072	}
2073	r = r->next;
2074	} while (!(done \|\| r == p));
2075	first = p->next->back;
2076	r = p->next;
2077	while (r->next != p)
2078	r = r->next;
2079	last = r->back;
2080	}
2081	done = (p == q);
2082	n = (long)(scale * ((long)(p->xcoord) - (long)(q->xcoord)) + 0.5);
2083	if (n < 3 && !q->tip)
2084	n = 3;
2085	if (extra) {
2086	n--;
2087	extra = false;
2088	}
2089	if ((long)(q->ycoord) == i && !done) {
2090	if (p->ycoord != q->ycoord)
2091	putc('+', outfile);
2092	else
2093	putc('-', outfile);
2094	if (!q->tip) {
2095	for (j = 1; j <= n - 2; j++)
2096	putc('-', outfile);
2097	if (q->index - spp >= 10)
2098	fprintf(outfile, "%2ld", q->index - spp);
2099	else
2100	fprintf(outfile, "-%ld", q->index - spp);
2101	extra = true;
2102	} else {
2103	for (j = 1; j < n; j++)
2104	putc('-', outfile);
2105	}
2106	} else if (!p->tip) {
2107	if ((long)(last->ycoord) > i && (long)(first->ycoord) < i &&
2108	i != (long)(p->ycoord)) {
2109	putc('!', outfile);
2110	for (j = 1; j < n; j++)
2111	putc(' ', outfile);
2112	} else {
2113	for (j = 1; j <= n; j++)
2114	putc(' ', outfile);
2115	}
2116	} else {
2117	for (j = 1; j <= n; j++)
2118	putc(' ', outfile);
2119	}
2120	if (p != q)
2121	p = q;
2122	} while (!done);
2123	if ((long)(p->ycoord) == i && p->tip) {
2124	for (j = 0; j < nmlngth; j++)
2125	putc(nayme[p->index - 1][j], outfile);
2126	}
2127	putc('\n', outfile);
2128	} /* promlk_drawline */
2129
2130
2131	void promlk_printree()
2132	{
2133	/* prints out diagram of the tree */
2134	long tipy;
2135	double scale;
2136	long i;
2137	node *p;
2138
2139	if (!treeprint)
2140	return;
2141	putc('\n', outfile);
2142	tipy = 1;
2143	promlk_coordinates(curtree.root, &tipy);
2144	p = curtree.root;
2145	while (!p->tip)
2146	p = p->next->back;
2147	scale = 1.0 / (long)(p->tyme - curtree.root->tyme + 1.000);
2148	putc('\n', outfile);
2149	for (i = 1; i <= tipy - down; i++)
2150	promlk_drawline(i, scale);
2151	putc('\n', outfile);
2152	} /* promlk_printree */
2153
2154
2155	void describe(node *p)
2156	{
2157	long i, num_sibs;
2158	node sib_ptr, sib_back_ptr;
2159	double v;
2160
2161	if (p == curtree.root)
2162	fprintf(outfile, " root ");
2163	else
2164	fprintf(outfile, "%4ld ", p->back->index - spp);
2165	if (p->tip) {
2166	for (i = 0; i < nmlngth; i++)
2167	putc(nayme[p->index - 1][i], outfile);
2168	} else
2169	fprintf(outfile, "%4ld ", p->index - spp);
2170	if (p != curtree.root) {
2171	fprintf(outfile, "%11.5f", (p->tyme - curtree.root->tyme));
2172	v = (p->tyme - curtree.nodep[p->back->index - 1]->tyme);
2173	fprintf(outfile, "%13.5f", v);
2174	}
2175	putc('\n', outfile);
2176	if (!p->tip) {
2177
2178	sib_ptr = p;
2179	num_sibs = count_sibs(p);
2180	for (i=0 ; i < num_sibs; i++) {
2181	sib_ptr = sib_ptr->next;
2182	sib_back_ptr = sib_ptr->back;
2183	describe(sib_back_ptr);
2184	}
2185	}
2186	} /* describe */
2187
2188
2189	void prot_reconstr(node *p, long n)
2190	{
2191	/* reconstruct and print out acid at site n+1 at node p */
2192	long i, j, k, first, num_sibs = 0;
2193	double f, sum, xx[20];
2194	node *q = NULL;
2195
2196	if (p->tip)
2197	putc(y[p->index-1][n], outfile);
2198	else {
2199	num_sibs = count_sibs(p);
2200	if ((ally[n] == 0) \|\| (location[ally[n]-1] == 0))
2201	putc('.', outfile);
2202	else {
2203	j = location[ally[n]-1] - 1;
2204	sum = 0;
2205	for (i = 0; i <= 19; i++) {
2206	f = p->protx[j][mx-1][i];
2207	if (!p->tip) {
2208	q = p;
2209	for (k = 0; k < num_sibs; k++) {
2210	q = q->next;
2211	f *= q->protx[j][mx-1][i];
2212	}
2213	}
2214	f = sqrt(f);
2215	xx[i] = f * freqaa[i];
2216	sum += xx[i];
2217	}
2218	for (i = 0; i <= 19; i++)
2219	xx[i] /= sum;
2220	first = 0;
2221	for (i = 0; i <= 19; i++)
2222	if (xx[i] > xx[first])
2223	first = i;
2224	if (xx[first] > 0.95)
2225	putc(aachar[first], outfile);
2226	else
2227	putc(tolower(aachar[first]), outfile);
2228	if (rctgry && rcategs > 1)
2229	mx = mp[n][mx - 1];
2230	else
2231	mx = 1;
2232	}
2233	}
2234	} /* prot_reconstr */
2235
2236
2237	void rectrav(node *p, long m, long n)
2238	{
2239	/* print out segment of reconstructed sequence for one branch */
2240	long num_sibs, i;
2241	node *sib_ptr;
2242
2243	putc(' ', outfile);
2244	if (p->tip) {
2245	for (i = 0; i < nmlngth; i++)
2246	putc(nayme[p->index-1][i], outfile);
2247	} else
2248	fprintf(outfile, "%4ld ", p->index - spp);
2249	fprintf(outfile, " ");
2250	mx = mx0;
2251	for (i = m; i <= n; i++) {
2252	if ((i % 10 == 0) && (i != m))
2253	putc(' ', outfile);
2254	prot_reconstr(p, i);
2255	}
2256	putc('\n', outfile);
2257	if (!p->tip) {
2258	num_sibs = count_sibs(p);
2259	sib_ptr = p;
2260	for (i = 0; i < num_sibs; i++) {
2261	sib_ptr = sib_ptr->next;
2262	rectrav(sib_ptr->back, m, n);
2263	}
2264	}
2265	mx1 = mx;
2266	} /* rectrav */
2267
2268
2269	void summarize()
2270	{
2271	long i, j, mm;
2272	double mode, sum;
2273	double like[maxcategs], nulike[maxcategs];
2274	double **marginal;
2275	long **mp;
2276
2277	mp = (long *)Malloc(sites sizeof(long *));
2278	for (i = 0; i <= sites-1; ++i)
2279	mp[i] = (long )Malloc(sizeof(long)rcategs);
2280	fprintf(outfile, "\nLn Likelihood = %11.5f\n\n", curtree.likelihood);
2281	fprintf(outfile, " Ancestor Node Node Height Length\n");
2282	fprintf(outfile, " -------- ---- ---- ------ ------\n");
2283	describe(curtree.root);
2284	putc('\n', outfile);
2285	if (rctgry && rcategs > 1) {
2286	for (i = 0; i < rcategs; i++)
2287	like[i] = 1.0;
2288	for (i = sites - 1; i >= 0; i--) {
2289	sum = 0.0;
2290	for (j = 0; j < rcategs; j++) {
2291	nulike[j] = (lambda1 + lambda * probcat[j]) * like[j];
2292	mp[i][j] = j + 1;
2293	for (k = 1; k <= rcategs; k++) {
2294	if (k != j + 1) {
2295	if (lambda * probcat[k - 1] * like[k - 1] > nulike[j]) {
2296	nulike[j] = lambda * probcat[k - 1] * like[k - 1];
2297	mp[i][j] = k;
2298	}
2299	}
2300	}
2301	if ((ally[i] > 0) && (location[ally[i]-1] > 0))
2302	nulike[j] *= contribution[location[ally[i] - 1] - 1][j];
2303	sum += nulike[j];
2304	}
2305	for (j = 0; j < rcategs; j++)
2306	nulike[j] /= sum;
2307	memcpy(like, nulike, rcategs * sizeof(double));
2308	}
2309	mode = 0.0;
2310	mx = 1;
2311	for (i = 1; i <= rcategs; i++) {
2312	if (probcat[i - 1] * like[i - 1] > mode) {
2313	mx = i;
2314	mode = probcat[i - 1] * like[i - 1];
2315	}
2316	}
2317	mx0 = mx;
2318	fprintf(outfile,
2319	"Combination of categories that contributes the most to the likelihood:\n\n");
2320	for (i = 1; i <= nmlngth + 3; i++)
2321	putc(' ', outfile);
2322	for (i = 1; i <= sites; i++) {
2323	fprintf(outfile, "%ld", mx);
2324	if (i % 10 == 0)
2325	putc(' ', outfile);
2326	if (i % 60 == 0 && i != sites) {
2327	putc('\n', outfile);
2328	for (j = 1; j <= nmlngth + 3; j++)
2329	putc(' ', outfile);
2330	}
2331	mx = mp[i - 1][mx - 1];
2332	}
2333	fprintf(outfile, "\n\n");
2334	marginal = (double *) Malloc( sitessizeof(double *));
2335	for (i = 0; i < sites; i++)
2336	marginal[i] = (double ) Malloc( rcategssizeof(double));
2337	for (i = 0; i < rcategs; i++)
2338	like[i] = 1.0;
2339	for (i = sites - 1; i >= 0; i--) {
2340	sum = 0.0;
2341	for (j = 0; j < rcategs; j++) {
2342	nulike[j] = (lambda1 + lambda * probcat[j]) * like[j];
2343	for (k = 1; k <= rcategs; k++) {
2344	if (k != j + 1)
2345	nulike[j] += lambda * probcat[k - 1] * like[k - 1];
2346	}
2347	if ((ally[i] > 0) && (location[ally[i]-1] > 0))
2348	nulike[j] *= contribution[location[ally[i] - 1] - 1][j];
2349	sum += nulike[j];
2350	}
2351	for (j = 0; j < rcategs; j++) {
2352	nulike[j] /= sum;
2353	marginal[i][j] = nulike[j];
2354	}
2355	memcpy(like, nulike, rcategs * sizeof(double));
2356	}
2357	for (i = 0; i < rcategs; i++)
2358	like[i] = 1.0;
2359	for (i = 0; i < sites; i++) {
2360	sum = 0.0;
2361	for (j = 0; j < rcategs; j++) {
2362	nulike[j] = (lambda1 + lambda * probcat[j]) * like[j];
2363	for (k = 1; k <= rcategs; k++) {
2364	if (k != j + 1)
2365	nulike[j] += lambda * probcat[k - 1] * like[k - 1];
2366	}
2367	marginal[i][j] = like[j] probcat[j];
2368	sum += nulike[j];
2369	}
2370	for (j = 0; j < rcategs; j++)
2371	nulike[j] /= sum;
2372	memcpy(like, nulike, rcategs * sizeof(double));
2373	sum = 0.0;
2374	for (j = 0; j < rcategs; j++)
2375	sum += marginal[i][j];
2376	for (j = 0; j < rcategs; j++)
2377	marginal[i][j] /= sum;
2378	}
2379	fprintf(outfile, "Most probable category at each site if > 0.95");
2380	fprintf(outfile, " probability (\".\" otherwise)\n\n");
2381	for (i = 1; i <= nmlngth + 3; i++)
2382	putc(' ', outfile);
2383	for (i = 0; i < sites; i++) {
2384	sum = 0.0;
2385	for (j = 0; j < rcategs; j++)
2386	if (marginal[i][j] > sum) {
2387	sum = marginal[i][j];
2388	mm = j;
2389	}
2390	if (sum >= 0.95)
2391	fprintf(outfile, "%ld", mm+1);
2392	else
2393	putc('.', outfile);
2394	if ((i+1) % 60 == 0) {
2395	if (i != 0) {
2396	putc('\n', outfile);
2397	for (j = 1; j <= nmlngth + 3; j++)
2398	putc(' ', outfile);
2399	}
2400	}
2401	else if ((i+1) % 10 == 0)
2402	putc(' ', outfile);
2403	}
2404	putc('\n', outfile);
2405	for (i = 0; i < sites; i++)
2406	free(marginal[i]);
2407	free(marginal);
2408	}
2409	putc('\n', outfile);
2410	putc('\n', outfile);
2411	putc('\n', outfile);
2412	if (hypstate) {
2413	fprintf(outfile, "Probable sequences at interior nodes:\n\n");
2414	fprintf(outfile, " node ");
2415	for (i = 0; (i < 13) && (i < ((sites + (sites-1)/10 - 39) / 2)); i++)
2416	putc(' ', outfile);
2417	fprintf(outfile, "Reconstructed sequence (caps if > 0.95)\n\n");
2418	if (!rctgry \|\| (rcategs == 1))
2419	mx0 = 1;
2420	for (i = 0; i < sites; i += 60) {
2421	k = i + 59;
2422	if (k >= sites)
2423	k = sites - 1;
2424	rectrav(curtree.root, i, k);
2425	putc('\n', outfile);
2426	mx0 = mx1;
2427	}
2428	}
2429	} /* summarize */
2430
2431
2432	void promlk_treeout(node *p)
2433	{
2434	/* write out file with representation of final tree */
2435	node *sib_ptr;
2436	long i, n, w, num_sibs;
2437	Char c;
2438	double x;
2439
2440	if (p->tip) {
2441	n = 0;
2442	for (i = 1; i <= nmlngth; i++) {
2443	if (nayme[p->index - 1][i - 1] != ' ')
2444	n = i;
2445	}
2446	for (i = 0; i < n; i++) {
2447	c = nayme[p->index - 1][i];
2448	if (c == ' ')
2449	c = '_';
2450	putc(c, outtree);
2451	}
2452	col += n;
2453	} else {
2454	sib_ptr = p;
2455	num_sibs = count_sibs(p);
2456	putc('(', outtree);
2457	col++;
2458
2459	for (i=0; i < (num_sibs - 1); i++) {
2460	sib_ptr = sib_ptr->next;
2461	promlk_treeout(sib_ptr->back);
2462	putc(',', outtree);
2463	col++;
2464	if (col > 55) {
2465	putc('\n', outtree);
2466	col = 0;
2467	}
2468	}
2469	sib_ptr = sib_ptr->next;
2470	promlk_treeout(sib_ptr->back);
2471	putc(')', outtree);
2472	col++;
2473	}
2474	if (p == curtree.root) {
2475	fprintf(outtree, ";\n");
2476	return;
2477	}
2478	x = (p->tyme - curtree.nodep[p->back->index - 1]->tyme);
2479	if (x > 0.0)
2480	w = (long)(0.4342944822 * log(x));
2481	else if (x == 0.0)
2482	w = 0;
2483	else
2484	w = (long)(0.4342944822 * log(-x)) + 1;
2485	if (w < 0)
2486	w = 0;
2487	fprintf(outtree, ":%*.5f", (int)(w + 7), x);
2488	col += w + 8;
2489	} /* promlk_treeout */
2490
2491
2492	void initpromlnode(node p, node grbg, node *q, long len, long nodei,
2493	long ntips, long parens, initops whichinit,
2494	pointarray treenode, pointarray nodep, Char *str,
2495	Char ch, FILE intree)
2496	{
2497	/* initializes a node */
2498	boolean minusread;
2499	double valyew, divisor;
2500
2501	switch (whichinit) {
2502	case bottom:
2503	gnu(grbg, p);
2504	(*p)->index = nodei;
2505	(*p)->tip = false;
2506	malloc_ppheno((*p), endsite, rcategs);
2507	nodep[(p)->index - 1] = (p);
2508	break;
2509	case nonbottom:
2510	gnu(grbg, p);
2511	malloc_ppheno(*p, endsite, rcategs);
2512	(*p)->index = nodei;
2513	break;
2514	case tip:
2515	match_names_to_data(str, nodep, p, spp);
2516	break;
2517	case iter:
2518	(*p)->initialized = false;
2519	(*p)->v = initialv;
2520	(*p)->iter = true;
2521	if ((*p)->back != NULL)
2522	(*p)->back->iter = true;
2523	break;
2524	case length:
2525	processlength(&valyew, &divisor, ch, &minusread, intree, parens);
2526	(*p)->v = valyew / divisor;
2527	(*p)->iter = false;
2528	if ((*p)->back != NULL) {
2529	(p)->back->v = (p)->v;
2530	(*p)->back->iter = false;
2531	}
2532	break;
2533	case unittrwt:
2534	curtree.nodep[spp]->iter = false;
2535	break;
2536	default: /* cases hslength, hsnolength, treewt */
2537	break; /* should never occur */
2538	}
2539	} /* initpromlnode */
2540
2541
2542	void tymetrav(node p, double x)
2543	{
2544	/* set up times of nodes */
2545	node sib_ptr, q;
2546	long i, num_sibs;
2547	double xmax;
2548
2549	xmax = 0.0;
2550	if (!p->tip) {
2551	sib_ptr = p;
2552	num_sibs = count_sibs(p);
2553	for (i=0; i < num_sibs; i++) {
2554	sib_ptr = sib_ptr->next;
2555	tymetrav(sib_ptr->back, x);
2556	if (xmax > (*x))
2557	xmax = (*x);
2558	}
2559	} else
2560	(*x) = 0.0;
2561	p->tyme = xmax;
2562	if (!p->tip) {
2563	q = p;
2564	while (q->next != p) {
2565	q = q->next;
2566	q->tyme = p->tyme;
2567	}
2568	}
2569	(*x) = p->tyme - p->v;
2570	} /* tymetrav */
2571
2572
2573	void free_all_protx (long nonodes, pointarray treenode)
2574	{
2575	/* used in proml */
2576	long i, j, k;
2577	node *p;
2578
2579	/* Zero thru spp are tips, */
2580	for (i = 0; i < spp; i++) {
2581	for (j = 0; j < endsite; j++)
2582	free(treenode[i]->protx[j]);
2583	free(treenode[i]->protx);
2584	}
2585
2586	/* The rest are rings (i.e. triads) */
2587	for (i = spp; i < nonodes; i++) {
2588	if (treenode[i] != NULL) {
2589	p = treenode[i];
2590	for (j = 1; j <= 3; j++) {
2591	for (k = 0; k < endsite; k++)
2592	free(p->protx[k]);
2593	free(p->protx);
2594	p = p->next;
2595	}
2596	}
2597	}
2598	} /* free_all_protx */
2599
2600
2601	void maketree()
2602	{
2603	/* constructs a binary tree from the pointers in curtree.nodep,
2604	adds each node at location which yields highest likelihood
2605	then rearranges the tree for greatest likelihood */
2606
2607	long i, j, k, l;
2608	long numtrees = 0, num_sibs;
2609	double bestlike, gotlike, x;
2610	node item, nufork, dummy, q, *root=NULL;
2611	boolean dummy_haslengths, dummy_first, goteof;
2612	long nextnode;
2613	long grcategs;
2614	pointarray dummy_treenode=NULL;
2615
2616	grcategs = (categs > rcategs) ? categs : rcategs;
2617
2618	prot_inittable();
2619
2620	if (!usertree) {
2621	for (i = 1; i <= spp; i++)
2622	enterorder[i - 1] = i;
2623	if (jumble)
2624	randumize(seed, enterorder);
2625	curtree.root = curtree.nodep[spp];
2626	curtree.root->back = NULL;
2627	for (i = 0; i < spp; i++)
2628	curtree.nodep[i]->back = NULL;
2629	for (i = spp; i < nonodes; i++) {
2630	q = curtree.nodep[i];
2631	q->back = NULL;
2632	while ((q = q->next) != curtree.nodep[i])
2633	q->back = NULL;
2634	}
2635	polishing = false;
2636	promlk_add(curtree.nodep[enterorder[0]-1], curtree.nodep[enterorder[1]-1],
2637	curtree.nodep[spp], false);
2638	if (progress) {
2639	printf("\nAdding species:\n");
2640	writename(0, 2, enterorder);
2641	#ifdef WIN32
2642	phyFillScreenColor();
2643	#endif
2644	}
2645	lastsp = false;
2646	smoothit = false;
2647	for (i = 3; i <= spp; i++) {
2648	bestyet = UNDEFINED;
2649	bestree.likelihood = bestyet;
2650	there = curtree.root;
2651	item = curtree.nodep[enterorder[i - 1] - 1];
2652	nufork = curtree.nodep[spp + i - 2];
2653	lastsp = (i == spp);
2654	addpreorder(curtree.root, item, nufork, true, true);
2655	promlk_add(there, item, nufork, false);
2656	like = prot_evaluate(curtree.root);
2657	rearrange(&curtree.root);
2658	if (curtree.likelihood > bestree.likelihood) {
2659	prot_copy_(&curtree, &bestree, nonodes, grcategs);
2660	}
2661	if (progress) {
2662	writename(i - 1, 1, enterorder);
2663	#ifdef WIN32
2664	phyFillScreenColor();
2665	#endif
2666	}
2667	if (lastsp && global) {
2668	if (progress) {
2669	printf("Doing global rearrangements\n");
2670	printf(" !");
2671	for (j = 1; j <= nonodes; j++)
2672	putchar('-');
2673	printf("!\n");
2674	}
2675	bestlike = bestyet;
2676	do {
2677	if (progress)
2678	printf(" ");
2679	gotlike = bestlike;
2680	for (j = 0; j < nonodes; j++) {
2681	bestyet = UNDEFINED;
2682	item = curtree.nodep[j];
2683	if (item != curtree.root) {
2684	nufork = curtree.nodep[curtree.nodep[j]->back->index - 1];
2685	promlk_re_move(&item, &nufork, false);
2686	there = curtree.root;
2687	addpreorder(curtree.root, item, nufork, true, true);
2688	promlk_add(there, item, nufork, false);
2689	}
2690	if (progress) {
2691	putchar('.');
2692	fflush(stdout);
2693	}
2694	}
2695	if (progress)
2696	putchar('\n');
2697	} while (bestlike < gotlike);
2698	}
2699	}
2700	if (njumble > 1 && lastsp) {
2701	for (i = 0; i < spp; i++ )
2702	promlk_re_move(&curtree.nodep[i], &dummy, false);
2703	if (jumb == 1 \|\| bestree2.likelihood < bestree.likelihood)
2704	prot_copy_(&bestree, &bestree2, nonodes, grcategs);
2705	}
2706	if (jumb == njumble) {
2707	if (njumble > 1)
2708	prot_copy_(&bestree2, &curtree, nonodes, grcategs);
2709	else
2710	prot_copy_(&bestree, &curtree, nonodes, grcategs);
2711	fprintf(outfile, "\n\n");
2712	treevaluate();
2713	curtree.likelihood = prot_evaluate(curtree.root);
2714	promlk_printree();
2715	summarize();
2716	if (trout) {
2717	col = 0;
2718	promlk_treeout(curtree.root);
2719	}
2720	}
2721	} else {
2722	openfile(&intree, INTREE, "input tree file", "r", progname, intreename);
2723	numtrees = countsemic(&intree);
2724	if (numtrees > 2)
2725	initseed(&inseed, &inseed0, seed);
2726	l0gl = (double )Malloc(numtrees sizeof(double));
2727	l0gf = (double *)Malloc(numtrees sizeof(double *));
2728	for (i=0;i<numtrees;++i)
2729	l0gf[i] = (double )Malloc(endsite sizeof(double));
2730	if (treeprint) {
2731	fprintf(outfile, "User-defined tree");
2732	if (numtrees > 1)
2733	putc('s', outfile);
2734	fprintf(outfile, ":\n\n");
2735	}
2736	fprintf(outfile, "\n\n");
2737	which = 1;
2738	while (which <= numtrees) {
2739
2740	/* These initializations required each time through the loop
2741	since multiple trees require re-initialization */
2742	dummy_haslengths = true;
2743	nextnode = 0;
2744	dummy_first = true;
2745	goteof = false;
2746
2747	treeread(intree, &root, dummy_treenode, &goteof, &dummy_first,
2748	curtree.nodep, &nextnode,
2749	&dummy_haslengths, &grbg, initpromlnode);
2750
2751	nonodes = nextnode;
2752
2753	root = curtree.nodep[root->index - 1];
2754	curtree.root = root;
2755
2756	if (lngths)
2757	tymetrav(curtree.root, &x);
2758
2759	if (goteof && (which <= numtrees)) {
2760	/* if we hit the end of the file prematurely */
2761	printf ("\n");
2762	printf ("ERROR: trees missing at end of file.\n");
2763	printf ("\tExpected number of trees:\t\t%ld\n", numtrees);
2764	printf ("\tNumber of trees actually in file:\t%ld.\n\n", which - 1);
2765	exxit(-1);
2766	}
2767	curtree.start = curtree.nodep[0]->back;
2768	treevaluate();
2769	promlk_printree();
2770	summarize();
2771	if (trout) {
2772	col = 0;
2773	promlk_treeout(curtree.root);
2774	}
2775	which++;
2776	}
2777
2778	FClose(intree);
2779	if (!auto_ && numtrees > 1 && weightsum > 1 )
2780	standev2(numtrees, maxwhich, 0, endsite, maxlogl, l0gl, l0gf,
2781	aliasweight, seed);
2782	}
2783	if (usertree) {
2784	free(l0gl);
2785	for (i=0; i < numtrees; i++)
2786	free(l0gf[i]);
2787	free(l0gf);
2788	}
2789	for (num_sibs = 0; num_sibs < max_num_sibs; num_sibs++) {
2790	for (j = 0; j < rcategs; j++) {
2791	for (k = 0; k < categs; k++) {
2792	for (l = 0; l < 20; l++) {
2793	free(pmatrices[num_sibs][j][k][l]);
2794	}
2795	free(pmatrices[num_sibs][j][k]);
2796	}
2797	free(pmatrices[num_sibs][j]);
2798	}
2799	free(pmatrices[num_sibs]);
2800	}
2801	if (jumb < njumble)
2802	return;
2803	free(contribution);
2804	free(mp);
2805	free_all_protx(nonodes2, curtree.nodep);
2806	if (!usertree \|\| reconsider) {
2807	free_all_protx(nonodes2, bestree.nodep);
2808	if (njumble > 1)
2809	free_all_protx(nonodes2, bestree2.nodep);
2810	}
2811	if (progress) {
2812	printf("\n\nOutput written to file \"%s\"\n\n", outfilename);
2813	if (trout)
2814	printf("Tree also written onto file \"%s\"\n", outtreename);
2815	putchar('\n');
2816	}
2817
2818	free(root);
2819	} /* maketree */
2820
2821
2822	void clean_up()
2823	{
2824	/* Free and/or close stuff */
2825	long i;
2826
2827	free (rrate);
2828	free (probcat);
2829	free (rate);
2830	/* Seems to require freeing every time... */
2831	for (i = 0; i < spp; i++) {
2832	free (y[i]);
2833	}
2834	free (y);
2835	free (nayme);
2836	free (enterorder);
2837	free (category);
2838	free (weight);
2839	free (alias);
2840	free (ally);
2841	free (location);
2842	free (aliasweight);
2843	free (probmat);
2844	free (eigmat);
2845	if (! (njumble <= 1))
2846	freetree2(bestree2.nodep, nonodes2);
2847	FClose(infile);
2848	FClose(outfile);
2849	FClose(outtree);
2850	#ifdef MAC
2851	fixmacfile(outfilename);
2852	fixmacfile(outtreename);
2853	#endif
2854	} /* clean_up */
2855
2856
2857	int main(int argc, Char *argv[])
2858	{ /* Protein Maximum Likelihood with molecular clock */
2859
2860	#ifdef MAC
2861	argc = 1; /* macsetup("Promlk", ""); */
2862	argv[0] = "Promlk";
2863	#endif
2864	init(argc,argv);
2865	progname = argv[0];
2866	openfile(&infile, INFILE, "input file", "r", argv[0], infilename);
2867	openfile(&outfile, OUTFILE, "output file", "w", argv[0], outfilename);
2868
2869	ibmpc = IBMCRT;
2870	ansi = ANSICRT;
2871	datasets = 1;
2872	mulsets = false;
2873	firstset = true;
2874	doinit();
2875
2876	if (trout)
2877	openfile(&outtree,OUTTREE,"output tree file","w",argv[0],outtreename);
2878	if (ctgry)
2879	openfile(&catfile,CATFILE,"categories file","r",argv[0],catfilename);
2880	if (weights \|\| justwts)
2881	openfile(&weightfile,WEIGHTFILE,"weights file","r",argv[0],weightfilename);
2882	for (ith = 1; ith <= datasets; ith++) {
2883	if (datasets > 1) {
2884	fprintf(outfile, "Data set # %ld:\n\n", ith);
2885	if (progress)
2886	printf("\nData set # %ld:\n", ith);
2887	}
2888	getinput();
2889
2890	if (ith == 1)
2891	firstset = false;
2892	for (jumb = 1; jumb <= njumble; jumb++){
2893	max_num_sibs = 0;
2894	maketree();
2895	}
2896	}
2897
2898	clean_up();
2899	printf("Done.\n\n");
2900	#ifdef WIN32
2901	phyRestoreConsoleAttributes();
2902	#endif
2903	return 0;
2904	} /* Protein Maximum Likelihood with molecular clock */
2905

Note: See TracBrowser for help on using the repository browser.

Download in other formats: