source: trunk/GDE/PHYLIP/proml.c

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