forester/archive/RIO/others/phylip_mod/src/proml.c

   1
   2 #include "phylip.h"
   3 #include "seq.h"
   4
   5 /* version 3.6. (c) Copyright 1993-2004 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
  12 typedef long vall[maxcategs];
  13 typedef double contribarr[maxcategs];
  14
  15 #ifndef OLDC
  16 /* function prototypes */
  17 void   init_protmats(void);
  18 void   getoptions(void);
  19 void   makeprotfreqs(void);
  20 void   allocrest(void);
  21 void   doinit(void);
  22 void   inputoptions(void);
  23 void   input_protdata(long);
  24 void   makeweights(void);
  25 void   prot_makevalues(long, pointarray, long, long, sequence, steptr);
  26 void   prot_inittable(void);
  27
  28 void   alloc_pmatrix(long);
  29 void   getinput(void);
  30 void   inittravtree(node *);
  31 void   prot_nuview(node *);
  32 void   prot_slopecurv(node *, double, double *, double *, double *);
  33 void   makenewv(node *);
  34 void   update(node *);
  35 void   smooth(node *);
  36 void   make_pmatrix(double **, double **, double **, long, double,
  37             double, double *, double **);
  38 double prot_evaluate(node *, boolean);
  39
  40 void   treevaluate(void);
  41 void   promlcopy(tree *, tree *, long, long);
  42 void   proml_re_move(node **, node **);
  43 void   insert_(node *, node *, boolean);
  44 void   addtraverse(node *, node *, boolean);
  45 void   rearrange(node *, node *);
  46 void   proml_coordinates(node *, double, long *, double *);
  47 void   proml_printree(void);
  48 void   sigma(node *, double *, double *, double *);
  49 void   describe(node *);
  50
  51 void   prot_reconstr(node *, long);
  52 void   rectrav(node *, long, long);
  53 void   summarize(void);
  54 void   initpromlnode(node **, node **, node *, long, long, long *, long *,
  55             initops, pointarray, pointarray, Char *, Char *, FILE *);
  56 void   dnaml_treeout(node *);
  57 void   buildnewtip(long, tree *);
  58 void   buildsimpletree(tree *);
  59 void   free_all_protx (long, pointarray);
  60 void   maketree(void);
  61 void   clean_up(void);
  62 void   globrearrange(void);
  63 void   proml_unroot(node* root, node** nodep, long nonodes) ;
  64 void   reallocsites(void);
  65 void   prot_freetable(void);
  66 void   free_pmatrix(long sib);
  67 void   alloclrsaves(void);
  68 void   freelrsaves(void);
  69 void   resetlrsaves(void);
  70 /* function prototypes */
  71 #endif
  72
  73
  74 long rcategs;
  75 boolean haslengths;
  76 long oldendsite=0;
  77
  78 Char infilename[100], outfilename[100], intreename[100], outtreename[100],
  79      catfilename[100], weightfilename[100];
  80 double *rate, *rrate, *probcat;
  81 long nonodes2, sites, weightsum, categs,
  82   datasets, ith, njumble, jumb;
  83 long inseed, inseed0, parens;
  84 boolean global, jumble, weights, trout, usertree, inserting = false,
  85   ctgry, rctgry, auto_, hypstate, progress, mulsets, justwts, firstset,
  86   improve, smoothit, polishing, lngths, gama, invar, usepmb, usepam, usejtt;
  87 tree curtree, bestree, bestree2, priortree;
  88 node *qwhere, *grbg, *addwhere;
  89 double cv, alpha, lambda, invarfrac, bestyet;
  90 long *enterorder;
  91 steptr aliasweight;
  92 contribarr *contribution, like, nulike, clai;
  93 double **term, **slopeterm, **curveterm;
  94 longer seed;
  95 char *progname;
  96 char aachar[26]="ARNDCQEGHILKMFPSTWYVBZX?*-";
  97 node **lrsaves;
  98
  99 /* Local variables for maketree, propagated globally for c version: */
 100 long k, nextsp, numtrees, maxwhich, mx, mx0, mx1, shimotrees;
 101 double dummy, maxlogl;
 102 boolean succeeded, smoothed;
 103 double **l0gf;
 104 double *l0gl;
 105 double **tbl;
 106 Char ch, ch2;
 107 long col;
 108 vall *mp;
 109
 110
 111 /* Variables introduced to allow for protein probability calculations   */
 112 long   max_num_sibs;            /* maximum number of siblings used in a */
 113                                 /* nuview calculation.  determines size */
 114                                 /* final size of pmatrices              */
 115 double *eigmat;                 /* eig matrix variable                  */
 116 double **probmat;               /* prob matrix variable                 */
 117 double ****dpmatrix;            /* derivative of pmatrix                */
 118 double ****ddpmatrix;           /* derivative of xpmatrix               */
 119 double *****pmatrices;          /* matrix of probabilities of protien   */
 120                                 /* conversion.  The 5 subscripts refer  */
 121                                 /* to sibs, rcategs, categs, final and  */
 122                                 /* initial states, respectively.        */
 123 double freqaa[20];              /* amino acid frequencies               */
 124
 125 /* this JTT matrix decomposition thanks to Elisabeth Tillier */
 126 static double jtteigmat[] =
 127 {0.0,        -0.7031123, -0.6484345, -0.6086499, -0.5514432,
 128 -0.772664, -0.8643413, -1.0620756, -0.9965552, -1.1671808,
 129 -1.2222418,-0.4589201, -1.3103714, -1.4048038, -0.3170582,
 130 -0.347935, -1.5311677, -1.6021194, -1.7991454, -1.8911888};
 131
 132 static double jttprobmat[20][20] =
 133 {{0.076999996, 0.051000003, 0.043000004, 0.051999998, 0.019999996, 0.041,
 134   0.061999994, 0.073999997, 0.022999999, 0.052000004, 0.090999997, 0.058999988,
 135   0.024000007, 0.04, 0.050999992, 0.069, 0.059000006, 0.014000008, 0.032000004,
 136   0.066000005},
 137  {0.015604455, -0.068062363, 0.020106264, 0.070723273, 0.011702977, 0.009674053,
 138   0.074000798, -0.169750458, 0.005560808, -0.008208636, -0.012305869,
 139  -0.063730179, -0.005674643, -0.02116828, 0.104586169, 0.016480839, 0.016765139,
 140   0.005936994, 0.006046367, -0.0082877},
 141  {-0.049778281, -0.007118197, 0.003801272, 0.070749616, 0.047506147,
 142    0.006447017, 0.090522425, -0.053620432, -0.008508175, 0.037170603,
 143    0.051805545, 0.015413608, 0.019939916, -0.008431976, -0.143511376,
 144   -0.052486072, -0.032116542, -0.000860626, -0.02535993, 0.03843545},
 145  {-0.028906423, 0.092952047, -0.009615343, -0.067870117, 0.031970392,
 146    0.048338335, -0.054396304, -0.135916654, 0.017780083, 0.000129242,
 147    0.031267424, 0.116333586, 0.007499746, -0.032153596, 0.033517051,
 148   -0.013719269, -0.00347293, -0.003291821, -0.02158326, -0.008862168},
 149  {0.037181176, -0.023106564, -0.004482225, -0.029899635, 0.118139633,
 150  -0.032298569, -0.04683198, 0.05566988, -0.012622847, 0.002023096,
 151  -0.043921088, -0.04792557, -0.003452711, -0.037744513, 0.020822974,
 152   0.036580187, 0.02331425, -0.004807711, -0.017504496, 0.01086673},
 153  {0.044754061, -0.002503471, 0.019452517, -0.015611487, -0.02152807,
 154  -0.013131425, -0.03465365, -0.047928912, 0.020608851, 0.067843095,
 155  -0.122130014, 0.002521499, 0.013021646, -0.082891087, -0.061590119,
 156   0.016270856, 0.051468938, 0.002079063, 0.081019713, 0.082927944},
 157  {0.058917882, 0.007320741, 0.025278141, 0.000357541, -0.002831285,
 158  -0.032453034, -0.010177288, -0.069447924, -0.034467324, 0.011422358,
 159  -0.128478324, 0.04309667, -0.015319944, 0.113302422, -0.035052393,
 160   0.046885372, 0.06185183, 0.00175743, -0.06224497, 0.020282093},
 161  {-0.014562092, 0.022522921, -0.007094389, 0.03480089, -0.000326144,
 162  -0.124039037, 0.020577906, -0.005056454, -0.081841576, -0.004381786,
 163   0.030826152, 0.091261631, 0.008878828, -0.02829487, 0.042718836,
 164  -0.011180886, -0.012719227, -0.000753926, 0.048062375, -0.009399129},
 165  {0.033789571, -0.013512235, 0.088010984, 0.017580292, -0.006608005,
 166  -0.037836971, -0.061344686, -0.034268357, 0.018190209, -0.068484614,
 167   0.120024744, -0.00319321, -0.001349477, -0.03000546, -0.073063759,
 168   0.081912399, 0.0635245, 0.000197, -0.002481798, -0.09108114},
 169  {-0.113947615, 0.019230545, 0.088819683, 0.064832765, 0.001801467,
 170  -0.063829682, -0.072001633, 0.018429333, 0.057465965, 0.043901014,
 171  -0.048050874, -0.001705918, 0.022637173, 0.017404665, 0.043877902,
 172  -0.017089594, -0.058489485, 0.000127498, -0.029357194, 0.025943972},
 173  {0.01512923, 0.023603725, 0.006681954, 0.012360216, -0.000181447,
 174  -0.023011838, -0.008960024, -0.008533239, 0.012569835, 0.03216118,
 175   0.061986403, -0.001919083, -0.1400832, -0.010669741, -0.003919454,
 176  -0.003707024, -0.026806029, -0.000611603, -0.001402648, 0.065312824},
 177  {-0.036405351, 0.020816769, 0.011408213, 0.019787053, 0.038897829,
 178    0.017641789, 0.020858533, -0.006067252, 0.028617353, -0.064259496,
 179   -0.081676567, 0.024421823, -0.028751676, 0.07095096, -0.024199434,
 180   -0.007513119, -0.028108766, -0.01198095, 0.111761119, -0.076198809},
 181  {0.060831772, 0.144097327, -0.069151377, 0.023754576, -0.003322955,
 182  -0.071618574, 0.03353154, -0.02795295, 0.039519769, -0.023453968,
 183  -0.000630308, -0.098024591, 0.017672997, 0.003813378, -0.009266499,
 184  -0.011192111, 0.016013873, -0.002072968, -0.010022044, -0.012526904},
 185  {-0.050776604, 0.092833081, 0.044069596, 0.050523021, -0.002628417,
 186    0.076542572, -0.06388631, -0.00854892, -0.084725311, 0.017401063,
 187   -0.006262541, -0.094457679, -0.002818678, -0.0044122, -0.002883973,
 188    0.028729685, -0.004961596, -0.001498627, 0.017994575, -0.000232779},
 189  {-0.01894566, -0.007760205, -0.015160993, -0.027254587, 0.009800903,
 190   -0.013443561, -0.032896517, -0.022734138, -0.001983861, 0.00256111,
 191    0.024823166, -0.021256768, 0.001980052, 0.028136263, -0.012364384,
 192   -0.013782446, -0.013061091, 0.111173981, 0.021702122, 0.00046654},
 193  {-0.009444193, -0.042106824, -0.02535015, -0.055125574, 0.006369612,
 194   -0.02945416, -0.069922064, -0.067221068, -0.003004999, 0.053624311,
 195    0.128862984, -0.057245803, 0.025550508, 0.087741073, -0.001119043,
 196   -0.012036202, -0.000913488, -0.034864475, 0.050124813, 0.055534723},
 197  {0.145782464, -0.024348311, -0.031216873, 0.106174443, 0.00202862,
 198   0.02653866, -0.113657267, -0.00755018, 0.000307232, -0.051241158,
 199   0.001310685, 0.035275877, 0.013308898, 0.002957626, -0.002925034,
 200  -0.065362319, -0.071844582, 0.000475894, -0.000112419, 0.034097762},
 201  {0.079840455, 0.018769331, 0.078685899, -0.084329807, -0.00277264,
 202  -0.010099754, 0.059700608, -0.019209715, -0.010442992, -0.042100476,
 203  -0.006020556, -0.023061786, 0.017246106, -0.001572858, -0.006703785,
 204   0.056301316, -0.156787357, -0.000303638, 0.001498195, 0.051363455},
 205  {0.049628261, 0.016475144, 0.094141653, -0.04444633, 0.005206131,
 206  -0.001827555, 0.02195624, 0.013066683, -0.010415582, -0.022338403,
 207   0.007837197, -0.023397671, -0.002507095, 0.005177694, 0.017109561,
 208  -0.202340113, 0.069681441, 0.000120736, 0.002201146, 0.004670849},
 209  {0.089153689, 0.000233354, 0.010826822, -0.004273519, 0.001440618,
 210   0.000436077, 0.001182351, -0.002255508, -0.000700465, 0.150589876,
 211  -0.003911914, -0.00050154, -0.004564983, 0.00012701, -0.001486973,
 212  -0.018902754, -0.054748555, 0.000217377, -0.000319302, -0.162541651}};
 213
 214 /* this PMB matrix decomposition due to Elisabeth Tillier */
 215 static double pmbeigmat[20] =
 216 {0.0000001586972220,-1.8416770496147100, -1.6025046986139100,-1.5801012515121300,
 217 -1.4987794099715900,-1.3520794233801900,-1.3003469390479700,-1.2439503327631300,
 218 -1.1962574080244200,-1.1383730501367500,-1.1153278910708000,-0.4934843510654760,
 219 -0.5419014550215590,-0.9657997830826700,-0.6276075673757390,-0.6675927795018510,
 220 -0.6932641383465870,-0.8897872681859630,-0.8382698977371710,-0.8074694642446040};
 221
 222 static double pmbprobmat[20][20] =
 223 {{0.0771762457248147,0.0531913844998640,0.0393445076407294,0.0466756566755510,
 224 0.0286348361997465,0.0312327748383639,0.0505410248721427,0.0767106611472993,
 225 0.0258916271688597,0.0673140562194124,0.0965705469252199,0.0515979465932174,
 226 0.0250628079438675,0.0503492018628350,0.0399908189418273,0.0641898881894471,
 227 0.0517539616710987,0.0143507440546115,0.0357994592438322,0.0736218495862984},
 228 {0.0368263046116572,-0.0006728917107827,0.0008590805287740,-0.0002764255356960,
 229 0.0020152937187455,0.0055743720652960,0.0003213317669367,0.0000449190281568,
 230 -0.0004226254397134,0.1805040629634510,-0.0272246813586204,0.0005904606533477,
 231 -0.0183743200073889,-0.0009194625608688,0.0008173657533167,-0.0262629806302238,
 232 0.0265738757209787,0.0002176606241904,0.0021315644838566,-0.1823229927207580},
 233 {-0.0194800075560895,0.0012068088610652,-0.0008803318319596,-0.0016044273960017,
 234 -0.0002938633803197,-0.0535796754602196,0.0155163896648621,-0.0015006360762140,
 235 0.0021601372013703,0.0268513218744797,-0.1085292493742730,0.0149753083138452,
 236 0.1346457366717310,-0.0009371698759829,0.0013501708044116,0.0346352293103622,
 237 -0.0276963770242276,0.0003643142783940,0.0002074817333067,-0.0174108903914110},
 238 {0.0557839400850153,0.0023271577185437,0.0183481103396687,0.0023339480096311,
 239 0.0002013267015151,-0.0227406863569852,0.0098644845475047,0.0064721276774396,
 240 0.0001389408104210,-0.0473713878768274,-0.0086984445005797,0.0026913674934634,
 241 0.0283724052562196,0.0001063665179457,0.0027442574779383,-0.1875312134708470,
 242 0.1279864877057640,0.0005103347834563,0.0003155113168637,0.0081451082759554},
 243 {0.0037510125027265,0.0107095920636885,0.0147305410328404,-0.0112351252180332,
 244 -0.0001500408626446,-0.1523450933729730,0.0611532413339872,-0.0005496748939503,
 245 0.0048714378736644,-0.0003826320053999,0.0552010244407311,0.0482555671001955,
 246 -0.0461664995115847,-0.0021165008617978,-0.0004574454232187,0.0233755883688949,
 247 -0.0035484915422384,0.0009090698422851,0.0013840637687758,-0.0073895139302231},
 248 {-0.0111512564930024,0.1025460064723080,0.0396772456883791,-0.0298408501361294,
 249 -0.0001656742634733,-0.0079876311843289,0.0712644184507945,-0.0010780604625230,
 250 -0.0035880882043592,0.0021070399334252,0.0016716329894279,-0.1810123023850110,
 251 0.0015141703608724,-0.0032700852781804,0.0035503782441679,0.0118634302028026,
 252 0.0044561606458028,-0.0001576678495964,0.0023470722225751,-0.0027457045397157},
 253 {0.1474525743949170,-0.0054432538500293,0.0853848892349828,-0.0137787746207348,
 254 -0.0008274830358513,0.0042248844582553,0.0019556229305563,-0.0164191435175148,
 255 -0.0024501858854849,0.0120908948084233,-0.0381456105972653,0.0101271614855119,
 256 -0.0061945941321859,0.0178841099895867,-0.0014577779202600,-0.0752120602555032,
 257 -0.1426985695849920,0.0002862275078983,-0.0081191734261838,0.0313401149422531},
 258 {0.0542034611735289,-0.0078763926211829,0.0060433542506096,0.0033396210615510,
 259 0.0013965072374079,0.0067798903832256,-0.0135291136622509,-0.0089982442731848,
 260 -0.0056744537593887,-0.0766524225176246,0.1881210263933930,-0.0065875518675173,
 261 0.0416627569300375,-0.0953804133524747,-0.0012559228448735,0.0101622644292547,
 262 -0.0304742453119050,0.0011702318499737,0.0454733434783982,-0.1119239362388150},
 263 {0.1069409037912470,0.0805064400880297,-0.1127352030714600,0.1001181253523260,
 264 -0.0021480427488769,-0.0332884841459003,-0.0679837575848452,-0.0043812841356657,
 265 0.0153418716846395,-0.0079441315103188,-0.0121766182046363,-0.0381127991037620,
 266 -0.0036338726532673,0.0195324059593791,-0.0020165963699984,-0.0061222685010268,
 267 -0.0253761448771437,-0.0005246410999057,-0.0112205170502433,0.0052248485517237},
 268 {-0.0325247648326262,0.0238753651653669,0.0203684886605797,0.0295666232678825,
 269 -0.0003946714764213,-0.0157242718469554,-0.0511737848084862,0.0084725632040180,
 270 -0.0167068828528921,0.0686962159427527,-0.0659702890616198,-0.0014289912494271,
 271 -0.0167000964093416,-0.1276689083678200,0.0036575057830967,-0.0205958145531018,
 272 0.0000368919612829,0.0014413626622426,0.1064360941926030,0.0863372661517408},
 273 {-0.0463777468104402,0.0394712148670596,0.1118686750747160,0.0440711686389031,
 274 -0.0026076286506751,-0.0268454015202516,-0.1464943067133240,-0.0137514051835380,
 275 -0.0094395514284145,-0.0144124844774228,0.0249103379323744,-0.0071832157138676,
 276 0.0035592787728526,0.0415627419826693,0.0027040097365669,0.0337523666612066,
 277 0.0316121324137152,-0.0011350177559026,-0.0349998884574440,-0.0302651879823361},
 278 {0.0142360925194728,0.0413145623127025,0.0324976427846929,0.0580930922002398,
 279 -0.0586974207121084,0.0202001168873069,0.0492204086749069,0.1126593173463060,
 280 0.0116620013776662,-0.0780333711712066,-0.1109786767320410,0.0407775100936731,
 281 -0.0205013161312652,-0.0653458585025237,0.0347351829703865,0.0304448983224773,
 282 0.0068813748197884,-0.0189002309261882,-0.0334507528405279,-0.0668143558699485},
 283 {-0.0131548829657936,0.0044244322828034,-0.0050639951827271,-0.0038668197633889,
 284 -0.1536822386530220,0.0026336969165336,0.0021585651200470,-0.0459233839062969,
 285 0.0046854727140565,0.0393815434593599,0.0619554007991097,0.0027456299925622,
 286 0.0117574347936383,0.0373018612990383,0.0024818527553328,-0.0133956606027299,
 287 -0.0020457128424105,0.0154178819990401,0.0246524142683911,0.0275363065682921},
 288 {-0.1542307272455030,0.0364861558267547,-0.0090880407008181,0.0531673937889863,
 289 0.0157585615170580,0.0029986538457297,0.0180194047699875,0.0652152443589317,
 290 0.0266842840376180,0.0388457366405908,0.0856237634510719,0.0126955778952183,
 291 0.0099593861698250,-0.0013941794862563,0.0294065511237513,-0.1151906949298290,
 292 -0.0852991447389655,0.0028699120202636,-0.0332087026659522,0.0006811857297899},
 293 {0.0281300736924501,-0.0584072081898638,-0.0178386569847853,-0.0536470338171487,
 294 -0.0186881656029960,-0.0240008730656106,-0.0541064820498883,0.2217137098936020,
 295 -0.0260500001542033,0.0234505236798375,0.0311127151218573,-0.0494139126682672,
 296 0.0057093465049849,0.0124937286655911,-0.0298322975915689,0.0006520211333102,
 297 -0.0061018680727128,-0.0007081999479528,-0.0060523759094034,0.0215845995364623},
 298 {0.0295321046399105,-0.0088296411830544,-0.0065057049917325,-0.0053478115612781,
 299 -0.0100646496794634,-0.0015473619084872,0.0008539960632865,-0.0376381933046211,
 300 -0.0328135588935604,0.0672161874239480,0.0667626853916552,-0.0026511651464901,
 301 0.0140451514222062,-0.0544836996133137,0.0427485157912094,0.0097455780205802,
 302 0.0177309072915667,-0.0828759701187452,-0.0729504795471370,0.0670731961252313},
 303 {0.0082646581043963,-0.0319918630534466,-0.0188454445200422,-0.0374976353856606,
 304 0.0037131290686848,-0.0132507796987883,-0.0306958830735725,-0.0044119395527308,
 305 -0.0140786756619672,-0.0180512599925078,-0.0208243802903953,-0.0232202769398931,
 306 -0.0063135878270273,0.0110442171178168,0.1824538048228460,-0.0006644614422758,
 307 -0.0069909097436659,0.0255407650654681,0.0099119399501151,-0.0140911517070698},
 308 {0.0261344441524861,-0.0714454044548650,0.0159436926233439,0.0028462736216688,
 309 -0.0044572637889080,-0.0089474834434532,-0.0177570282144517,-0.0153693244094452,
 310 0.1160919467206400,0.0304911481385036,0.0047047513411774,-0.0456535116423972,
 311 0.0004491494948617,-0.0767108879444462,-0.0012688533741441,0.0192445965934123,
 312 0.0202321954782039,0.0281039933233607,-0.0590403018490048,0.0364080426546883},
 313 {0.0115826306265004,0.1340228176509380,-0.0236200652949049,-0.1284484655137340,
 314 -0.0004742338006503,0.0127617346949511,-0.0428560878860394,0.0060030732454125,
 315 0.0089182609926781,0.0085353834972860,0.0048464809638033,0.0709740071429510,
 316 0.0029940462557054,-0.0483434904493132,-0.0071713680727884,-0.0036840391887209,
 317 0.0031454003250096,0.0246243550241551,-0.0449551277644180,0.0111449232769393},
 318 {0.0140356721886765,-0.0196518236826680,0.0030517022326582,0.0582672093364850,
 319 -0.0000973895685457,0.0021704767224292,0.0341806268602705,-0.0152035987563018,
 320 -0.0903198657739177,0.0259623214586925,0.0155832497882743,-0.0040543568451651,
 321 0.0036477631918247,-0.0532892744763217,-0.0142569373662724,0.0104500681408622,
 322 0.0103483945857315,0.0679534422398752,-0.0768068882938636,0.0280289727046158}}
 323 ;
 324
 325
 326 static double pameigmat[] = {0.0, -0.2350753691875762, -0.2701991863800379,
 327          -0.2931612442853115, -0.4262492032364507, -0.5395980482561625,
 328          -0.7141172690079523, -0.7392844756151318, -0.7781761342200766,
 329          -0.810032066366362, -0.875299712761124, -1.048227332164386,
 330          -1.109594097332267, -1.298616073142234, -1.342036228188581,
 331          -1.552599145527578, -1.658762802054814, -1.74893445623765,
 332          -1.933280832903272, -2.206353522613025};
 333
 334 static double pamprobmat[20][20] =
 335  {{0.087683339901135, 0.04051291829598762, 0.04087846315185977,
 336    0.04771603459744777, 0.03247095396561266, 0.03784612688594957,
 337    0.0504933695604875, 0.0898249006830755, 0.03285885059543713,
 338    0.0357514442352119, 0.0852464099207521, 0.07910313444070642,
 339    0.01488243946396588, 0.04100101908956829, 0.05158026947089499,
 340    0.06975497205982451, 0.05832757042475474, 0.00931264523877807,
 341    0.03171540880870517, 0.06303972920984541},
 342   {0.01943453646811026, -0.004492574160484092, 0.007694891061220776,
 343    0.01278399096887701, 0.0106157418450234, 0.007542140341575122,
 344    0.01326994069032819, 0.02615565199894889, 0.003123125764490066,
 345    0.002204507682495444, -0.004782898215768979, 0.01204241965177619,
 346    0.0007847400096924341, -0.03043626073172116, 0.01221202591902536,
 347    0.01100527004684405, 0.01116495631339549, -0.0925364931988571,
 348    -0.02622065387931562, 0.00843494142432107},
 349   {0.01855357100209072, 0.01493642835763868, 0.0127983090766285,
 350    0.0200533250704364, -0.1681898360107787, 0.01551657969909255,
 351    0.02128060163107209, 0.03100633591848964, 0.00845480845269879,
 352    0.000927149370785571, 0.00937207565817036, 0.03490557769673472,
 353    0.00300443019551563, -0.02590837220264415, 0.01329376859943192,
 354    0.006854110889741407, 0.01102593860528263, 0.003360844186685888,
 355    -0.03459712356647764, 0.003351477369404443},
 356   {0.02690642688200102, 0.02131745801890152, 0.0143626616005213,
 357    0.02405101425725929, 0.05041008641436849, 0.01430925051050233,
 358    0.02362114036816964, 0.04688381789373886, 0.005250115453626377,
 359    -0.02040112168595516, -0.0942720776915669, 0.03773004996758644,
 360    -0.00822831940782616, -0.1164872809439224, 0.02286281877257392,
 361    0.02849551240669926, 0.01468856796295663, 0.02377110964207936,
 362    -0.094380545436577, -0.02089068498518036},
 363   {0.00930172577225213, 0.01493463068441099, 0.020186920775608,
 364    0.02892154953912524, -0.01224593358361567, 0.01404228329986624,
 365    0.02671186617119041, 0.04537535161795231, 0.02229995804098249,
 366    -0.04635704133961575, -0.1966910360247138, 0.02796648065439046,
 367    -0.02263484732621436, 0.0440490503242072, 0.01148782948302166,
 368    0.01989170531824069, 0.001306805142981245, -0.005676690969116321,
 369    0.07680476281625202, -0.07967537039721849},
 370   {0.06602274245435476, -0.0966661981471856, -0.005241648783844579,
 371    0.00859135188171146, -0.007762129660943368, -0.02888965572526196,
 372    0.003592291525888222, 0.1668410669287673, -0.04082039290551406,
 373    0.005233775047553415, -0.01758244726137135, -0.1493955762326898,
 374    -0.00855819137835548, 0.004211419253492328, 0.01929306335052688,
 375    0.03008056746359405, 0.0190444422412472, 0.005577189741419315,
 376    0.0000874156155112068, 0.02634091459108298},
 377   {0.01933897472880726, 0.05874583569377844, -0.02293534606228405,
 378    -0.07206314017962175, -0.004580681581546643, -0.0628814337610561,
 379    -0.0850783812795136, 0.07988417636610614, -0.0852798990133397,
 380    0.01649047166155952, -0.05416647263757423, 0.1089834536254064,
 381    0.005093403979413865, 0.02520300254161142, 0.0005951431406455604,
 382    0.02441251821224675, 0.02796099482240553, -0.002574933994926502,
 383    -0.007172237553012804, 0.03002455129086954},
 384   {0.04041118479094272, -0.002476225672095412, -0.01494505811263243,
 385    -0.03759443758599911, -0.00892246902492875, -0.003634714029239211,
 386    -0.03085671837973749, -0.126176309029931, 0.005814031139083794,
 387    0.01313561962646063, -0.04760487162503322, -0.0490563712725484,
 388    -0.005082243450421558, -0.01213634309383557, 0.1806666927079249,
 389    0.02111663336185495, 0.02963486860587087, -0.0000175020101657785,
 390    0.01197155383597686, 0.0357526792184636},
 391   {-0.01184769557720525, 0.01582776076338872, -0.006570708266564639,
 392    -0.01471915653734024, 0.00894343616503608, 0.00562664968033149,
 393    -0.01465878888356943, 0.05365282692645818, 0.00893509735776116,
 394    -0.05879312944436473, 0.0806048683392995, -0.007722897986905326,
 395    -0.001819943882718859, 0.0942535573077267, 0.07483883782251654,
 396    0.004354639673913651, -0.02828804845740341, -0.001318222184691827,
 397    -0.07613149604246563, -0.1251675867732172},
 398   {0.00834167031558193, -0.01509357596974962, 0.007098172811092488,
 399    0.03127677418040319, 0.001992448468465455, 0.00915441566808454,
 400    0.03430175973499201, -0.0730648147535803, -0.001402707145575659,
 401    0.04780949194330815, -0.1115035603461273, -0.01292297197609604,
 402    -0.005056270550868528, 0.1112053349612027, -0.03801929822379964,
 403    -0.001191241001736563, 0.01872874622910247, 0.0005314214903865993,
 404    -0.0882576318311789, 0.07607183599610171},
 405   {-0.01539460099727769, 0.04988596184297883, -0.01187240760647617,
 406    -0.06987843637091853, -0.002490472846497859, 0.01009857892494956,
 407    -0.07473588067847209, 0.0906009925879084, 0.1243612446505172,
 408    0.02152806401345371, -0.03504879644860233, -0.06680752427613573,
 409    -0.005574485153629651, 0.001518282948127752, -0.01999168507510701,
 410    -0.01478606199529457, -0.02203749419458996, -0.00132680708294333,
 411    -0.01137505997867614, 0.05332658773667142},
 412   {-0.06104378736432388, 0.0869446603393548, -0.03298331234537257,
 413    0.03128515657456024, 0.003906358569208259, 0.03578694104193928,
 414    0.06241936133189683, 0.06182827284921748, -0.05566564263245907,
 415    0.02640868588189002, -0.01349751243059039, -0.05507866642582638,
 416    -0.006671347738489326, -0.001470096466016046, 0.05185743641479938,
 417    -0.07494697511168257, -0.1175185439057584, -0.001188074094105709,
 418    0.00937934805737347, 0.05024773745437657},
 419   {-0.07252555582124737, -0.116554459356382, 0.003605361887406413,
 420    -0.00836518656029184, 0.004615715410745561, 0.005105376617651312,
 421    -0.00944938657024391, 0.05602449420950007, 0.02722719610561933,
 422    0.01959357494748446, -0.0258655103753962, 0.1440733975689835,
 423    0.01446782819722976, 0.003718896062070054, 0.05825843045655135,
 424    -0.06230154142733073, -0.07833704962300169, 0.003160836143568724,
 425    -0.001169873777936648, 0.03471745590503304},
 426   {-0.03204352258752698, 0.01019272923862322, 0.04509668708733181,
 427    0.05756522429120813, -0.0004601149081726732, -0.0984718150777423,
 428    -0.01107826100664925, -0.005680277810520585, 0.01962359392320817,
 429    0.01550006899131986, 0.05143956925922197, 0.02462476682588468,
 430    -0.0888843861002653, -0.00171553583659411, 0.01606331750661664,
 431    0.001176847743518958, -0.02070972978912828, -0.000341523293579971,
 432    -0.002654732745607882, 0.02075709428885848},
 433   {0.03595199666430258, -0.02800219615234468, -0.04341570015493925,
 434    -0.0748275906176658, 0.0001051403676377422, 0.1137431321746627,
 435    0.005852087565974318, 0.003443037513847801, -0.02481931657706633,
 436    -0.003651181839831423, 0.03195794176786321, 0.04135411406392523,
 437    -0.07562030263210619, 0.001769332364699, -0.01984381173403915,
 438    -0.005029750745010152, 0.02649253902476472, 0.000518085571702734,
 439    0.001062936684474851, 0.01295950668914449},
 440   {-0.16164552322896, -0.0006050035060464324, 0.0258380054414968,
 441    0.003188424740960557, -0.0002058911341821877, 0.03157555987384681,
 442    -0.01678913462596107, 0.03096216145389774, -0.0133791110666919,
 443    0.1125249625204277, -0.00769017706442472, -0.02653938062180483,
 444    -0.002555329863523985, -0.00861833362947954, 0.01775148884754278,
 445    0.02529310679774722, 0.0826243417011238, -0.0001036728183032624,
 446    0.001963562313294209, -0.0935900561309786},
 447   {0.1652394174588469, -0.002814245280784351, -0.0328982001821263,
 448    -0.02000104712964131, 0.0002208121995725443, -0.02733462178511839,
 449    0.02648078162927627, -0.01788316626401427, 0.01630747623755998,
 450    0.1053849023838147, -0.005447706553811218, 0.01810876922536839,
 451    -0.001808914710282444, -0.007687912115607397, -0.01332593672114388,
 452    -0.02110750894891371, -0.07456116592983384, 0.000219072589592394,
 453    0.001270886972191055, -0.1083616930749109},
 454   {0.02453279389716254, -0.005820072356487439, 0.100260287284095,
 455    0.01277522280305745, -0.003184943445296999, 0.05814689527984152,
 456    -0.0934012278200201, -0.03017986487349484, -0.03136625380994165,
 457    0.00988668352785117, -0.00358900410973142, -0.02017443675004764,
 458    0.000915384582922184, -0.001460963415183106, -0.01370112443251124,
 459    0.1130040979284457, -0.1196161771323699, -0.0005800211204222045,
 460    -0.0006153403201024954, 0.00416806428223025},
 461   {-0.0778089244252535, -0.007055161182430869, -0.0349307504860869,
 462    -0.0811915584276571, -0.004689825871599125, -0.03726108871471753,
 463    0.1072225647141469, -0.00917015113070944, 0.01381628985996913,
 464    -0.00123227881492089, 0.001815954515275675, 0.005708744099349901,
 465    -0.0001448985044877925, -0.001306578795561384, -0.006992743514185243,
 466    0.1744720240732789, -0.05353628497814023, -0.0007613684227234787,
 467    -0.0003550282315997644, 0.01340106423804634},
 468   {-0.0159527329868513, -0.007622151568160798, -0.1389875105184963,
 469    0.1165051999914764, -0.002217810389087748, 0.01550003226513692,
 470    -0.07427664222230566, -0.003371438498619264, 0.01385754771325365,
 471    0.004759020167383304, 0.001624078805220564, 0.02011638303109029,
 472    -0.001717827082842178, -0.0007424036708598594, -0.003978884451898934,
 473    0.0866418927301209, -0.01280817739158123, -0.00023039242454603,
 474    0.002309205802479111, 0.0005926106991001195}};
 475
 476
 477 void init_protmats()
 478 {
 479   long l, m;
 480
 481   eigmat = (double *) Malloc (20 * sizeof(double));
 482   for (l = 0; l <= 19; l++)
 483     if (usejtt)
 484       eigmat[l] = jtteigmat[l];
 485     else {
 486       if (usepmb)
 487         eigmat[l] = pmbeigmat[l];
 488       else
 489         eigmat[l] = pameigmat[l];
 490       }
 491   probmat = (double **) Malloc (20 * sizeof(double *));
 492   for (l = 0; l <= 19; l++)
 493     for (m= 0; m <= 19; m++)
 494       if (usejtt)
 495         probmat[l] = jttprobmat[l];
 496       else {
 497         if (usepmb)
 498           probmat[l] = pmbprobmat[l];
 499         else
 500           probmat[l] = pamprobmat[l];
 501         }
 502 }  /* init_protmats */
 503
 504
 505 void getoptions()
 506 {
 507   /* interactively set options */
 508   long i, loopcount, loopcount2;
 509   Char ch;
 510   boolean didchangecat, didchangercat;
 511   double probsum;
 512
 513   fprintf(outfile, "\nAmino acid sequence Maximum Likelihood");
 514   fprintf(outfile, " method, version %s\n\n",VERSION);
 515   putchar('\n');
 516   ctgry = false;
 517   didchangecat = false;
 518   rctgry = false;
 519   didchangercat = false;
 520   categs = 1;
 521   rcategs = 1;
 522   auto_ = false;
 523   gama = false;
 524   global = false;
 525   hypstate = false;
 526   improve = false;
 527   invar = false;
 528   jumble = false;
 529   njumble = 1;
 530   lngths = false;
 531   lambda = 1.0;
 532   outgrno = 1;
 533   outgropt = false;
 534   trout = true;
 535   usertree = false;
 536   weights = false;
 537   printdata = false;
 538   progress = true;
 539   treeprint = true;
 540   usejtt = true;
 541   usepmb = false;
 542   usepam = false;
 543   interleaved = true;
 544   loopcount = 0;
 545   for (;;){
 546     cleerhome();
 547     printf("Amino acid sequence Maximum Likelihood");
 548     printf(" method, version %s\n\n",VERSION);
 549     printf("Settings for this run:\n");
 550     printf("  U                 Search for best tree?  %s\n",
 551            (usertree ? "No, use user trees in input file" : "Yes"));
 552     if (usertree) {
 553       printf("  L          Use lengths from user trees?  %s\n",
 554                (lngths ? "Yes" : "No"));
 555     }
 556     printf("  P    JTT, PMB or PAM probability model?  %s\n",
 557       usejtt ? "Jones-Taylor-Thornton" :
 558       usepmb ? "Henikoff/Tillier PMB" : "Dayhoff PAM");
 559     printf("  C                One category of sites?");
 560     if (!ctgry || categs == 1)
 561       printf("  Yes\n");
 562     else
 563       printf("  %ld categories of sites\n", categs);
 564     printf("  R           Rate variation among sites?");
 565     if (!rctgry)
 566       printf("  constant rate of change\n");
 567     else {
 568       if (gama)
 569         printf("  Gamma distributed rates\n");
 570       else {
 571         if (invar)
 572           printf("  Gamma+Invariant sites\n");
 573         else
 574           printf("  user-defined HMM of rates\n");
 575       }
 576       printf("  A   Rates at adjacent sites correlated?");
 577       if (!auto_)
 578         printf("  No, they are independent\n");
 579       else
 580         printf("  Yes, mean block length =%6.1f\n", 1.0 / lambda);
 581     }
 582     printf("  W                       Sites weighted?  %s\n",
 583            (weights ? "Yes" : "No"));
 584     if (!usertree) {
 585       printf("  S        Speedier but rougher analysis?  %s\n",
 586              (improve ? "No, not rough" : "Yes"));
 587       printf("  G                Global rearrangements?  %s\n",
 588              (global ? "Yes" : "No"));
 589     }
 590     if (!usertree) {
 591       printf("  J   Randomize input order of sequences?");
 592       if (jumble)
 593         printf("  Yes (seed =%8ld,%3ld times)\n", inseed0, njumble);
 594       else
 595         printf("  No. Use input order\n");
 596     }
 597     printf("  O                        Outgroup root?  %s%3ld\n",
 598            (outgropt ? "Yes, at sequence number" :
 599                        "No, use as outgroup species"),outgrno);
 600     printf("  M           Analyze multiple data sets?");
 601     if (mulsets)
 602       printf("  Yes, %2ld %s\n", datasets,
 603                (justwts ? "sets of weights" : "data sets"));
 604     else
 605       printf("  No\n");
 606     printf("  I          Input sequences interleaved?  %s\n",
 607            (interleaved ? "Yes" : "No, sequential"));
 608     printf("  0   Terminal type (IBM PC, ANSI, none)?  %s\n",
 609            (ibmpc ? "IBM PC" : ansi  ? "ANSI" : "(none)"));
 610     printf("  1    Print out the data at start of run  %s\n",
 611            (printdata ? "Yes" : "No"));
 612     printf("  2  Print indications of progress of run  %s\n",
 613            (progress ? "Yes" : "No"));
 614     printf("  3                        Print out tree  %s\n",
 615            (treeprint ? "Yes" : "No"));
 616     printf("  4       Write out trees onto tree file?  %s\n",
 617            (trout ? "Yes" : "No"));
 618     printf("  5   Reconstruct hypothetical sequences?  %s\n",
 619            (hypstate ? "Yes" : "No"));
 620     printf("\n  Y to accept these or type the letter for one to change\n");
 621 #ifdef WIN32
 622     phyFillScreenColor();
 623 #endif
 624     scanf("%c%*[^\n]", &ch);
 625     getchar();
 626     if (ch == '\n')
 627       ch = ' ';
 628     uppercase(&ch);
 629     if (ch == 'Y')
 630       break;
 631     if (strchr("UPLCRAWSGJOMI012345",ch) != NULL){
 632       switch (ch) {
 633
 634       case 'C':
 635         ctgry = !ctgry;
 636         if (ctgry) {
 637           printf("\nSitewise user-assigned categories:\n\n");
 638           initcatn(&categs);
 639           if (rate){
 640             free(rate);
 641           }
 642           rate    = (double *) Malloc(categs * sizeof(double));
 643           didchangecat = true;
 644           initcategs(categs, rate);
 645         }
 646         break;
 647
 648       case 'P':
 649         if (usejtt) {
 650           usejtt = false;
 651           usepmb = true;
 652         } else {
 653             if (usepmb) {
 654               usepmb = false;
 655               usepam = true;
 656             } else {
 657               usepam = false;
 658               usejtt = true;
 659             }
 660           }
 661         break;
 662
 663       case 'R':
 664         if (!rctgry) {
 665           rctgry = true;
 666           gama = true;
 667         } else {
 668           if (gama) {
 669             gama = false;
 670             invar = true;
 671           } else {
 672             if (invar)
 673               invar = false;
 674             else
 675               rctgry = false;
 676           }
 677         }
 678         break;
 679
 680       case 'A':
 681         auto_ = !auto_;
 682         if (auto_)
 683           initlambda(&lambda);
 684         break;
 685
 686       case 'W':
 687         weights = !weights;
 688         break;
 689
 690       case 'S':
 691         improve = !improve;
 692         break;
 693
 694       case 'G':
 695         global = !global;
 696         break;
 697
 698       case 'J':
 699         jumble = !jumble;
 700         if (jumble)
 701           initjumble(&inseed, &inseed0, seed, &njumble);
 702         else njumble = 1;
 703         break;
 704
 705       case 'L':
 706         lngths = !lngths;
 707         break;
 708
 709       case 'O':
 710         outgropt = !outgropt;
 711         if (outgropt)
 712           initoutgroup(&outgrno, spp);
 713         break;
 714
 715       case 'U':
 716         usertree = !usertree;
 717         break;
 718
 719       case 'M':
 720         mulsets = !mulsets;
 721         if (mulsets) {
 722           printf("Multiple data sets or multiple weights?");
 723           loopcount2 = 0;
 724           do {
 725             printf(" (type D or W)\n");
 726 #ifdef WIN32
 727             phyFillScreenColor();
 728 #endif
 729             scanf("%c%*[^\n]", &ch2);
 730             getchar();
 731             if (ch2 == '\n')
 732                 ch2 = ' ';
 733             uppercase(&ch2);
 734             countup(&loopcount2, 10);
 735           } while ((ch2 != 'W') && (ch2 != 'D'));
 736           justwts = (ch2 == 'W');
 737           if (justwts)
 738             justweights(&datasets);
 739           else
 740             initdatasets(&datasets);
 741           if (!jumble) {
 742             jumble = true;
 743             initjumble(&inseed, &inseed0, seed, &njumble);
 744           }
 745         }
 746         break;
 747
 748       case 'I':
 749         interleaved = !interleaved;
 750         break;
 751
 752       case '0':
 753         initterminal(&ibmpc, &ansi);
 754         break;
 755
 756       case '1':
 757         printdata = !printdata;
 758         break;
 759
 760       case '2':
 761         progress = !progress;
 762         break;
 763
 764       case '3':
 765         treeprint = !treeprint;
 766         break;
 767
 768       case '4':
 769         trout = !trout;
 770         break;
 771
 772       case '5':
 773         hypstate = !hypstate;
 774         break;
 775       }
 776     } else
 777       printf("Not a possible option!\n");
 778     countup(&loopcount, 100);
 779   }
 780   if (gama || invar) {
 781     loopcount = 0;
 782     do {
 783       printf(
 784 "\nCoefficient of variation of substitution rate among sites (must be positive)\n");
 785       printf(
 786   " In gamma distribution parameters, this is 1/(square root of alpha)\n");
 787 #ifdef WIN32
 788       phyFillScreenColor();
 789 #endif
 790       scanf("%lf%*[^\n]", &cv);
 791       getchar();
 792       countup(&loopcount, 10);
 793     } while (cv <= 0.0);
 794     alpha = 1.0 / (cv * cv);
 795   }
 796   if (!rctgry)
 797     auto_ = false;
 798   if (rctgry) {
 799     printf("\nRates in HMM");
 800     if (invar)
 801       printf(" (including one for invariant sites)");
 802     printf(":\n");
 803     initcatn(&rcategs);
 804     if (probcat){
 805       free(probcat);
 806       free(rrate);
 807     }
 808     probcat = (double *) Malloc(rcategs * sizeof(double));
 809     rrate   = (double *) Malloc(rcategs * sizeof(double));
 810     didchangercat = true;
 811     if (gama)
 812       initgammacat(rcategs, alpha, rrate, probcat);
 813     else {
 814       if (invar) {
 815         loopcount = 0;
 816         do {
 817           printf("Fraction of invariant sites?\n");
 818           scanf("%lf%*[^\n]", &invarfrac);
 819           getchar();
 820           countup (&loopcount, 10);
 821         } while ((invarfrac <= 0.0) || (invarfrac >= 1.0));
 822         initgammacat(rcategs-1, alpha, rrate, probcat);
 823         for (i = 0; i < rcategs-1; i++)
 824           probcat[i] = probcat[i]*(1.0-invarfrac);
 825         probcat[rcategs-1] = invarfrac;
 826         rrate[rcategs-1] = 0.0;
 827       } else {
 828         initcategs(rcategs, rrate);
 829         initprobcat(rcategs, &probsum, probcat);
 830       }
 831     }
 832   }
 833   if (!didchangercat){
 834     rrate      = (double *) Malloc(rcategs*sizeof(double));
 835     probcat    = (double *) Malloc(rcategs*sizeof(double));
 836     rrate[0]   = 1.0;
 837     probcat[0] = 1.0;
 838   }
 839   if (!didchangecat) {
 840     rate       = (double *) Malloc(categs*sizeof(double));
 841     rate[0]    = 1.0;
 842   }
 843   init_protmats();
 844 }  /* getoptions */
 845
 846
 847 void makeprotfreqs()
 848 {
 849   /* calculate amino acid frequencies based on eigmat */
 850   long i, mineig;
 851
 852   mineig = 0;
 853   for (i = 0; i <= 19; i++)
 854     if (fabs(eigmat[i]) < fabs(eigmat[mineig]))
 855       mineig = i;
 856   memcpy(freqaa, probmat[mineig], 20 * sizeof(double));
 857   for (i = 0; i <= 19; i++)
 858     freqaa[i] = fabs(freqaa[i]);
 859 } /* makeprotfreqs */
 860
 861 void reallocsites()
 862 {
 863   long i;
 864   for (i = 0; i < spp; i++)
 865     y[i] = (Char *) Malloc(sites*sizeof(Char));
 866
 867   free(category);
 868   free(weight);
 869   free(alias);
 870   free(ally);
 871   free(location);
 872   free(aliasweight);
 873
 874   category    = (long *) Malloc(sites*sizeof(long));
 875   weight      = (long *) Malloc(sites*sizeof(long));
 876   alias       = (long *) Malloc(sites*sizeof(long));
 877   ally        = (long *) Malloc(sites*sizeof(long));
 878   location    = (long *) Malloc(sites*sizeof(long));
 879   aliasweight = (long *) Malloc(sites*sizeof(long));
 880   for (i = 0; i < sites; i++)
 881       category[i] = 1;
 882   for (i = 0; i < sites; i++)
 883     weight[i] = 1;
 884   makeweights();
 885 }
 886
 887 void allocrest()
 888 {
 889   long i;
 890
 891   y = (Char **) Malloc(spp*sizeof(Char *));
 892   for (i = 0; i < spp; i++)
 893     y[i] = (Char *) Malloc(sites*sizeof(Char));
 894   nayme       = (naym *) Malloc(spp*sizeof(naym));
 895   enterorder  = (long *) Malloc(spp*sizeof(long));
 896   category    = (long *) Malloc(sites*sizeof(long));
 897   weight      = (long *) Malloc(sites*sizeof(long));
 898   alias       = (long *) Malloc(sites*sizeof(long));
 899   ally        = (long *) Malloc(sites*sizeof(long));
 900   location    = (long *) Malloc(sites*sizeof(long));
 901   aliasweight = (long *) Malloc(sites*sizeof(long));
 902 }  /* allocrest */
 903
 904
 905 void doinit()
 906 { /* initializes variables */
 907   inputnumbers(&spp, &sites, &nonodes2, 1);
 908   getoptions();
 909   if (!usertree)
 910     nonodes2--;
 911   makeprotfreqs();
 912   if (printdata)
 913     fprintf(outfile, "%2ld species, %3ld  sites\n", spp, sites);
 914   alloctree(&curtree.nodep, nonodes2, usertree);
 915   allocrest();
 916   if (usertree)
 917     return;
 918   alloctree(&bestree.nodep, nonodes2, 0);
 919   alloctree(&priortree.nodep, nonodes2, 0);
 920   if (njumble <= 1)
 921     return;
 922   alloctree(&bestree2.nodep, nonodes2, 0);
 923 }  /* doinit */
 924
 925
 926 void inputoptions()
 927 {
 928   long i;
 929
 930   if (!firstset) {
 931     samenumsp(&sites, ith);
 932     reallocsites();
 933   }
 934   if (firstset) {
 935     for (i = 0; i < sites; i++)
 936       category[i] = 1;
 937     for (i = 0; i < sites; i++)
 938       weight[i] = 1;
 939   }
 940   if (justwts || weights)
 941     inputweights(sites, weight, &weights);
 942   weightsum = 0;
 943   for (i = 0; i < sites; i++)
 944     weightsum += weight[i];
 945   if ((ctgry && categs > 1) && (firstset || !justwts)) {
 946     inputcategs(0, sites, category, categs, "ProML");
 947     if (printdata)
 948       printcategs(outfile, sites, category, "Site categories");
 949   }
 950   if (weights && printdata)
 951     printweights(outfile, 0, sites, weight, "Sites");
 952   fprintf(outfile, "%s model of amino acid change\n\n",
 953           (usejtt ? "Jones-Taylor-Thornton" :
 954            usepmb ? "Henikoff/Tillier PMB" : "Dayhoff PAM"));
 955 }  /* inputoptions */
 956
 957
 958 void input_protdata(long chars)
 959 {
 960   /* input the names and sequences for each species */
 961   /* used by proml */
 962   long i, j, k, l, basesread, basesnew;
 963   Char charstate;
 964   boolean allread, done;
 965
 966   if (printdata)
 967     headings(chars, "Sequences", "---------");
 968   basesread = 0;
 969   basesnew = 0;
 970   allread = false;
 971   while (!(allread)) {
 972     /* eat white space -- if the separator line has spaces on it*/
 973     do {
 974       charstate = gettc(infile);
 975     } while (charstate == ' ' || charstate == '\t');
 976     ungetc(charstate, infile);
 977     if (eoln(infile))
 978       scan_eoln(infile);
 979     i = 1;
 980     while (i <= spp) {
 981       if ((interleaved && basesread == 0) || !interleaved)
 982         initname(i - 1);
 983       j = (interleaved) ? basesread : 0;
 984       done = false;
 985       while (!done && !eoff(infile)) {
 986         if (interleaved)
 987           done = true;
 988         while (j < chars && !(eoln(infile) || eoff(infile))) {
 989           charstate = gettc(infile);
 990           if (charstate == '\n' || charstate == '\t')
 991             charstate = ' ';
 992           if (charstate == ' ' || (charstate >= '0' && charstate <= '9'))
 993             continue;
 994           uppercase(&charstate);
 995           if ((strchr("ABCDEFGHIKLMNPQRSTVWXYZ*?-", charstate)) == NULL) {
 996         printf("ERROR: bad amino acid: %c at position %ld of species %ld\n",
 997                    charstate, j+1, i);
 998             if (charstate == '.') {
 999           printf("       Periods (.) may not be used as gap characters.\n");
1000           printf("       The correct gap character is (-)\n");
1001             }
1002             exxit(-1);
1003           }
1004           j++;
1005           y[i - 1][j - 1] = charstate;
1006         }
1007         if (interleaved)
1008           continue;
1009         if (j < chars)
1010           scan_eoln(infile);
1011         else if (j == chars)
1012           done = true;
1013       }
1014       if (interleaved && i == 1)
1015         basesnew = j;
1016
1017       scan_eoln(infile);
1018
1019       if ((interleaved && j != basesnew) ||
1020           (!interleaved && j != chars)) {
1021         printf("ERROR: SEQUENCES OUT OF ALIGNMENT AT POSITION %ld.\n", j);
1022         exxit(-1);
1023       }
1024       i++;
1025     }
1026
1027     if (interleaved) {
1028       basesread = basesnew;
1029       allread = (basesread == chars);
1030     } else
1031       allread = (i > spp);
1032   }
1033   if (!printdata)
1034     return;
1035   for (i = 1; i <= ((chars - 1) / 60 + 1); i++) {
1036     for (j = 1; j <= spp; j++) {
1037       for (k = 0; k < nmlngth; k++)
1038         putc(nayme[j - 1][k], outfile);
1039       fprintf(outfile, "   ");
1040       l = i * 60;
1041       if (l > chars)
1042         l = chars;
1043       for (k = (i - 1) * 60 + 1; k <= l; k++) {
1044         if (j > 1 && y[j - 1][k - 1] == y[0][k - 1])
1045           charstate = '.';
1046         else
1047           charstate = y[j - 1][k - 1];
1048         putc(charstate, outfile);
1049         if (k % 10 == 0 && k % 60 != 0)
1050           putc(' ', outfile);
1051       }
1052       putc('\n', outfile);
1053     }
1054     putc('\n', outfile);
1055   }
1056   putc('\n', outfile);
1057 }  /* input_protdata */
1058
1059
1060 void makeweights()
1061 {
1062   /* make up weights vector to avoid duplicate computations */
1063   long i;
1064
1065   for (i = 1; i <= sites; i++) {
1066     alias[i - 1] = i;
1067     ally[i - 1] = 0;
1068     aliasweight[i - 1] = weight[i - 1];
1069     location[i - 1] = 0;
1070   }
1071   sitesort2   (sites, aliasweight);
1072   sitecombine2(sites, aliasweight);
1073   sitescrunch2(sites, 1, 2, aliasweight);
1074   for (i = 1; i <= sites; i++) {
1075     if (aliasweight[i - 1] > 0)
1076       endsite = i;
1077   }
1078   for (i = 1; i <= endsite; i++) {
1079     location[alias[i - 1] - 1] = i;
1080     ally[alias[i - 1] - 1] = alias[i - 1];
1081   }
1082   term = (double **) Malloc(endsite * sizeof(double *));
1083   for (i = 0; i < endsite; i++)
1084      term[i] = (double *) Malloc(rcategs * sizeof(double));
1085   slopeterm = (double **) Malloc(endsite * sizeof(double *));
1086   for (i = 0; i < endsite; i++)
1087      slopeterm[i] = (double *) Malloc(rcategs * sizeof(double));
1088   curveterm = (double **) Malloc(endsite * sizeof(double *));
1089   for (i = 0; i < endsite; i++)
1090      curveterm[i] = (double *) Malloc(rcategs * sizeof(double));
1091   mp = (vall *) Malloc(sites*sizeof(vall));
1092   contribution = (contribarr *) Malloc(endsite*sizeof(contribarr));
1093 }  /* makeweights */
1094
1095
1096 void prot_makevalues(long categs, pointarray treenode, long endsite,
1097                         long spp, sequence y, steptr alias)
1098 {
1099   /* set up fractional likelihoods at tips   */
1100   /* a version of makevalues2 found in seq.c */
1101   /* used by proml                             */
1102   long i, j, k, l;
1103   long b;
1104
1105   for (k = 0; k < endsite; k++) {
1106     j = alias[k];
1107     for (i = 0; i < spp; i++) {
1108       for (l = 0; l < categs; l++) {
1109         memset(treenode[i]->protx[k][l], 0, sizeof(double)*20);
1110         switch (y[i][j - 1]) {
1111
1112         case 'A':
1113           treenode[i]->protx[k][l][0] = 1.0;
1114           break;
1115
1116         case 'R':
1117           treenode[i]->protx[k][l][(long)arginine   - (long)alanine] = 1.0;
1118           break;
1119
1120         case 'N':
1121           treenode[i]->protx[k][l][(long)asparagine - (long)alanine] = 1.0;
1122           break;
1123
1124         case 'D':
1125           treenode[i]->protx[k][l][(long)aspartic   - (long)alanine] = 1.0;
1126           break;
1127
1128         case 'C':
1129           treenode[i]->protx[k][l][(long)cysteine   - (long)alanine] = 1.0;
1130           break;
1131
1132         case 'Q':
1133           treenode[i]->protx[k][l][(long)glutamine  - (long)alanine] = 1.0;
1134           break;
1135
1136         case 'E':
1137           treenode[i]->protx[k][l][(long)glutamic   - (long)alanine] = 1.0;
1138           break;
1139
1140         case 'G':
1141           treenode[i]->protx[k][l][(long)glycine    - (long)alanine] = 1.0;
1142           break;
1143
1144         case 'H':
1145           treenode[i]->protx[k][l][(long)histidine  - (long)alanine] = 1.0;
1146           break;
1147
1148         case 'I':
1149           treenode[i]->protx[k][l][(long)isoleucine - (long)alanine] = 1.0;
1150           break;
1151
1152         case 'L':
1153           treenode[i]->protx[k][l][(long)leucine    - (long)alanine] = 1.0;
1154           break;
1155
1156         case 'K':
1157           treenode[i]->protx[k][l][(long)lysine     - (long)alanine] = 1.0;
1158           break;
1159
1160         case 'M':
1161           treenode[i]->protx[k][l][(long)methionine - (long)alanine] = 1.0;
1162           break;
1163
1164         case 'F':
1165           treenode[i]->protx[k][l][(long)phenylalanine - (long)alanine] = 1.0;
1166           break;
1167
1168         case 'P':
1169           treenode[i]->protx[k][l][(long)proline    - (long)alanine] = 1.0;
1170           break;
1171
1172         case 'S':
1173           treenode[i]->protx[k][l][(long)serine     - (long)alanine] = 1.0;
1174           break;
1175
1176         case 'T':
1177           treenode[i]->protx[k][l][(long)threonine  - (long)alanine] = 1.0;
1178           break;
1179
1180         case 'W':
1181           treenode[i]->protx[k][l][(long)tryptophan - (long)alanine] = 1.0;
1182           break;
1183
1184         case 'Y':
1185           treenode[i]->protx[k][l][(long)tyrosine   - (long)alanine] = 1.0;
1186           break;
1187
1188         case 'V':
1189           treenode[i]->protx[k][l][(long)valine     - (long)alanine] = 1.0;
1190           break;
1191
1192         case 'B':
1193           treenode[i]->protx[k][l][(long)asparagine - (long)alanine] = 1.0;
1194           treenode[i]->protx[k][l][(long)aspartic   - (long)alanine] = 1.0;
1195           break;
1196
1197         case 'Z':
1198           treenode[i]->protx[k][l][(long)glutamine  - (long)alanine] = 1.0;
1199           treenode[i]->protx[k][l][(long)glutamic   - (long)alanine] = 1.0;
1200           break;
1201
1202         case 'X':                /* unknown aa                            */
1203           for (b = 0; b <= 19; b++)
1204             treenode[i]->protx[k][l][b] = 1.0;
1205           break;
1206
1207         case '?':                /* unknown aa                            */
1208           for (b = 0; b <= 19; b++)
1209             treenode[i]->protx[k][l][b] = 1.0;
1210           break;
1211
1212         case '*':                /* stop codon symbol                    */
1213           for (b = 0; b <= 19; b++)
1214             treenode[i]->protx[k][l][b] = 1.0;
1215           break;
1216
1217         case '-':                /* deletion event-absent data or aa */
1218           for (b = 0; b <= 19; b++)
1219             treenode[i]->protx[k][l][b] = 1.0;
1220           break;
1221         }
1222       }
1223     }
1224   }
1225 }  /* prot_makevalues */
1226
1227
1228 void free_pmatrix(long sib)
1229 {
1230   long j,k,l;
1231
1232   for (j = 0; j < rcategs; j++) {
1233     for (k = 0; k < categs; k++) {
1234       for (l = 0; l < 20; l++)
1235         free(pmatrices[sib][j][k][l]);
1236       free(pmatrices[sib][j][k]);
1237     }
1238     free(pmatrices[sib][j]);
1239   }
1240   free(pmatrices[sib]);
1241 }
1242
1243 void alloc_pmatrix(long sib)
1244 {
1245   /* Allocate memory for a new pmatrix.  Called iff num_sibs>max_num_sibs */
1246   long j, k, l;
1247   double ****temp_matrix;
1248
1249   temp_matrix = (double ****) Malloc (rcategs * sizeof(double ***));
1250   for (j = 0; j < rcategs; j++) {
1251     temp_matrix[j] = (double ***) Malloc(categs * sizeof(double **));
1252     for (k = 0; k < categs; k++) {
1253       temp_matrix[j][k] = (double **) Malloc(20 * sizeof (double *));
1254       for (l = 0; l < 20; l++)
1255         temp_matrix[j][k][l] = (double *) Malloc(20 * sizeof(double));
1256     }
1257   }
1258   pmatrices[sib] = temp_matrix;
1259   max_num_sibs++;
1260 }  /* alloc_pmatrix */
1261
1262 void prot_freetable()
1263 {
1264   long i,j,k,l;
1265   for (j = 0; j < rcategs; j++) {
1266     for (k = 0; k < categs; k++) {
1267       for (l = 0; l < 20; l++)
1268         free(ddpmatrix[j][k][l]);
1269       free(ddpmatrix[j][k]);
1270     }
1271     free(ddpmatrix[j]);
1272   }
1273   free(ddpmatrix);
1274
1275   for (j = 0; j < rcategs; j++) {
1276     for (k = 0; k < categs; k++) {
1277       for (l = 0; l < 20; l++)
1278         free(dpmatrix[j][k][l]);
1279       free(dpmatrix[j][k]);
1280     }
1281     free(dpmatrix[j]);
1282   }
1283   free(dpmatrix);
1284
1285
1286   for (j = 0; j < rcategs; j++)
1287     free(tbl[j]);
1288   free(tbl);
1289
1290   for ( i = 0 ; i < max_num_sibs ; i++ )
1291     free_pmatrix(i);
1292   free(pmatrices);
1293 }
1294
1295 void prot_inittable()
1296 {
1297   /* Define a lookup table. Precompute values and print them out in tables */
1298   /* Allocate memory for the pmatrices, dpmatices and ddpmatrices          */
1299   long i, j, k, l;
1300   double sumrates;
1301
1302   /* Allocate memory for pmatrices, the array of pointers to pmatrices     */
1303
1304   pmatrices = (double *****) Malloc ( spp * sizeof(double ****));
1305
1306   /* Allocate memory for the first 2 pmatrices, the matrix of conversion   */
1307   /* probabilities, but only once per run (aka not on the second jumble.   */
1308
1309     alloc_pmatrix(0);
1310     alloc_pmatrix(1);
1311
1312   /*  Allocate memory for one dpmatrix, the first derivative matrix        */
1313
1314   dpmatrix = (double ****) Malloc( rcategs * sizeof(double ***));
1315   for (j = 0; j < rcategs; j++) {
1316     dpmatrix[j] = (double ***) Malloc( categs * sizeof(double **));
1317     for (k = 0; k < categs; k++) {
1318       dpmatrix[j][k] = (double **) Malloc( 20 * sizeof(double *));
1319       for (l = 0; l < 20; l++)
1320         dpmatrix[j][k][l] = (double *) Malloc( 20 * sizeof(double));
1321     }
1322   }
1323
1324   /*  Allocate memory for one ddpmatrix, the second derivative matrix      */
1325   ddpmatrix = (double ****) Malloc( rcategs * sizeof(double ***));
1326   for (j = 0; j < rcategs; j++) {
1327     ddpmatrix[j] = (double ***) Malloc( categs * sizeof(double **));
1328     for (k = 0; k < categs; k++) {
1329       ddpmatrix[j][k] = (double **) Malloc( 20 * sizeof(double *));
1330       for (l = 0; l < 20; l++)
1331         ddpmatrix[j][k][l] = (double *) Malloc( 20 * sizeof(double));
1332     }
1333   }
1334
1335   /* Allocate memory and assign values to tbl, the matrix of possible rates*/
1336
1337   tbl = (double **) Malloc( rcategs * sizeof(double *));
1338   for (j = 0; j < rcategs; j++)
1339     tbl[j] = (double *) Malloc( categs * sizeof(double));
1340
1341   for (j = 0; j < rcategs; j++)
1342     for (k = 0; k < categs; k++)
1343       tbl[j][k] = rrate[j]*rate[k];
1344
1345   sumrates = 0.0;
1346   for (i = 0; i < endsite; i++) {
1347     for (j = 0; j < rcategs; j++)
1348       sumrates += aliasweight[i] * probcat[j]
1349         * tbl[j][category[alias[i] - 1] - 1];
1350   }
1351   sumrates /= (double)sites;
1352   for (j = 0; j < rcategs; j++)
1353     for (k = 0; k < categs; k++) {
1354       tbl[j][k] /= sumrates;
1355     }
1356
1357   if(jumb > 1)
1358     return;
1359
1360   if (gama) {
1361     fprintf(outfile, "\nDiscrete approximation to gamma distributed rates\n");
1362     fprintf(outfile,
1363     " Coefficient of variation of rates = %f  (alpha = %f)\n",
1364       cv, alpha);
1365   }
1366   if (rcategs > 1) {
1367     fprintf(outfile, "\nStates in HMM   Rate of change    Probability\n\n");
1368     for (i = 0; i < rcategs; i++)
1369       if (probcat[i] < 0.0001)
1370         fprintf(outfile, "%9ld%16.3f%20.6f\n", i+1, rrate[i], probcat[i]);
1371       else if (probcat[i] < 0.001)
1372           fprintf(outfile, "%9ld%16.3f%19.5f\n", i+1, rrate[i], probcat[i]);
1373         else if (probcat[i] < 0.01)
1374             fprintf(outfile, "%9ld%16.3f%18.4f\n", i+1, rrate[i], probcat[i]);
1375           else
1376             fprintf(outfile, "%9ld%16.3f%17.3f\n", i+1, rrate[i], probcat[i]);
1377     putc('\n', outfile);
1378     if (auto_)
1379       fprintf(outfile,
1380      "Expected length of a patch of sites having the same rate = %8.3f\n",
1381              1/lambda);
1382     putc('\n', outfile);
1383   }
1384   if (categs > 1) {
1385     fprintf(outfile, "\nSite category   Rate of change\n\n");
1386     for (k = 0; k < categs; k++)
1387       fprintf(outfile, "%9ld%16.3f\n", k+1, rate[k]);
1388   }
1389   if ((rcategs  > 1) || (categs >> 1))
1390     fprintf(outfile, "\n\n");
1391 }  /* prot_inittable */
1392
1393
1394 void getinput()
1395 {
1396   /* reads the input data */
1397   if (!justwts || firstset)
1398     inputoptions();
1399   if (!justwts || firstset)
1400     input_protdata(sites);
1401   if ( !firstset ) freelrsaves();
1402   makeweights();
1403   alloclrsaves();
1404   setuptree2(curtree);
1405   if (!usertree) {
1406     setuptree2(bestree);
1407     setuptree2(priortree);
1408     if (njumble > 1)
1409       setuptree2(bestree2);
1410   }
1411   prot_allocx(nonodes2, rcategs, curtree.nodep, usertree);
1412   if (!usertree) {
1413     prot_allocx(nonodes2, rcategs, bestree.nodep, 0);
1414     prot_allocx(nonodes2, rcategs, priortree.nodep, 0);
1415     if (njumble > 1)
1416       prot_allocx(nonodes2, rcategs, bestree2.nodep, 0);
1417   }
1418   prot_makevalues(rcategs, curtree.nodep, endsite, spp, y, alias);
1419 }  /* getinput */
1420
1421
1422 void inittravtree(node *p)
1423 {
1424   /* traverse tree to set initialized and v to initial values */
1425   node* q;
1426
1427   p->initialized = false;
1428   p->back->initialized = false;
1429   if ((!lngths) || p->iter) {
1430     p->v = initialv;
1431     p->back->v = initialv;
1432   }
1433
1434   if ( !p->tip ) {
1435     q = p->next;
1436     while ( q != p ) {
1437       inittravtree(q->back);
1438       q = q->next;
1439     }
1440   }
1441 }  /* inittravtree */
1442
1443
1444 void prot_nuview(node *p)
1445 {
1446   long i, j, k, l, m, num_sibs, sib_index;
1447   node *sib_ptr, *sib_back_ptr;
1448   psitelike prot_xx, x2;
1449   double lw, prod7;
1450   double **pmat;
1451   double maxx;
1452   double correction;
1453
1454   /* Figure out how many siblings the current node has  */
1455   /* and be sure that pmatrices is large enough         */
1456   num_sibs = count_sibs(p);
1457   for (i = 0; i < num_sibs; i++)
1458     if (pmatrices[i] == NULL)
1459       alloc_pmatrix(i);
1460
1461   /* Recursive calls, should be called for all children */
1462   sib_ptr = p;
1463   for (i=0 ; i < num_sibs; i++) {
1464     sib_ptr      = sib_ptr->next;
1465     sib_back_ptr = sib_ptr->back;
1466
1467     if (!sib_back_ptr->tip &&
1468         !sib_back_ptr->initialized)
1469       prot_nuview(sib_back_ptr);
1470   }
1471
1472   /* Make pmatrices for all possible combinations of category, rcateg      */
1473   /* and sib                                                               */
1474   sib_ptr = p;                                /* return to p */
1475   for (sib_index=0; sib_index < num_sibs; sib_index++) {
1476     sib_ptr      = sib_ptr->next;
1477     sib_back_ptr = sib_ptr->back;
1478
1479     lw = sib_back_ptr->v;
1480
1481     for (j = 0; j < rcategs; j++)
1482       for (k = 0; k < categs; k++)
1483         make_pmatrix(pmatrices[sib_index][j][k], NULL, NULL, 0, lw,
1484                                         tbl[j][k], eigmat, probmat);
1485   }
1486
1487   for (i = 0; i < endsite; i++) {
1488     maxx = 0;
1489     correction = 0;
1490
1491     k = category[alias[i]-1] - 1;
1492     for (j = 0; j < rcategs; j++) {
1493
1494       /* initialize to 1 all values of prot_xx */
1495       for (m = 0; m <= 19; m++)
1496         prot_xx[m] = 1;
1497
1498       sib_ptr = p;                        /* return to p */
1499       /* loop through all sibs and calculate likelihoods for all possible*/
1500       /* amino acid combinations                                         */
1501       for (sib_index=0; sib_index < num_sibs; sib_index++) {
1502         sib_ptr      = sib_ptr->next;
1503         sib_back_ptr = sib_ptr->back;
1504
1505         if ( j == 0)
1506           correction += sib_back_ptr->underflows[i];
1507
1508         memcpy(x2, sib_back_ptr->protx[i][j], sizeof(psitelike));
1509         pmat = pmatrices[sib_index][j][k];
1510         for (m = 0; m <= 19; m++) {
1511           prod7 = 0;
1512           for (l = 0; l <= 19; l++)
1513             prod7 += (pmat[m][l] * x2[l]);
1514           prot_xx[m] *= prod7;
1515           if ( prot_xx[m] > maxx && sib_index == (num_sibs - 1))
1516             maxx = prot_xx[m];
1517         }
1518       }
1519       /* And the final point of this whole function: */
1520       memcpy(p->protx[i][j], prot_xx, sizeof(psitelike));
1521     }
1522     p->underflows[i] = 0;
1523     if ( maxx < MIN_DOUBLE )
1524       fix_protx(p,i,maxx,rcategs);
1525     p->underflows[i] += correction;
1526   }
1527
1528   p->initialized = true;
1529 }  /* prot_nuview */
1530
1531
1532 void prot_slopecurv(node *p,double y,double *like,double *slope,double *curve)
1533 {
1534   /* compute log likelihood, slope and curvature at node p */
1535   long i, j, k, l, m, lai;
1536   double sum, sumc, sumterm, lterm, sumcs, sumcc, sum2, slope2, curve2;
1537   double frexm = 0;                        /* frexm = freqaa[m]*x1[m] */
1538                                         /* frexml = frexm*x2[l]    */
1539   double prod4m, prod5m, prod6m;        /* elements of prod4-5 for */
1540                                         /* each m                   */
1541   double **pmat, **dpmat, **ddpmat;        /* local pointers to global*/
1542                                         /* matrices                   */
1543   double prod4, prod5, prod6;
1544   contribarr thelike, nulike, nuslope, nucurve,
1545     theslope, thecurve, clai, cslai, cclai;
1546   node *q;
1547   psitelike x1, x2;
1548
1549   q = p->back;
1550   sum = 0.0;
1551   for (j = 0; j < rcategs; j++) {
1552     for (k = 0; k < categs; k++) {
1553       make_pmatrix(pmatrices[0][j][k], dpmatrix[j][k], ddpmatrix[j][k],
1554                                         2, y, tbl[j][k], eigmat, probmat);
1555     }
1556   }
1557   for (i = 0; i < endsite; i++) {
1558     k = category[alias[i]-1] - 1;
1559     for (j = 0; j < rcategs; j++) {
1560       memcpy(x1, p->protx[i][j], sizeof(psitelike));
1561       memcpy(x2, q->protx[i][j], sizeof(psitelike));
1562       pmat = pmatrices[0][j][k];
1563       dpmat = dpmatrix[j][k];
1564       ddpmat = ddpmatrix[j][k];
1565       prod4 = 0.0;
1566       prod5 = 0.0;
1567       prod6 = 0.0;
1568       for (m = 0; m <= 19; m++) {
1569         prod4m = 0.0;
1570         prod5m = 0.0;
1571         prod6m = 0.0;
1572         frexm = x1[m] * freqaa[m];
1573         for (l = 0; l <= 19; l++) {
1574           prod4m += x2[l] * pmat[m][l];
1575           prod5m += x2[l] * dpmat[m][l];
1576           prod6m += x2[l] * ddpmat[m][l];
1577         }
1578         prod4 += frexm * prod4m;
1579         prod5 += frexm * prod5m;
1580         prod6 += frexm * prod6m;
1581       }
1582       term[i][j] = prod4;
1583       slopeterm[i][j] = prod5;
1584       curveterm[i][j] = prod6;
1585     }
1586     sumterm = 0.0;
1587     for (j = 0; j < rcategs; j++)
1588       sumterm += probcat[j] * term[i][j];
1589     if (sumterm <= 0.0)
1590         sumterm = 0.000000001;        /* ? shouldn't get here ?? */
1591     lterm = log(sumterm) + p->underflows[i] + q->underflows[i];
1592     for (j = 0; j < rcategs; j++) {
1593       term[i][j] = term[i][j] / sumterm;
1594       slopeterm[i][j] = slopeterm[i][j] / sumterm;
1595       curveterm[i][j] = curveterm[i][j] / sumterm;
1596     }
1597     sum += (aliasweight[i] * lterm);
1598   }
1599   for (i = 0; i < rcategs; i++) {
1600     thelike[i] = 1.0;
1601     theslope[i] = 0.0;
1602     thecurve[i] = 0.0;
1603   }
1604   for (i = 0; i < sites; i++) {
1605     sumc = 0.0;
1606     sumcs = 0.0;
1607     sumcc = 0.0;
1608     for (k = 0; k < rcategs; k++) {
1609       sumc += probcat[k] * thelike[k];
1610       sumcs += probcat[k] * theslope[k];
1611       sumcc += probcat[k] * thecurve[k];
1612     }
1613     sumc *= lambda;
1614     sumcs *= lambda;
1615     sumcc *= lambda;
1616     if ((ally[i] > 0) && (location[ally[i]-1] > 0)) {
1617       lai = location[ally[i] - 1];
1618       memcpy(clai, term[lai - 1], rcategs*sizeof(double));
1619       memcpy(cslai, slopeterm[lai - 1], rcategs*sizeof(double));
1620       memcpy(cclai, curveterm[lai - 1], rcategs*sizeof(double));
1621       if (weight[i] > 1) {
1622         for (j = 0; j < rcategs; j++) {
1623           if (clai[j] > 0.0)
1624             clai[j] = exp(weight[i]*log(clai[j]));
1625           else clai[j] = 0.0;
1626           if (cslai[j] > 0.0)
1627             cslai[j] = exp(weight[i]*log(cslai[j]));
1628           else cslai[j] = 0.0;
1629           if (cclai[j] > 0.0)
1630             cclai[j] = exp(weight[i]*log(cclai[j]));
1631           else cclai[j] = 0.0;
1632         }
1633       }
1634       for (j = 0; j < rcategs; j++) {
1635         nulike[j]  = ((1.0 - lambda) * thelike[j]  + sumc) *  clai[j];
1636         nuslope[j] = ((1.0 - lambda) * theslope[j] + sumcs) * clai[j]
1637                    + ((1.0 - lambda) * thelike[j]  + sumc) *  cslai[j];
1638         nucurve[j] = ((1.0 - lambda) * thecurve[j] + sumcc) * clai[j]
1639              + 2.0 * ((1.0 - lambda) * theslope[j] + sumcs) * cslai[j]
1640                    + ((1.0 - lambda) * thelike[j]  + sumc) *  cclai[j];
1641       }
1642     } else {
1643       for (j = 0; j < rcategs; j++) {
1644         nulike[j]  = ((1.0 - lambda) * thelike[j]  + sumc);
1645         nuslope[j] = ((1.0 - lambda) * theslope[j] + sumcs);
1646         nucurve[j] = ((1.0 - lambda) * thecurve[j] + sumcc);
1647       }
1648     }
1649     memcpy(thelike, nulike, rcategs*sizeof(double));
1650     memcpy(theslope, nuslope, rcategs*sizeof(double));
1651     memcpy(thecurve, nucurve, rcategs*sizeof(double));
1652   }
1653   sum2 = 0.0;
1654   slope2 = 0.0;
1655   curve2 = 0.0;
1656   for (i = 0; i < rcategs; i++) {
1657     sum2 += probcat[i] * thelike[i];
1658     slope2 += probcat[i] * theslope[i];
1659     curve2 += probcat[i] * thecurve[i];
1660   }
1661   sum += log(sum2);
1662   (*like) = sum;
1663   (*slope) = slope2 / sum2;
1664   (*curve) = (curve2 - slope2 * slope2 / sum2) / sum2;
1665 } /* prot_slopecurv */
1666
1667
1668 void makenewv(node *p)
1669 {
1670   /* Newton-Raphson algorithm improvement of a branch length */
1671   long it, ite;
1672   double y, yold=0, yorig, like, slope, curve, oldlike=0;
1673   boolean done, firsttime, better;
1674   node *q;
1675
1676   q = p->back;
1677   y = p->v;
1678   yorig = y;
1679   done = false;
1680   firsttime = true;
1681   it = 1;
1682   ite = 0;
1683   while ((it < iterations) && (ite < 20) && (!done)) {
1684     prot_slopecurv(p, y, &like, &slope, &curve);
1685     better = false;
1686     if (firsttime) {
1687       yold = y;
1688       oldlike = like;
1689       firsttime = false;
1690       better = true;
1691     } else {
1692       if (like > oldlike) {
1693         yold = y;
1694         oldlike = like;
1695         better = true;
1696         it++;
1697       }
1698     }
1699     if (better) {
1700       y = y + slope/fabs(curve);
1701       if (y < epsilon)
1702         y = epsilon;
1703     } else {
1704       if (fabs(y - yold) < epsilon)
1705         ite = 20;
1706       y = (y + (7 * yold)) / 8;
1707     }
1708     ite++;
1709     done = fabs(y-yold) < epsilon;
1710   }
1711   smoothed = (fabs(yold-yorig) < epsilon) && (yorig > 1000.0*epsilon);
1712   p->v = yold;
1713   q->v = yold;
1714   curtree.likelihood = oldlike;
1715 }  /* makenewv */
1716
1717
1718 void update(node *p)
1719 {
1720   if (!p->tip && !p->initialized)
1721     prot_nuview(p);
1722   if (!p->back->tip && !p->back->initialized)
1723     prot_nuview(p->back);
1724   if ((!usertree) || (usertree && !lngths) || p->iter) {
1725     makenewv(p);
1726     if ( smoothit ) {
1727       inittrav(p);
1728       inittrav(p->back);
1729     }
1730     else if ( inserting && !p->tip ) {
1731       p->next->initialized = false;
1732       p->next->next->initialized = false;
1733     }
1734   }
1735 }  /* update */
1736
1737
1738 void smooth(node *p)
1739 {
1740   long i, num_sibs;
1741   node *sib_ptr;
1742
1743   smoothed = false;
1744   update(p);
1745   if (p->tip)
1746     return;
1747
1748   num_sibs = count_sibs(p);
1749   sib_ptr  = p;
1750
1751   for (i=0; i < num_sibs; i++) {
1752     sib_ptr = sib_ptr->next;
1753
1754     if (polishing || (smoothit && !smoothed)) {
1755       smooth(sib_ptr->back);
1756       p->initialized = false;
1757       sib_ptr->initialized = false;
1758     }
1759   }
1760 }  /* smooth */
1761
1762
1763 void make_pmatrix(double **matrix, double **dmat, double **ddmat,
1764                         long derivative, double lz, double rat,
1765                         double *eigmat, double **probmat)
1766 {
1767   /* Computes the R matrix such that matrix[m][l] is the joint probability */
1768   /* of m and l.                                                           */
1769   /* Computes a P matrix such that matrix[m][l] is the conditional         */
1770   /* probability of m given l.  This is accomplished by dividing all terms */
1771   /* in the R matrix by freqaa[m], the frequency of l.                     */
1772
1773   long k, l, m;                 /* (l) original character state */
1774                                 /* (m) final    character state */
1775                                 /* (k) lambda counter           */
1776   double p0, p1, p2, q;
1777   double elambdat[20], delambdat[20], ddelambdat[20];
1778                                 /* exponential term for matrix  */
1779                                 /* and both derivative matrices */
1780   for (k = 0; k <= 19; k++) {
1781     elambdat[k] = exp(lz * rat * eigmat[k]);
1782     if(derivative != 0) {
1783         delambdat[k] = (elambdat[k] * rat * eigmat[k]);
1784         ddelambdat[k] = (delambdat[k] * rat * eigmat[k]);
1785     }
1786    }
1787   for (m = 0; m <= 19; m++) {
1788     for (l = 0; l <= 19; l++) {
1789       p0 = 0.0;
1790       p1 = 0.0;
1791       p2 = 0.0;
1792       for (k = 0; k <= 19; k++) {
1793         q = probmat[k][m] * probmat[k][l];
1794         p0 += (q * elambdat[k]);
1795         if(derivative !=0) {
1796           p1 += (q * delambdat[k]);
1797           p2 += (q * ddelambdat[k]);
1798         }
1799       }
1800       matrix[m][l] = p0 / freqaa[m];
1801       if(derivative != 0) {
1802         dmat[m][l] = p1 / freqaa[m];
1803         ddmat[m][l] = p2 / freqaa[m];
1804       }
1805     }
1806   }
1807 }  /* make_pmatrix */
1808
1809
1810 double prot_evaluate(node *p, boolean saveit)
1811 {
1812   contribarr tterm;
1813   double sum, sum2, sumc, y, prod4, prodl, frexm, sumterm, lterm;
1814   double **pmat;
1815   long i, j, k, l, m, lai;
1816   node *q;
1817   psitelike x1, x2;
1818
1819   sum = 0.0;
1820   q = p->back;
1821   y = p->v;
1822   for (j = 0; j < rcategs; j++)
1823     for (k = 0; k < categs; k++)
1824       make_pmatrix(pmatrices[0][j][k],NULL,NULL,0,y,tbl[j][k],eigmat,probmat);
1825   for (i = 0; i < endsite; i++) {
1826     k = category[alias[i]-1] - 1;
1827     for (j = 0; j < rcategs; j++) {
1828       memcpy(x1, p->protx[i][j], sizeof(psitelike));
1829       memcpy(x2, q->protx[i][j], sizeof(psitelike));
1830       prod4 = 0.0;
1831       pmat = pmatrices[0][j][k];
1832       for (m = 0; m <= 19; m++) {
1833         prodl = 0.0;
1834         for (l = 0; l <= 19; l++)
1835           prodl += (pmat[m][l] * x2[l]);
1836         frexm = x1[m] * freqaa[m];
1837         prod4 += (prodl * frexm);
1838       }
1839       tterm[j] = prod4;
1840     }
1841     sumterm = 0.0;
1842     for (j = 0; j < rcategs; j++)
1843       sumterm += probcat[j] * tterm[j];
1844     if (sumterm < 0.0)
1845         sumterm = 0.00000001;        /* ??? */
1846     lterm = log(sumterm) + p->underflows[i] + q->underflows[i];
1847     for (j = 0; j < rcategs; j++)
1848       clai[j] = tterm[j] / sumterm;
1849     memcpy(contribution[i], clai, rcategs*sizeof(double));
1850     if (saveit && !auto_ && usertree && (which <= shimotrees))
1851       l0gf[which - 1][i] = lterm;
1852     sum += aliasweight[i] * lterm;
1853   }
1854   for (j = 0; j < rcategs; j++)
1855     like[j] = 1.0;
1856   for (i = 0; i < sites; i++) {
1857     sumc = 0.0;
1858     for (k = 0; k < rcategs; k++)
1859       sumc += probcat[k] * like[k];
1860     sumc *= lambda;
1861     if ((ally[i] > 0) && (location[ally[i]-1] > 0)) {
1862       lai = location[ally[i] - 1];
1863       memcpy(clai, contribution[lai - 1], rcategs*sizeof(double));
1864       for (j = 0; j < rcategs; j++)
1865         nulike[j] = ((1.0 - lambda) * like[j] + sumc) * clai[j];
1866     } else {
1867       for (j = 0; j < rcategs; j++)
1868         nulike[j] = ((1.0 - lambda) * like[j] + sumc);
1869     }
1870     memcpy(like, nulike, rcategs*sizeof(double));
1871   }
1872   sum2 = 0.0;
1873   for (i = 0; i < rcategs; i++)
1874     sum2 += probcat[i] * like[i];
1875   sum += log(sum2);
1876   curtree.likelihood = sum;
1877   if (!saveit || auto_ || !usertree)
1878     return sum;
1879   if(which <= shimotrees)
1880     l0gl[which - 1] = sum;
1881   if (which == 1) {
1882     maxwhich = 1;
1883     maxlogl = sum;
1884     return sum;
1885   }
1886   if (sum > maxlogl) {
1887     maxwhich = which;
1888     maxlogl = sum;
1889   }
1890   return sum;
1891 }  /* prot_evaluate */
1892
1893
1894 void treevaluate()
1895 {
1896   /* evaluate a user tree */
1897   long i;
1898
1899   inittravtree(curtree.start);
1900   polishing = true;
1901   smoothit = true;
1902   for (i = 1; i <= smoothings * 4; i++)
1903     smooth (curtree.start);
1904   dummy = prot_evaluate(curtree.start, true);
1905 }  /* treevaluate */
1906
1907
1908 void promlcopy(tree *a, tree *b, long nonodes, long categs)
1909 {
1910   /* copy tree a to tree b */
1911   long i, j=0;
1912   node *p, *q;
1913
1914   for (i = 0; i < spp; i++) {
1915     prot_copynode(a->nodep[i], b->nodep[i], categs);
1916     if (a->nodep[i]->back) {
1917       if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1])
1918         b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1];
1919       else if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1]->next
1920 )
1921         b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next;
1922       else
1923         b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next->next;
1924     }
1925     else b->nodep[i]->back = NULL;
1926   }
1927   for (i = spp; i < nonodes; i++) {
1928     p = a->nodep[i];
1929     q = b->nodep[i];
1930     for (j = 1; j <= 3; j++) {
1931       prot_copynode(p, q, categs);
1932       if (p->back) {
1933         if (p->back == a->nodep[p->back->index - 1])
1934           q->back = b->nodep[p->back->index - 1];
1935         else if (p->back == a->nodep[p->back->index - 1]->next)
1936           q->back = b->nodep[p->back->index - 1]->next;
1937         else
1938           q->back = b->nodep[p->back->index - 1]->next->next;
1939       }
1940       else
1941         q->back = NULL;
1942       p = p->next;
1943       q = q->next;
1944     }
1945   }
1946   b->likelihood = a->likelihood;
1947   b->start = a->start;               /* start used in dnaml only */
1948   b->root = a->root;                 /* root used in dnamlk only */
1949 }  /* promlcopy */
1950
1951
1952 void proml_re_move(node **p, node **q)
1953 {
1954   /* remove p and record in q where it was */
1955   long i;
1956
1957   /** assumes bifurcations */
1958   *q = (*p)->next->back;
1959   hookup(*q, (*p)->next->next->back);
1960   (*p)->next->back = NULL;
1961   (*p)->next->next->back = NULL;
1962   (*q)->v += (*q)->back->v;
1963   (*q)->back->v = (*q)->v;
1964   if ( smoothit ) {
1965     inittrav((*q));
1966     inittrav((*q)->back);
1967     inittrav((*p)->back);
1968   }
1969   if ( smoothit )  {
1970     for ( i = 0 ; i < smoothings ; i++ ) {
1971       smooth(*q);
1972       smooth((*q)->back);
1973     }
1974   }
1975   else
1976     smooth(*q);
1977 }  /* proml_re_move */
1978
1979
1980 void insert_(node *p, node *q, boolean dooinit)
1981 {
1982   /* Insert q near p */
1983   long i, j, num_sibs;
1984   node *r, *sib_ptr;
1985
1986   r = p->next->next;
1987   hookup(r, q->back);
1988   hookup(p->next, q);
1989   q->v = 0.5 * q->v;
1990   q->back->v = q->v;
1991   r->v = q->v;
1992   r->back->v = r->v;
1993   p->initialized = false;
1994   if (dooinit) {
1995     inittrav(p);
1996     inittrav(q);
1997     inittrav(q->back);
1998   }
1999   i = 1;
2000   inserting = true;
2001   while (i <= smoothings) {
2002     smooth(p);
2003     if (!p->tip) {
2004       num_sibs = count_sibs(p);
2005       sib_ptr  = p;
2006       for (j=0; j < num_sibs; j++) {
2007         smooth(sib_ptr->next->back);
2008         sib_ptr = sib_ptr->next;
2009       }
2010     }
2011     i++;
2012   }
2013   inserting = false;
2014 }  /* insert_ */
2015
2016
2017 void addtraverse(node *p, node *q, boolean contin)
2018 {
2019   /* try adding p at q, proceed recursively through tree */
2020   long i, num_sibs;
2021   double like, vsave = 0;
2022   node *qback = NULL, *sib_ptr;
2023
2024   if (!smoothit) {
2025     vsave = q->v;
2026     qback = q->back;
2027   }
2028   insert_(p, q, false);
2029   like = prot_evaluate(p, false);
2030   if (like > bestyet || bestyet == UNDEFINED) {
2031     bestyet = like;
2032     if (smoothit) {
2033       addwhere = q;
2034       promlcopy(&curtree, &bestree, nonodes2, rcategs);
2035     }
2036     else
2037       qwhere = q;
2038     succeeded = true;
2039   }
2040   if (smoothit)
2041     promlcopy(&priortree, &curtree, nonodes2, rcategs);
2042   else {
2043     hookup (q, qback);
2044     q->v = vsave;
2045     q->back->v = vsave;
2046     curtree.likelihood = bestyet;
2047   }
2048   if (!q->tip && contin) {
2049     num_sibs = count_sibs(q);
2050     if (q == curtree.start)
2051       num_sibs++;
2052     sib_ptr  = q;
2053     for (i=0; i < num_sibs; i++) {
2054       addtraverse(p, sib_ptr->next->back, contin);
2055       sib_ptr = sib_ptr->next;
2056     }
2057   }
2058 }  /* addtraverse */
2059
2060
2061 void globrearrange()
2062 {
2063   /* does global rearrangements */
2064   tree globtree;
2065   tree oldtree;
2066   int i,j,k,l,num_sibs,num_sibs2;
2067   node *where,*sib_ptr,*sib_ptr2;
2068   double oldbestyet = curtree.likelihood;
2069   int success = false;
2070
2071   alloctree(&globtree.nodep,nonodes2,0);
2072   alloctree(&oldtree.nodep,nonodes2,0);
2073   setuptree2(globtree);
2074   setuptree2(oldtree);
2075   prot_allocx(nonodes2, rcategs, globtree.nodep, 0);
2076   prot_allocx(nonodes2, rcategs, oldtree.nodep, 0);
2077   promlcopy(&curtree,&globtree,nonodes2,rcategs);
2078   promlcopy(&curtree,&oldtree,nonodes2,rcategs);
2079   bestyet = curtree.likelihood;
2080   for ( i = spp ; i < nonodes2 ; i++ ) {
2081     num_sibs = count_sibs(curtree.nodep[i]);
2082     sib_ptr  = curtree.nodep[i];
2083     if ( (i - spp) % (( nonodes2 / 72 ) + 1 ) == 0 )
2084       putchar('.');
2085     fflush(stdout);
2086     for ( j = 0 ; j <= num_sibs ; j++ ) {
2087       proml_re_move(&sib_ptr,&where);
2088       promlcopy(&curtree,&priortree,nonodes2,rcategs);
2089       qwhere = where;
2090
2091       if (where->tip) {
2092         promlcopy(&oldtree,&curtree,nonodes2,rcategs);
2093         promlcopy(&oldtree,&bestree,nonodes2,rcategs);
2094         sib_ptr=sib_ptr->next;
2095         continue;
2096       }
2097       else num_sibs2 = count_sibs(where);
2098       sib_ptr2 = where;
2099       for ( k = 0 ; k < num_sibs2 ; k++ ) {
2100         addwhere = NULL;
2101         addtraverse(sib_ptr,sib_ptr2->back,true);
2102         if ( !smoothit ) {
2103           if (succeeded && qwhere != where && qwhere != where->back) {
2104             insert_(sib_ptr,qwhere,true);
2105             smoothit = true;
2106             for (l = 1; l<=smoothings; l++) {
2107               smooth (where);
2108               smooth (where->back);
2109             }
2110             smoothit = false;
2111             success = true;
2112             promlcopy(&curtree,&globtree,nonodes2,rcategs);
2113             promlcopy(&priortree,&curtree,nonodes2,rcategs);
2114           }
2115         }
2116         else if ( addwhere && where != addwhere && where->back != addwhere
2117               && bestyet > globtree.likelihood) {
2118             promlcopy(&bestree,&globtree,nonodes2,rcategs);
2119             success = true;
2120         }
2121         sib_ptr2 = sib_ptr2->next;
2122       }
2123       promlcopy(&oldtree,&curtree,nonodes2,rcategs);
2124       promlcopy(&oldtree,&bestree,nonodes2,rcategs);
2125       sib_ptr = sib_ptr->next;
2126     }
2127   }
2128   promlcopy(&globtree,&curtree,nonodes2,rcategs);
2129   promlcopy(&globtree,&bestree,nonodes2,rcategs);
2130   if (success && globtree.likelihood > oldbestyet)  {
2131     succeeded = true;
2132   }
2133   else  {
2134     succeeded = false;
2135   }
2136   bestyet = globtree.likelihood;
2137   prot_freex(nonodes2,oldtree.nodep);
2138   prot_freex(nonodes2,globtree.nodep);
2139   freetree2(globtree.nodep,nonodes2);
2140   freetree2(oldtree.nodep,nonodes2);
2141 } /* globrearrange */
2142
2143
2144 void freelrsaves()
2145 {
2146   long i,j;
2147   for ( i = 0 ; i < NLRSAVES ; i++ ) {
2148     for (j = 0; j < oldendsite; j++)
2149       free(lrsaves[i]->protx[j]);
2150     free(lrsaves[i]->protx);
2151     free(lrsaves[i]->underflows);
2152     free(lrsaves[i]);
2153   }
2154   free(lrsaves);
2155 }
2156
2157
2158 void alloclrsaves()
2159 {
2160   long i,j;
2161   lrsaves = Malloc(NLRSAVES * sizeof(node*));
2162   oldendsite = endsite;
2163   for ( i = 0 ; i < NLRSAVES ; i++ ) {
2164     lrsaves[i] = Malloc(sizeof(node));
2165     lrsaves[i]->protx = Malloc(endsite*sizeof(ratelike));
2166     lrsaves[i]->underflows = Malloc(endsite * sizeof (double));
2167     for (j = 0; j < endsite; j++)
2168       lrsaves[i]->protx[j]  = (pratelike)Malloc(rcategs*sizeof(psitelike));
2169   }
2170 } /* alloclrsaves */
2171
2172
2173 void rearrange(node *p, node *pp)
2174 {
2175   /* rearranges the tree locally moving pp around near p */
2176   long i, num_sibs;
2177   node *q, *r, *sib_ptr;
2178   node *rnb, *rnnb;
2179
2180   if (!p->tip && !p->back->tip) {
2181     curtree.likelihood = bestyet;
2182     if (p->back->next != pp)
2183       r = p->back->next;
2184     else
2185       r = p->back->next->next;
2186     /* assumes bifurcations? */
2187     if (!smoothit) {
2188       rnb = r->next->back;
2189       rnnb = r->next->next->back;
2190       prot_copynode(r,lrsaves[0],categs);
2191       prot_copynode(r->next,lrsaves[1],categs);
2192       prot_copynode(r->next->next,lrsaves[2],categs);
2193       prot_copynode(p->next,lrsaves[3],categs);
2194       prot_copynode(p->next->next,lrsaves[4],categs);
2195     }
2196     else
2197       promlcopy(&curtree, &bestree, nonodes2, rcategs);
2198     proml_re_move(&r, &q);
2199     if (smoothit)
2200       promlcopy(&curtree, &priortree, nonodes2, rcategs);
2201     else
2202       qwhere = q;
2203     num_sibs = count_sibs (p);
2204     sib_ptr  = p;
2205     for (i=0; i < num_sibs; i++) {
2206       sib_ptr = sib_ptr->next;
2207       addtraverse(r, sib_ptr->back, false);
2208     }
2209     if (smoothit)
2210       promlcopy(&bestree, &curtree, nonodes2, rcategs);
2211     else {
2212       if (qwhere == q) {
2213         hookup(rnb,r->next);
2214         hookup(rnnb,r->next->next);
2215         prot_copynode(lrsaves[0],r,categs);
2216         prot_copynode(lrsaves[1],r->next,categs);
2217         prot_copynode(lrsaves[2],r->next->next,categs);
2218         prot_copynode(lrsaves[3],p->next,categs);
2219         prot_copynode(lrsaves[4],p->next->next,categs);
2220         rnb->v = r->next->v;
2221         rnnb->v = r->next->next->v;
2222         r->back->v = r->v;
2223         curtree.likelihood = bestyet;
2224       }
2225       else {
2226         insert_(r, qwhere, true);
2227         smoothit = true;
2228         for (i = 1; i<=smoothings; i++) {
2229           smooth(r);
2230           smooth(r->back);
2231         }
2232         smoothit = false;
2233         promlcopy(&curtree, &bestree, nonodes2, rcategs);
2234       }
2235     }
2236   }
2237   if (!p->tip) {
2238     num_sibs = count_sibs(p);
2239     if (p == curtree.start)
2240       num_sibs++;
2241     sib_ptr  = p;
2242     for (i=0; i < num_sibs; i++) {
2243       sib_ptr = sib_ptr->next;
2244       rearrange(sib_ptr->back, p);
2245     }
2246   }
2247 }  /* rearrange */
2248
2249
2250 void proml_coordinates(node *p, double lengthsum, long *tipy,
2251                         double *tipmax)
2252 {
2253   /* establishes coordinates of nodes */
2254   node *q, *first, *last;
2255   double xx;
2256
2257   if (p->tip) {
2258     p->xcoord = (long)(over * lengthsum + 0.5);
2259     p->ycoord = (*tipy);
2260     p->ymin = (*tipy);
2261     p->ymax = (*tipy);
2262     (*tipy) += down;
2263     if (lengthsum > (*tipmax))
2264       (*tipmax) = lengthsum;
2265     return;
2266   }
2267   q = p->next;
2268   do {
2269     xx = q->v;
2270     if (xx > 100.0)
2271       xx = 100.0;
2272     proml_coordinates(q->back, lengthsum + xx, tipy,tipmax);
2273     q = q->next;
2274   } while ((p == curtree.start || p != q) &&
2275            (p != curtree.start || p->next != q));
2276   first = p->next->back;
2277   q = p;
2278   while (q->next != p)
2279     q = q->next;
2280   last = q->back;
2281   p->xcoord = (long)(over * lengthsum + 0.5);
2282   if (p == curtree.start)
2283     p->ycoord = p->next->next->back->ycoord;
2284   else
2285     p->ycoord = (first->ycoord + last->ycoord) / 2;
2286   p->ymin = first->ymin;
2287   p->ymax = last->ymax;
2288 }  /* proml_coordinates */
2289
2290
2291 void proml_printree()
2292 {
2293   /* prints out diagram of the tree2 */
2294   long tipy;
2295   double scale, tipmax;
2296   long i;
2297
2298   if (!treeprint)
2299     return;
2300   putc('\n', outfile);
2301   tipy = 1;
2302   tipmax = 0.0;
2303   proml_coordinates(curtree.start, 0.0, &tipy, &tipmax);
2304   scale = 1.0 / (long)(tipmax + 1.000);
2305   for (i = 1; i <= (tipy - down); i++)
2306     drawline2(i, scale, curtree);
2307   putc('\n', outfile);
2308 }  /* proml_printree */
2309
2310
2311 void sigma(node *p, double *sumlr, double *s1, double *s2)
2312 {
2313   /* compute standard deviation */
2314   double tt, aa, like, slope, curv;
2315
2316   prot_slopecurv(p, p->v, &like, &slope, &curv);
2317   tt = p->v;
2318   p->v = epsilon;
2319   p->back->v = epsilon;
2320   aa = prot_evaluate(p, false);
2321   p->v = tt;
2322   p->back->v = tt;
2323   (*sumlr) = prot_evaluate(p, false) - aa;
2324   if (curv < -epsilon) {
2325     (*s1) = p->v + (-slope - sqrt(slope * slope -  3.841 * curv)) / curv;
2326     (*s2) = p->v + (-slope + sqrt(slope * slope -  3.841 * curv)) / curv;
2327   }
2328   else {
2329     (*s1) = -1.0;
2330     (*s2) = -1.0;
2331   }
2332 }  /* sigma */
2333
2334
2335 void describe(node *p)
2336 {
2337   /* print out information for one branch */
2338   long i, num_sibs;
2339   node *q, *sib_ptr;
2340   double sumlr, sigma1, sigma2;
2341
2342   if (!p->tip && !p->initialized)
2343     prot_nuview(p);
2344   if (!p->back->tip && !p->back->initialized)
2345     prot_nuview(p->back);
2346   q = p->back;
2347   if (q->tip) {
2348     fprintf(outfile, " ");
2349     for (i = 0; i < nmlngth; i++)
2350       putc(nayme[q->index-1][i], outfile);
2351     fprintf(outfile, "    ");
2352   } else
2353     fprintf(outfile, "  %4ld          ", q->index - spp);
2354   if (p->tip) {
2355     for (i = 0; i < nmlngth; i++)
2356       putc(nayme[p->index-1][i], outfile);
2357   } else
2358     fprintf(outfile, "%4ld      ", p->index - spp);
2359   fprintf(outfile, "%15.5f", q->v);
2360   if (!usertree || (usertree && !lngths) || p->iter) {
2361     sigma(q, &sumlr, &sigma1, &sigma2);
2362     if (sigma1 <= sigma2)
2363       fprintf(outfile, "     (     zero,    infinity)");
2364     else {
2365       fprintf(outfile, "     (");
2366       if (sigma2 <= 0.0)
2367         fprintf(outfile, "     zero");
2368       else
2369         fprintf(outfile, "%9.5f", sigma2);
2370       fprintf(outfile, ",%12.5f", sigma1);
2371       putc(')', outfile);
2372       }
2373     if (sumlr > 1.9205)
2374       fprintf(outfile, " *");
2375     if (sumlr > 2.995)
2376       putc('*', outfile);
2377     }
2378   putc('\n', outfile);
2379   if (!p->tip) {
2380     num_sibs = count_sibs(p);
2381     sib_ptr  = p;
2382     for (i=0; i < num_sibs; i++) {
2383       sib_ptr = sib_ptr->next;
2384       describe(sib_ptr->back);
2385     }
2386   }
2387 }  /* describe */
2388
2389
2390 void prot_reconstr(node *p, long n)
2391 {
2392   /* reconstruct and print out acid at site n+1 at node p */
2393   long i, j, k, first, num_sibs = 0;
2394   double f, sum, xx[20];
2395   node *q = NULL;
2396
2397   if (p->tip)
2398     putc(y[p->index-1][n], outfile);
2399   else {
2400     num_sibs = count_sibs(p);
2401     if ((ally[n] == 0) || (location[ally[n]-1] == 0))
2402       putc('.', outfile);
2403     else {
2404       j = location[ally[n]-1] - 1;
2405       sum = 0;
2406       for (i = 0; i <= 19; i++) {
2407         f = p->protx[j][mx-1][i];
2408         if (!p->tip) {
2409           q = p;
2410           for (k = 0; k < num_sibs; k++) {
2411             q = q->next;
2412             f *= q->protx[j][mx-1][i];
2413           }
2414         }
2415         f = sqrt(f);
2416         xx[i] = f * freqaa[i];
2417         sum += xx[i];
2418       }
2419       for (i = 0; i <= 19; i++)
2420         xx[i] /= sum;
2421       first = 0;
2422       for (i = 0; i <= 19; i++)
2423         if (xx[i] > xx[first])
2424           first = i;
2425       if (xx[first] > 0.95)
2426         putc(aachar[first], outfile);
2427       else
2428         putc(tolower(aachar[first]), outfile);
2429       if (rctgry && rcategs > 1)
2430         mx = mp[n][mx - 1];
2431       else
2432         mx = 1;
2433     }
2434   }
2435 } /* prot_reconstr */
2436
2437
2438 void rectrav(node *p, long m, long n)
2439 {
2440   /* print out segment of reconstructed sequence for one branch */
2441   long i;
2442
2443   putc(' ', outfile);
2444   if (p->tip) {
2445     for (i = 0; i < nmlngth; i++)
2446       putc(nayme[p->index-1][i], outfile);
2447   } else
2448     fprintf(outfile, "%4ld      ", p->index - spp);
2449   fprintf(outfile, "  ");
2450   mx = mx0;
2451   for (i = m; i <= n; i++) {
2452     if ((i % 10 == 0) && (i != m))
2453       putc(' ', outfile);
2454     prot_reconstr(p, i);
2455   }
2456   putc('\n', outfile);
2457   if (!p->tip) {
2458     rectrav(p->next->back, m, n);
2459     rectrav(p->next->next->back, m, n);
2460   }
2461   mx1 = mx;
2462 }  /* rectrav */
2463
2464
2465 void summarize()
2466 {
2467   /* print out branch length information and node numbers */
2468   long i, j, mm, num_sibs;
2469   double mode, sum;
2470   double like[maxcategs],nulike[maxcategs];
2471   double **marginal;
2472   node   *sib_ptr;
2473
2474   if (!treeprint)
2475     return;
2476   fprintf(outfile, "\nremember: ");
2477   if (outgropt)
2478     fprintf(outfile, "(although rooted by outgroup) ");
2479   fprintf(outfile, "this is an unrooted tree!\n\n");
2480   fprintf(outfile, "Ln Likelihood = %11.5f\n", curtree.likelihood);
2481   fprintf(outfile, "\n Between        And            Length");
2482   if (!(usertree && lngths && haslengths))
2483     fprintf(outfile, "      Approx. Confidence Limits");
2484   fprintf(outfile, "\n");
2485   fprintf(outfile, " -------        ---            ------");
2486   if (!(usertree && lngths && haslengths))
2487     fprintf(outfile, "      ------- ---------- ------");
2488   fprintf(outfile, "\n\n");
2489   for (i = spp; i < nonodes2; i++) {
2490     /* So this works with arbitrary multifurcations */
2491     if (curtree.nodep[i]) {
2492       num_sibs = count_sibs (curtree.nodep[i]);
2493       sib_ptr  = curtree.nodep[i];
2494       for (j = 0; j < num_sibs; j++) {
2495         sib_ptr->initialized = false;
2496         sib_ptr              = sib_ptr->next;
2497       }
2498     }
2499   }
2500
2501   describe(curtree.start->back);
2502
2503   /* So this works with arbitrary multifurcations */
2504   num_sibs = count_sibs(curtree.start);
2505   sib_ptr  = curtree.start;
2506   for (i=0; i < num_sibs; i++) {
2507     sib_ptr = sib_ptr->next;
2508     describe(sib_ptr->back);
2509   }
2510
2511   fprintf(outfile, "\n");
2512   if (!(usertree && lngths && haslengths)) {
2513     fprintf(outfile, "     *  = significantly positive, P < 0.05\n");
2514     fprintf(outfile, "     ** = significantly positive, P < 0.01\n\n");
2515   }
2516   dummy = prot_evaluate(curtree.start, false);
2517   if (rctgry && rcategs > 1) {
2518     for (i = 0; i < rcategs; i++)
2519       like[i] = 1.0;
2520     for (i = sites - 1; i >= 0; i--) {
2521       sum = 0.0;
2522       for (j = 0; j < rcategs; j++) {
2523         nulike[j] = (1.0 - lambda + lambda * probcat[j]) * like[j];
2524         mp[i][j] = j + 1;
2525         for (k = 1; k <= rcategs; k++) {
2526           if (k != j + 1) {
2527             if (lambda * probcat[k - 1] * like[k - 1] > nulike[j]) {
2528               nulike[j] = lambda * probcat[k - 1] * like[k - 1];
2529               mp[i][j] = k;
2530             }
2531           }
2532         }
2533         if ((ally[i] > 0) && (location[ally[i]-1] > 0))
2534           nulike[j] *= contribution[location[ally[i] - 1] - 1][j];
2535         sum += nulike[j];
2536       }
2537       for (j = 0; j < rcategs; j++)
2538         nulike[j] /= sum;
2539       memcpy(like, nulike, rcategs * sizeof(double));
2540     }
2541     mode = 0.0;
2542     mx = 1;
2543     for (i = 1; i <= rcategs; i++) {
2544       if (probcat[i - 1] * like[i - 1] > mode) {
2545         mx = i;
2546         mode = probcat[i - 1] * like[i - 1];
2547       }
2548     }
2549     mx0 = mx;
2550     fprintf(outfile,
2551  "Combination of categories that contributes the most to the likelihood:\n\n");
2552     for (i = 1; i <= nmlngth + 3; i++)
2553       putc(' ', outfile);
2554     for (i = 1; i <= sites; i++) {
2555       fprintf(outfile, "%ld", mx);
2556       if (i % 10 == 0)
2557         putc(' ', outfile);
2558       if (i % 60 == 0 && i != sites) {
2559         putc('\n', outfile);
2560         for (j = 1; j <= nmlngth + 3; j++)
2561           putc(' ', outfile);
2562       }
2563       mx = mp[i - 1][mx - 1];
2564     }
2565     fprintf(outfile, "\n\n");
2566     marginal = (double **) Malloc(sites*sizeof(double *));
2567     for (i = 0; i < sites; i++)
2568       marginal[i] = (double *) Malloc(rcategs*sizeof(double));
2569     for (i = 0; i < rcategs; i++)
2570       like[i] = 1.0;
2571     for (i = sites - 1; i >= 0; i--) {
2572       sum = 0.0;
2573       for (j = 0; j < rcategs; j++) {
2574         nulike[j] = (1.0 - lambda + lambda * probcat[j]) * like[j];
2575         for (k = 1; k <= rcategs; k++) {
2576           if (k != j + 1)
2577               nulike[j] += lambda * probcat[k - 1] * like[k - 1];
2578         }
2579         if ((ally[i] > 0) && (location[ally[i]-1] > 0))
2580           nulike[j] *= contribution[location[ally[i] - 1] - 1][j];
2581         sum += nulike[j];
2582       }
2583       for (j = 0; j < rcategs; j++) {
2584         nulike[j] /= sum;
2585         marginal[i][j] = nulike[j];
2586       }
2587       memcpy(like, nulike, rcategs * sizeof(double));
2588     }
2589     for (i = 0; i < rcategs; i++)
2590       like[i] = 1.0;
2591     for (i = 0; i < sites; i++) {
2592       sum = 0.0;
2593       for (j = 0; j < rcategs; j++) {
2594         nulike[j] = (1.0 - lambda + lambda * probcat[j]) * like[j];
2595         for (k = 1; k <= rcategs; k++) {
2596           if (k != j + 1)
2597               nulike[j] += lambda * probcat[k - 1] * like[k - 1];
2598         }
2599         marginal[i][j] *= like[j] * probcat[j];
2600         sum += nulike[j];
2601       }
2602       for (j = 0; j < rcategs; j++)
2603         nulike[j] /= sum;
2604       memcpy(like, nulike, rcategs * sizeof(double));
2605       sum = 0.0;
2606       for (j = 0; j < rcategs; j++)
2607         sum += marginal[i][j];
2608       for (j = 0; j < rcategs; j++)
2609         marginal[i][j] /= sum;
2610     }
2611     fprintf(outfile, "Most probable category at each site if > 0.95");
2612     fprintf(outfile, " probability (\".\" otherwise)\n\n");
2613     for (i = 1; i <= nmlngth + 3; i++)
2614       putc(' ', outfile);
2615     for (i = 0; i < sites; i++) {
2616       sum = 0.0;
2617       for (j = 0; j < rcategs; j++)
2618         if (marginal[i][j] > sum) {
2619           sum = marginal[i][j];
2620           mm = j;
2621         }
2622         if (sum >= 0.95)
2623         fprintf(outfile, "%ld", mm+1);
2624       else
2625         putc('.', outfile);
2626       if ((i+1) % 60 == 0) {
2627         if (i != 0) {
2628           putc('\n', outfile);
2629           for (j = 1; j <= nmlngth + 3; j++)
2630             putc(' ', outfile);
2631         }
2632       }
2633       else if ((i+1) % 10 == 0)
2634         putc(' ', outfile);
2635     }
2636     putc('\n', outfile);
2637     for (i = 0; i < sites; i++)
2638       free(marginal[i]);
2639     free(marginal);
2640   }
2641   putc('\n', outfile);
2642   if (hypstate) {
2643     fprintf(outfile, "Probable sequences at interior nodes:\n\n");
2644     fprintf(outfile, "  node       ");
2645     for (i = 0; (i < 13) && (i < ((sites + (sites-1)/10 - 39) / 2)); i++)
2646       putc(' ', outfile);
2647     fprintf(outfile, "Reconstructed sequence (caps if > 0.95)\n\n");
2648     if (!rctgry || (rcategs == 1))
2649       mx0 = 1;
2650     for (i = 0; i < sites; i += 60) {
2651       k = i + 59;
2652       if (k >= sites)
2653         k = sites - 1;
2654       rectrav(curtree.start, i, k);
2655       rectrav(curtree.start->back, i, k);
2656       putc('\n', outfile);
2657       mx0 = mx1;
2658     }
2659   }
2660 }  /* summarize */
2661
2662
2663 void initpromlnode(node **p, node **grbg, node *q, long len, long nodei,
2664                         long *ntips, long *parens, initops whichinit,
2665                         pointarray treenode, pointarray nodep, Char *str,
2666                         Char *ch, FILE *intree)
2667 {
2668   /* initializes a node */
2669   boolean minusread;
2670   double valyew, divisor;
2671
2672   switch (whichinit) {
2673   case bottom:
2674     gnu(grbg, p);
2675     (*p)->index = nodei;
2676     (*p)->tip = false;
2677     malloc_ppheno((*p), endsite, rcategs);
2678     nodep[(*p)->index - 1] = (*p);
2679     break;
2680   case nonbottom:
2681     gnu(grbg, p);
2682     malloc_ppheno(*p, endsite, rcategs);
2683     (*p)->index = nodei;
2684     break;
2685   case tip:
2686     match_names_to_data(str, nodep, p, spp);
2687     break;
2688   case iter:
2689     (*p)->initialized = false;
2690     (*p)->v = initialv;
2691     (*p)->iter = true;
2692     if ((*p)->back != NULL){
2693       (*p)->back->iter = true;
2694       (*p)->back->v = initialv;
2695       (*p)->back->initialized = false;
2696     }
2697     break;
2698   case length:
2699     processlength(&valyew, &divisor, ch, &minusread, intree, parens);
2700     (*p)->v = valyew / divisor;
2701     (*p)->iter = false;
2702     if ((*p)->back != NULL) {
2703       (*p)->back->v = (*p)->v;
2704       (*p)->back->iter = false;
2705     }
2706     break;
2707   case hsnolength:
2708     haslengths = false;
2709     break;
2710   default:      /* cases hslength, treewt, unittrwt */
2711     break;      /* should never occur               */
2712   }
2713 } /* initpromlnode */
2714
2715
2716 void dnaml_treeout(node *p)
2717 {
2718   /* write out file with representation of final tree2 */
2719   /* Only works for bifurcations! */
2720   long i, n, w;
2721   Char c;
2722   double x;
2723   node *q;
2724   boolean inloop;
2725
2726   if (p->tip) {
2727     n = 0;
2728     for (i = 1; i <= nmlngth; i++) {
2729       if (nayme[p->index-1][i - 1] != ' ')
2730         n = i;
2731     }
2732     for (i = 0; i < n; i++) {
2733       c = nayme[p->index-1][i];
2734       if (c == ' ')
2735         c = '_';
2736       putc(c, outtree);
2737     }
2738     col += n;
2739   } else {
2740     putc('(', outtree);
2741     col++;
2742
2743     inloop = false;
2744     q = p->next;
2745     do  {
2746       if (inloop) {
2747         putc(',', outtree);
2748         col++;
2749         if (col > 45) {
2750           putc('\n', outtree);
2751           col = 0;
2752         }
2753       }
2754       inloop = true;
2755       dnaml_treeout(q->back);
2756       q = q->next;
2757     } while ((p == curtree.start || p != q) &&
2758              (p != curtree.start || p->next != q));
2759
2760     putc(')', outtree);
2761     col++;
2762   }
2763   x = p->v;
2764   if (x > 0.0)
2765     w = (long)(0.43429448222 * log(x));
2766   else if (x == 0.0)
2767     w = 0;
2768   else
2769     w = (long)(0.43429448222 * log(-x)) + 1;
2770   if (w < 0)
2771     w = 0;
2772   if (p == curtree.start)
2773     fprintf(outtree, ";\n");
2774   else {
2775     fprintf(outtree, ":%*.5f", (int)(w + 7), x);
2776     col += w + 8;
2777   }
2778 }  /* dnaml_treeout */
2779
2780
2781 void buildnewtip(long m, tree *tr)
2782 {
2783   node *p;
2784
2785   p = tr->nodep[nextsp + spp - 3];
2786   hookup(tr->nodep[m - 1], p);
2787   p->v = initialv;
2788   p->back->v = initialv;
2789 }  /* buildnewtip */
2790
2791
2792 void buildsimpletree(tree *tr)
2793 {
2794   hookup(tr->nodep[enterorder[0] - 1], tr->nodep[enterorder[1] - 1]);
2795   tr->nodep[enterorder[0] - 1]->v = 1.0;
2796   tr->nodep[enterorder[0] - 1]->back->v = 1.0;
2797   tr->nodep[enterorder[1] - 1]->v = 1.0;
2798   tr->nodep[enterorder[1] - 1]->back->v = 1.0;
2799   buildnewtip(enterorder[2], tr);
2800   insert_(tr->nodep[enterorder[2] - 1]->back,
2801           tr->nodep[enterorder[0] - 1], false);
2802 }  /* buildsimpletree */
2803
2804
2805 void free_all_protx (long nonodes, pointarray treenode)
2806 {
2807   /* used in proml */
2808   long i, j, k;
2809   node *p;
2810
2811   /* Zero thru spp are tips, */
2812   for (i = 0; i < spp; i++) {
2813     for (j = 0; j < endsite; j++)
2814       free(treenode[i]->protx[j]);
2815     free(treenode[i]->protx);
2816   }
2817
2818   /* The rest are rings (i.e. triads) */
2819   for (i = spp; i < nonodes; i++) {
2820     if (treenode[i] != NULL) {
2821       p = treenode[i];
2822       do {
2823         for (k = 0; k < endsite; k++)
2824           free(p->protx[k]);
2825         free(p->protx);
2826         p = p->next;
2827       } while (p != treenode[i]);
2828     }
2829   }
2830 }  /* free_all_protx */
2831
2832 void proml_unroot(node* root, node** nodep, long nonodes)
2833 {
2834   node *r,*q,*tmpnode;
2835   double newl;
2836   long i;
2837   long numsibs;
2838
2839   numsibs = count_sibs(root);
2840
2841   if ( numsibs > 2 ) {
2842     q = root;
2843     r = root;
2844     while (!(q->next == root))
2845       q = q->next;
2846     q->next = root->next;
2847     root = q;
2848     for(i=0 ; i < endsite ; i++){
2849       free(r->protx[i]);
2850       r->protx[i] = NULL;
2851     }
2852     free(r->protx);
2853     r->protx = NULL;
2854     chucktreenode(&grbg, r);
2855     curtree.nodep[spp] = q;
2856   } else if ( root->next->next->next == root) {
2857     newl = root->next->oldlen + root->next->next->oldlen;
2858     root->next->back->oldlen = newl;
2859     root->next->next->back->oldlen = newl;
2860
2861     newl = root->next->v + root->next->next->v;
2862     root->next->back->v = newl;
2863     root->next->next->back->v = newl;
2864
2865     root->next->back->back=root->next->next->back;
2866     root->next->next->back->back = root->next->back;
2867     while ( root->index != nonodes ) {
2868       tmpnode = nodep[ root->index ];
2869       nodep[root->index] = root;
2870       root->index++;
2871       root->next->index++;
2872       root->next->next->index++;
2873       nodep[root->index - 2] = tmpnode;
2874       tmpnode->index--;
2875       tmpnode->next->index--;
2876       tmpnode->next->next->index--;
2877     }
2878     nodep[nonodes -1] = NULL;
2879     for(i=0 ; i < endsite ; i++){
2880       free(root->protx[i]);
2881       free(root->next->protx[i]);
2882       free(root->next->next->protx[i]);
2883       root->protx[i] = NULL;
2884       root->next->protx[i] = NULL;
2885       root->next->next->protx[i] = NULL;
2886     }
2887     free(root->protx);
2888     free(root->next->protx);
2889     free(root->next->next->protx);
2890
2891     chucktreenode(&grbg,root->next->next);
2892     chucktreenode(&grbg,root->next);
2893     chucktreenode(&grbg,root);
2894
2895   }
2896 }
2897
2898
2899
2900 void maketree()
2901 {
2902   long i, j;
2903   boolean dummy_first, goteof;
2904   pointarray dummy_treenode=NULL;
2905   long nextnode;
2906   node *root;
2907
2908   prot_inittable();
2909
2910   if (usertree) {
2911     openfile(&intree,INTREE,"input tree file", "r",progname,intreename);
2912     numtrees = countsemic(&intree);
2913     if(numtrees > MAXSHIMOTREES)
2914       shimotrees = MAXSHIMOTREES;
2915     else
2916       shimotrees = numtrees;
2917     if (numtrees > 2)
2918       initseed(&inseed, &inseed0, seed);
2919     l0gl = (double *) Malloc(shimotrees * sizeof(double));
2920     l0gf = (double **) Malloc(shimotrees * sizeof(double *));
2921     for (i=0; i < shimotrees; ++i)
2922       l0gf[i] = (double *) Malloc(endsite * sizeof(double));
2923     if (treeprint) {
2924       fprintf(outfile, "User-defined tree");
2925       if (numtrees > 1)
2926         putc('s', outfile);
2927       fprintf(outfile, ":\n\n");
2928     }
2929     which = 1;
2930
2931     /* This taken out of tree read, used to be [spp-1], but referring
2932        to [0] produces output identical to what the pre-modified dnaml
2933        produced. */
2934
2935     while (which <= numtrees) {
2936
2937       /* These initializations required each time through the loop
2938          since multiple trees require re-initialization */
2939       haslengths = true;
2940       nextnode         = 0;
2941       dummy_first      = true;
2942       goteof           = false;
2943
2944       treeread(intree, &root, dummy_treenode, &goteof, &dummy_first,
2945                       curtree.nodep, &nextnode, &haslengths, &grbg,
2946                       initpromlnode,false,nonodes2);
2947       proml_unroot(root,curtree.nodep,nonodes2);
2948       if (goteof && (which <= numtrees)) {
2949         /* if we hit the end of the file prematurely */
2950         printf ("\n");
2951         printf ("ERROR: trees missing at end of file.\n");
2952         printf ("\tExpected number of trees:\t\t%ld\n", numtrees);
2953         printf ("\tNumber of trees actually in file:\t%ld.\n\n", which - 1);
2954         exxit (-1);
2955       }
2956
2957       curtree.start = curtree.nodep[0]->back;
2958       if ( outgropt )
2959         curtree.start = curtree.nodep[outgrno - 1]->back;
2960
2961       treevaluate();
2962       proml_printree();
2963       summarize();
2964       if (trout) {
2965         col = 0;
2966         dnaml_treeout(curtree.start);
2967       }
2968       if(which < numtrees){
2969         prot_freex_notip(nextnode, curtree.nodep);
2970         gdispose(curtree.start, &grbg, curtree.nodep);
2971       } else nonodes2 = nextnode;
2972       which++;
2973     }
2974     FClose(intree);
2975     putc('\n', outfile);
2976     if (!auto_ && numtrees > 1 && weightsum > 1 )
2977       standev2(numtrees, maxwhich, 0, endsite-1, maxlogl,
2978                l0gl, l0gf, aliasweight, seed);
2979   } else {
2980     /* If there's no input user tree, */
2981     for (i = 1; i <= spp; i++)
2982       enterorder[i - 1] = i;
2983     if (jumble)
2984       randumize(seed, enterorder);
2985     if (progress) {
2986       printf("\nAdding species:\n");
2987       writename(0, 3, enterorder);
2988 #ifdef WIN32
2989       phyFillScreenColor();
2990 #endif
2991     }
2992     nextsp = 3;
2993     polishing = false;
2994     buildsimpletree(&curtree);
2995     curtree.start = curtree.nodep[enterorder[0] - 1]->back;
2996     smoothit = improve;
2997     nextsp = 4;
2998     while (nextsp <= spp) {
2999       buildnewtip(enterorder[nextsp - 1], &curtree);
3000       bestyet = UNDEFINED;
3001       if (smoothit)
3002         promlcopy(&curtree, &priortree, nonodes2, rcategs);
3003       addtraverse(curtree.nodep[enterorder[nextsp - 1] - 1]->back,
3004                   curtree.start, true);
3005       if (smoothit)
3006         promlcopy(&bestree, &curtree, nonodes2, rcategs);
3007       else {
3008         insert_(curtree.nodep[enterorder[nextsp - 1] - 1]->back, qwhere, true);
3009         smoothit = true;
3010         for (i = 1; i<=smoothings; i++) {
3011           smooth(curtree.start);
3012           smooth(curtree.start->back);
3013         }
3014         smoothit = false;
3015         promlcopy(&curtree, &bestree, nonodes2, rcategs);
3016         bestyet = curtree.likelihood;
3017       }
3018       if (progress) {
3019         writename(nextsp - 1, 1, enterorder);
3020 #ifdef WIN32
3021         phyFillScreenColor();
3022 #endif
3023       }
3024       if (global && nextsp == spp && progress) {
3025         printf("Doing global rearrangements\n");
3026         printf("  !");
3027         for (j = spp ; j < nonodes2 ; j++)
3028           if ( (j - spp) % (( nonodes2 / 72 ) + 1 ) == 0 )
3029             putchar('-');
3030         printf("!\n");
3031 #ifdef WIN32
3032         phyFillScreenColor();
3033 #endif
3034       }
3035       succeeded = true;
3036       while (succeeded) {
3037         succeeded = false;
3038         if (global && nextsp == spp && progress) {
3039           printf("   ");
3040           fflush(stdout);
3041         }
3042         if (global && nextsp == spp)
3043           globrearrange();
3044         else
3045           rearrange(curtree.start, curtree.start->back);
3046         if (global && nextsp == spp && progress)
3047           putchar('\n');
3048       }
3049       nextsp++;
3050     }
3051     if (global && progress) {
3052       putchar('\n');
3053       fflush(stdout);
3054 #ifdef WIN32
3055       phyFillScreenColor();
3056 #endif
3057     }
3058     promlcopy(&curtree, &bestree, nonodes2, rcategs);
3059     if (njumble > 1) {
3060       if (jumb == 1)
3061         promlcopy(&bestree, &bestree2, nonodes2, rcategs);
3062       else
3063         if (bestree2.likelihood < bestree.likelihood)
3064           promlcopy(&bestree, &bestree2, nonodes2, rcategs);
3065     }
3066     if (jumb == njumble) {
3067       if (njumble > 1)
3068         promlcopy(&bestree2, &curtree, nonodes2, rcategs);
3069       curtree.start = curtree.nodep[outgrno - 1]->back;
3070       for (i = 0; i < nonodes2; i++) {
3071         if (i < spp)
3072           curtree.nodep[i]->initialized = false;
3073         else {
3074           curtree.nodep[i]->initialized = false;
3075           curtree.nodep[i]->next->initialized = false;
3076           curtree.nodep[i]->next->next->initialized = false;
3077         }
3078       }
3079       treevaluate();
3080       proml_printree();
3081       summarize();
3082       if (trout) {
3083         col = 0;
3084         dnaml_treeout(curtree.start);
3085       }
3086     }
3087   }
3088   if (usertree) {
3089     free(l0gl);
3090     for (i=0; i < shimotrees; i++)
3091       free(l0gf[i]);
3092     free(l0gf);
3093   }
3094   prot_freetable();
3095   if (jumb < njumble)
3096     return;
3097   free(contribution);
3098   free(mp);
3099   for (i=0; i < endsite; i++)
3100      free(term[i]);
3101   free(term);
3102   for (i=0; i < endsite; i++)
3103      free(slopeterm[i]);
3104   free(slopeterm);
3105   for (i=0; i < endsite; i++)
3106      free(curveterm[i]);
3107   free(curveterm);
3108   free_all_protx(nonodes2, curtree.nodep);
3109   if (!usertree) {
3110     free_all_protx(nonodes2, bestree.nodep);
3111     free_all_protx(nonodes2, priortree.nodep);
3112     if (njumble > 1)
3113       free_all_protx(nonodes2, bestree2.nodep);
3114   }
3115   if (progress) {
3116     printf("\n\nOutput written to file \"%s\"\n\n", outfilename);
3117     if (trout)
3118       printf("Tree also written onto file \"%s\"\n", outtreename);
3119     putchar('\n');
3120   }
3121 }  /* maketree */
3122
3123
3124 void clean_up()
3125 {
3126   /* Free and/or close stuff */
3127   long i;
3128
3129     free (rrate);
3130     free (probcat);
3131     free (rate);
3132     /* Seems to require freeing every time... */
3133     for (i = 0; i < spp; i++) {
3134       free (y[i]);
3135     }
3136     free (y);
3137     free (nayme);
3138     free (enterorder);
3139     free (category);
3140     free (weight);
3141     free (alias);
3142     free (ally);
3143     free (location);
3144     free (aliasweight);
3145     free (probmat);
3146     free (eigmat);
3147
3148   FClose(infile);
3149   FClose(outfile);
3150   FClose(outtree);
3151 #ifdef MAC
3152   fixmacfile(outfilename);
3153   fixmacfile(outtreename);
3154 #endif
3155 }   /* clean_up */
3156
3157
3158 int main(int argc, Char *argv[])
3159 {  /* Protein Maximum Likelihood */
3160 #ifdef MAC
3161   argc = 1;             /* macsetup("ProML","");        */
3162   argv[0] = "ProML";
3163 #endif
3164   init(argc,argv);
3165   progname = argv[0];
3166   openfile(&infile,INFILE,"input file","r",argv[0],infilename);
3167   openfile(&outfile,OUTFILE,"output file","w",argv[0],outfilename);
3168   mulsets = false;
3169   datasets = 1;
3170   firstset = true;
3171   ibmpc = IBMCRT;
3172   ansi = ANSICRT;
3173   grbg = NULL;
3174   doinit();
3175   if (ctgry)
3176     openfile(&catfile,CATFILE,"categories file","r",argv[0],catfilename);
3177   if (weights || justwts)
3178     openfile(&weightfile,WEIGHTFILE,"weights file","r",argv[0],weightfilename);
3179   if (trout)
3180     openfile(&outtree,OUTTREE,"output tree file","w",argv[0],outtreename);
3181   for (ith = 1; ith <= datasets; ith++) {
3182     if (datasets > 1) {
3183       fprintf(outfile, "Data set # %ld:\n", ith);
3184       printf("\nData set # %ld:\n", ith);
3185     }
3186     getinput();
3187     if (ith == 1)
3188       firstset = false;
3189     for (jumb = 1; jumb <= njumble; jumb++) {
3190       max_num_sibs = 0;
3191       maketree();
3192     }
3193   }
3194
3195   clean_up();
3196   printf("Done.\n\n");
3197 #ifdef WIN32
3198   phyRestoreConsoleAttributes();
3199 #endif
3200   return 0;
3201 }  /* Protein Maximum Likelihood */
3202