Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

promlk.c

Visit:

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

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

Context Navigation

source: trunk/GDE/PHYLIP/promlk.c

Download in other formats: