#include "phylip.h" #include "seq.h" /* version 3.6. (c) Copyright 1986-2004 by the University of Washington and by Joseph Felsenstein. Written by Joseph Felsenstein and Lucas Mix. Permission is granted to copy and use this program provided no fee is charged for it and provided that this copyright notice is not removed. */ #define epsilon 0.0001 /* used in makenewv, getthree, update */ #define over 60 typedef long vall[maxcategs]; typedef double contribarr[maxcategs]; #ifndef OLDC /* function prototypes */ void init_protmats(void); void getoptions(void); void makeprotfreqs(void); void allocrest(void); void doinit(void); void inputoptions(void); void input_protdata(long); void makeweights(void); void prot_makevalues(long, pointarray, long, long, sequence, steptr); void getinput(void); void prot_inittable(void); void alloc_pmatrix(long); void make_pmatrix(double **, double **, double **, long, double, double, double *, double **); void prot_nuview(node *); void getthree(node *p, double thigh, double tlow); void makenewv(node *); void update(node *); void smooth(node *); void promlk_add(node *, node *, node *, boolean); void promlk_re_move(node **, node **, boolean); double prot_evaluate(node *); void tryadd(node *, node **, node **); void addpreorder(node *, node *, node *, boolean, boolean); void restoradd(node *, node *, node *, double); void tryrearr(node *, boolean *); void repreorder(node *, boolean *); void rearrange(node **); void nodeinit(node *); void initrav(node *); void travinit(node *); void travsp(node *); void treevaluate(void); void promlk_coordinates(node *, long *); void promlk_drawline(long, double); void promlk_printree(void); void describe(node *); void prot_reconstr(node *, long); void rectrav(node *, long, long); void summarize(void); void promlk_treeout(node *); void initpromlnode(node **, node **, node *, long, long, long *, long *, initops, pointarray, pointarray, Char *, Char *, FILE *); void tymetrav(node *, double *); void free_all_protx(long, pointarray); void maketree(void); void clean_up(void); void reallocsites(void); void prot_freetable(void); void free_pmatrix(long sib); /* function prototypes */ #endif double **tbl; Char infilename[100], outfilename[100], intreename[100], outtreename[100], catfilename[100], weightfilename[100]; double *rrate; long sites, weightsum, categs, datasets, ith, njumble, jumb, numtrees, shimotrees; /* sites = number of sites in actual sequences numtrees = number of user-defined trees */ long inseed, inseed0, mx, mx0, mx1; boolean global, jumble, lngths, trout, usertree, weights, rctgry, ctgry, auto_, progress, mulsets, firstset, hypstate, smoothit, polishing, justwts, gama, invar, usejtt, usepmb, usepam; tree curtree, bestree, bestree2; node *qwhere, *grbg; double sumrates, cv, alpha, lambda, lambda1, invarfrac; long *enterorder; steptr aliasweight; double *rate; longer seed; double *probcat; contribarr *contribution; char aachar[26]="ARNDCQEGHILKMFPSTWYVBZX?*-"; char *progname; long rcategs, nonodes2; /* Local variables for maketree, propagated globally for C version: */ long k, maxwhich, col; double like, bestyet, tdelta, lnlike, slope, curv, maxlogl; boolean lastsp, smoothed, succeeded; double *l0gl; double x[3], lnl[3]; double expon1i[maxcategs], expon1v[maxcategs], expon2i[maxcategs], expon2v[maxcategs]; node *there; double **l0gf; Char ch, ch2; long **mp; /* Variables introduced to allow for protein probability calculations */ long max_num_sibs; /* maximum number of siblings used in a */ /* nuview calculation. determines size */ /* final size of pmatrices */ double *eigmat; /* eig matrix variable */ double **probmat; /* prob matrix variable */ double ****dpmatrix; /* derivative of pmatrix */ double ****ddpmatrix; /* derivative of xpmatrix */ double *****pmatrices; /* matrix of probabilities of protien */ /* conversion. The 5 subscripts refer */ /* to sibs, rcategs, categs, final and */ /* initial states, respectively. */ double freqaa[20]; /* amino acid frequencies */ /* this jtt matrix decomposition due to Elisabeth Tillier */ static double jtteigmat[] = {0.0, -0.007031123, -0.006484345, -0.006086499, -0.005514432, -0.00772664, -0.008643413, -0.010620756, -0.009965552, -0.011671808, -0.012222418,-0.004589201, -0.013103714, -0.014048038, -0.003170582, -0.00347935, -0.015311677, -0.016021194, -0.017991454, -0.018911888}; static double jttprobmat[20][20] = {{0.076999996, 0.051000003, 0.043000004, 0.051999998, 0.019999996, 0.041, 0.061999994, 0.073999997, 0.022999999, 0.052000004, 0.090999997, 0.058999988, 0.024000007, 0.04, 0.050999992, 0.069, 0.059000006, 0.014000008, 0.032000004, 0.066000005}, {0.015604455, -0.068062363, 0.020106264, 0.070723273, 0.011702977, 0.009674053, 0.074000798, -0.169750458, 0.005560808, -0.008208636, -0.012305869, -0.063730179, -0.005674643, -0.02116828, 0.104586169, 0.016480839, 0.016765139, 0.005936994, 0.006046367, -0.0082877}, {-0.049778281, -0.007118197, 0.003801272, 0.070749616, 0.047506147, 0.006447017, 0.090522425, -0.053620432, -0.008508175, 0.037170603, 0.051805545, 0.015413608, 0.019939916, -0.008431976, -0.143511376, -0.052486072, -0.032116542, -0.000860626, -0.02535993, 0.03843545}, {-0.028906423, 0.092952047, -0.009615343, -0.067870117, 0.031970392, 0.048338335, -0.054396304, -0.135916654, 0.017780083, 0.000129242, 0.031267424, 0.116333586, 0.007499746, -0.032153596, 0.033517051, -0.013719269, -0.00347293, -0.003291821, -0.02158326, -0.008862168}, {0.037181176, -0.023106564, -0.004482225, -0.029899635, 0.118139633, -0.032298569, -0.04683198, 0.05566988, -0.012622847, 0.002023096, -0.043921088, -0.04792557, -0.003452711, -0.037744513, 0.020822974, 0.036580187, 0.02331425, -0.004807711, -0.017504496, 0.01086673}, {0.044754061, -0.002503471, 0.019452517, -0.015611487, -0.02152807, -0.013131425, -0.03465365, -0.047928912, 0.020608851, 0.067843095, -0.122130014, 0.002521499, 0.013021646, -0.082891087, -0.061590119, 0.016270856, 0.051468938, 0.002079063, 0.081019713, 0.082927944}, {0.058917882, 0.007320741, 0.025278141, 0.000357541, -0.002831285, -0.032453034, -0.010177288, -0.069447924, -0.034467324, 0.011422358, -0.128478324, 0.04309667, -0.015319944, 0.113302422, -0.035052393, 0.046885372, 0.06185183, 0.00175743, -0.06224497, 0.020282093}, {-0.014562092, 0.022522921, -0.007094389, 0.03480089, -0.000326144, -0.124039037, 0.020577906, -0.005056454, -0.081841576, -0.004381786, 0.030826152, 0.091261631, 0.008878828, -0.02829487, 0.042718836, -0.011180886, -0.012719227, -0.000753926, 0.048062375, -0.009399129}, {0.033789571, -0.013512235, 0.088010984, 0.017580292, -0.006608005, -0.037836971, -0.061344686, -0.034268357, 0.018190209, -0.068484614, 0.120024744, -0.00319321, -0.001349477, -0.03000546, -0.073063759, 0.081912399, 0.0635245, 0.000197, -0.002481798, -0.09108114}, {-0.113947615, 0.019230545, 0.088819683, 0.064832765, 0.001801467, -0.063829682, -0.072001633, 0.018429333, 0.057465965, 0.043901014, -0.048050874, -0.001705918, 0.022637173, 0.017404665, 0.043877902, -0.017089594, -0.058489485, 0.000127498, -0.029357194, 0.025943972}, {0.01512923, 0.023603725, 0.006681954, 0.012360216, -0.000181447, -0.023011838, -0.008960024, -0.008533239, 0.012569835, 0.03216118, 0.061986403, -0.001919083, -0.1400832, -0.010669741, -0.003919454, -0.003707024, -0.026806029, -0.000611603, -0.001402648, 0.065312824}, {-0.036405351, 0.020816769, 0.011408213, 0.019787053, 0.038897829, 0.017641789, 0.020858533, -0.006067252, 0.028617353, -0.064259496, -0.081676567, 0.024421823, -0.028751676, 0.07095096, -0.024199434, -0.007513119, -0.028108766, -0.01198095, 0.111761119, -0.076198809}, {0.060831772, 0.144097327, -0.069151377, 0.023754576, -0.003322955, -0.071618574, 0.03353154, -0.02795295, 0.039519769, -0.023453968, -0.000630308, -0.098024591, 0.017672997, 0.003813378, -0.009266499, -0.011192111, 0.016013873, -0.002072968, -0.010022044, -0.012526904}, {-0.050776604, 0.092833081, 0.044069596, 0.050523021, -0.002628417, 0.076542572, -0.06388631, -0.00854892, -0.084725311, 0.017401063, -0.006262541, -0.094457679, -0.002818678, -0.0044122, -0.002883973, 0.028729685, -0.004961596, -0.001498627, 0.017994575, -0.000232779}, {-0.01894566, -0.007760205, -0.015160993, -0.027254587, 0.009800903, -0.013443561, -0.032896517, -0.022734138, -0.001983861, 0.00256111, 0.024823166, -0.021256768, 0.001980052, 0.028136263, -0.012364384, -0.013782446, -0.013061091, 0.111173981, 0.021702122, 0.00046654}, {-0.009444193, -0.042106824, -0.02535015, -0.055125574, 0.006369612, -0.02945416, -0.069922064, -0.067221068, -0.003004999, 0.053624311, 0.128862984, -0.057245803, 0.025550508, 0.087741073, -0.001119043, -0.012036202, -0.000913488, -0.034864475, 0.050124813, 0.055534723}, {0.145782464, -0.024348311, -0.031216873, 0.106174443, 0.00202862, 0.02653866, -0.113657267, -0.00755018, 0.000307232, -0.051241158, 0.001310685, 0.035275877, 0.013308898, 0.002957626, -0.002925034, -0.065362319, -0.071844582, 0.000475894, -0.000112419, 0.034097762}, {0.079840455, 0.018769331, 0.078685899, -0.084329807, -0.00277264, -0.010099754, 0.059700608, -0.019209715, -0.010442992, -0.042100476, -0.006020556, -0.023061786, 0.017246106, -0.001572858, -0.006703785, 0.056301316, -0.156787357, -0.000303638, 0.001498195, 0.051363455}, {0.049628261, 0.016475144, 0.094141653, -0.04444633, 0.005206131, -0.001827555, 0.02195624, 0.013066683, -0.010415582, -0.022338403, 0.007837197, -0.023397671, -0.002507095, 0.005177694, 0.017109561, -0.202340113, 0.069681441, 0.000120736, 0.002201146, 0.004670849}, {0.089153689, 0.000233354, 0.010826822, -0.004273519, 0.001440618, 0.000436077, 0.001182351, -0.002255508, -0.000700465, 0.150589876, -0.003911914, -0.00050154, -0.004564983, 0.00012701, -0.001486973, -0.018902754, -0.054748555, 0.000217377, -0.000319302, -0.162541651}}; static double pameigmat[] = {0.0, -0.002350753691875762, -0.002701991863800379, -0.002931612442853115, -0.004262492032364507, -0.005395980482561625, -0.007141172690079523, -0.007392844756151318, -0.007781761342200766, -0.00810032066366362, -0.00875299712761124, -0.01048227332164386, -0.01109594097332267, -0.01298616073142234, -0.01342036228188581, -0.01552599145527578, -0.01658762802054814, -0.0174893445623765, -0.01933280832903272, -0.02206353522613025}; static double pamprobmat[20][20] = {{0.087683339901135, 0.04051291829598762, 0.04087846315185977, 0.04771603459744777, 0.03247095396561266, 0.03784612688594957, 0.0504933695604875, 0.0898249006830755, 0.03285885059543713, 0.0357514442352119, 0.0852464099207521, 0.07910313444070642, 0.01488243946396588, 0.04100101908956829, 0.05158026947089499, 0.06975497205982451, 0.05832757042475474, 0.00931264523877807, 0.03171540880870517, 0.06303972920984541}, {0.01943453646811026, -0.004492574160484092, 0.007694891061220776, 0.01278399096887701, 0.0106157418450234, 0.007542140341575122, 0.01326994069032819, 0.02615565199894889, 0.003123125764490066, 0.002204507682495444, -0.004782898215768979, 0.01204241965177619, 0.0007847400096924341, -0.03043626073172116, 0.01221202591902536, 0.01100527004684405, 0.01116495631339549, -0.0925364931988571, -0.02622065387931562, 0.00843494142432107}, {0.01855357100209072, 0.01493642835763868, 0.0127983090766285, 0.0200533250704364, -0.1681898360107787, 0.01551657969909255, 0.02128060163107209, 0.03100633591848964, 0.00845480845269879, 0.000927149370785571, 0.00937207565817036, 0.03490557769673472, 0.00300443019551563, -0.02590837220264415, 0.01329376859943192, 0.006854110889741407, 0.01102593860528263, 0.003360844186685888, -0.03459712356647764, 0.003351477369404443}, {0.02690642688200102, 0.02131745801890152, 0.0143626616005213, 0.02405101425725929, 0.05041008641436849, 0.01430925051050233, 0.02362114036816964, 0.04688381789373886, 0.005250115453626377, -0.02040112168595516, -0.0942720776915669, 0.03773004996758644, -0.00822831940782616, -0.1164872809439224, 0.02286281877257392, 0.02849551240669926, 0.01468856796295663, 0.02377110964207936, -0.094380545436577, -0.02089068498518036}, {0.00930172577225213, 0.01493463068441099, 0.020186920775608, 0.02892154953912524, -0.01224593358361567, 0.01404228329986624, 0.02671186617119041, 0.04537535161795231, 0.02229995804098249, -0.04635704133961575, -0.1966910360247138, 0.02796648065439046, -0.02263484732621436, 0.0440490503242072, 0.01148782948302166, 0.01989170531824069, 0.001306805142981245, -0.005676690969116321, 0.07680476281625202, -0.07967537039721849}, {0.06602274245435476, -0.0966661981471856, -0.005241648783844579, 0.00859135188171146, -0.007762129660943368, -0.02888965572526196, 0.003592291525888222, 0.1668410669287673, -0.04082039290551406, 0.005233775047553415, -0.01758244726137135, -0.1493955762326898, -0.00855819137835548, 0.004211419253492328, 0.01929306335052688, 0.03008056746359405, 0.0190444422412472, 0.005577189741419315, 0.0000874156155112068, 0.02634091459108298}, {0.01933897472880726, 0.05874583569377844, -0.02293534606228405, -0.07206314017962175, -0.004580681581546643, -0.0628814337610561, -0.0850783812795136, 0.07988417636610614, -0.0852798990133397, 0.01649047166155952, -0.05416647263757423, 0.1089834536254064, 0.005093403979413865, 0.02520300254161142, 0.0005951431406455604, 0.02441251821224675, 0.02796099482240553, -0.002574933994926502, -0.007172237553012804, 0.03002455129086954}, {0.04041118479094272, -0.002476225672095412, -0.01494505811263243, -0.03759443758599911, -0.00892246902492875, -0.003634714029239211, -0.03085671837973749, -0.126176309029931, 0.005814031139083794, 0.01313561962646063, -0.04760487162503322, -0.0490563712725484, -0.005082243450421558, -0.01213634309383557, 0.1806666927079249, 0.02111663336185495, 0.02963486860587087, -0.0000175020101657785, 0.01197155383597686, 0.0357526792184636}, {-0.01184769557720525, 0.01582776076338872, -0.006570708266564639, -0.01471915653734024, 0.00894343616503608, 0.00562664968033149, -0.01465878888356943, 0.05365282692645818, 0.00893509735776116, -0.05879312944436473, 0.0806048683392995, -0.007722897986905326, -0.001819943882718859, 0.0942535573077267, 0.07483883782251654, 0.004354639673913651, -0.02828804845740341, -0.001318222184691827, -0.07613149604246563, -0.1251675867732172}, {0.00834167031558193, -0.01509357596974962, 0.007098172811092488, 0.03127677418040319, 0.001992448468465455, 0.00915441566808454, 0.03430175973499201, -0.0730648147535803, -0.001402707145575659, 0.04780949194330815, -0.1115035603461273, -0.01292297197609604, -0.005056270550868528, 0.1112053349612027, -0.03801929822379964, -0.001191241001736563, 0.01872874622910247, 0.0005314214903865993, -0.0882576318311789, 0.07607183599610171}, {-0.01539460099727769, 0.04988596184297883, -0.01187240760647617, -0.06987843637091853, -0.002490472846497859, 0.01009857892494956, -0.07473588067847209, 0.0906009925879084, 0.1243612446505172, 0.02152806401345371, -0.03504879644860233, -0.06680752427613573, -0.005574485153629651, 0.001518282948127752, -0.01999168507510701, -0.01478606199529457, -0.02203749419458996, -0.00132680708294333, -0.01137505997867614, 0.05332658773667142}, {-0.06104378736432388, 0.0869446603393548, -0.03298331234537257, 0.03128515657456024, 0.003906358569208259, 0.03578694104193928, 0.06241936133189683, 0.06182827284921748, -0.05566564263245907, 0.02640868588189002, -0.01349751243059039, -0.05507866642582638, -0.006671347738489326, -0.001470096466016046, 0.05185743641479938, -0.07494697511168257, -0.1175185439057584, -0.001188074094105709, 0.00937934805737347, 0.05024773745437657}, {-0.07252555582124737, -0.116554459356382, 0.003605361887406413, -0.00836518656029184, 0.004615715410745561, 0.005105376617651312, -0.00944938657024391, 0.05602449420950007, 0.02722719610561933, 0.01959357494748446, -0.0258655103753962, 0.1440733975689835, 0.01446782819722976, 0.003718896062070054, 0.05825843045655135, -0.06230154142733073, -0.07833704962300169, 0.003160836143568724, -0.001169873777936648, 0.03471745590503304}, {-0.03204352258752698, 0.01019272923862322, 0.04509668708733181, 0.05756522429120813, -0.0004601149081726732, -0.0984718150777423, -0.01107826100664925, -0.005680277810520585, 0.01962359392320817, 0.01550006899131986, 0.05143956925922197, 0.02462476682588468, -0.0888843861002653, -0.00171553583659411, 0.01606331750661664, 0.001176847743518958, -0.02070972978912828, -0.000341523293579971, -0.002654732745607882, 0.02075709428885848}, {0.03595199666430258, -0.02800219615234468, -0.04341570015493925, -0.0748275906176658, 0.0001051403676377422, 0.1137431321746627, 0.005852087565974318, 0.003443037513847801, -0.02481931657706633, -0.003651181839831423, 0.03195794176786321, 0.04135411406392523, -0.07562030263210619, 0.001769332364699, -0.01984381173403915, -0.005029750745010152, 0.02649253902476472, 0.000518085571702734, 0.001062936684474851, 0.01295950668914449}, {-0.16164552322896, -0.0006050035060464324, 0.0258380054414968, 0.003188424740960557, -0.0002058911341821877, 0.03157555987384681, -0.01678913462596107, 0.03096216145389774, -0.0133791110666919, 0.1125249625204277, -0.00769017706442472, -0.02653938062180483, -0.002555329863523985, -0.00861833362947954, 0.01775148884754278, 0.02529310679774722, 0.0826243417011238, -0.0001036728183032624, 0.001963562313294209, -0.0935900561309786}, {0.1652394174588469, -0.002814245280784351, -0.0328982001821263, -0.02000104712964131, 0.0002208121995725443, -0.02733462178511839, 0.02648078162927627, -0.01788316626401427, 0.01630747623755998, 0.1053849023838147, -0.005447706553811218, 0.01810876922536839, -0.001808914710282444, -0.007687912115607397, -0.01332593672114388, -0.02110750894891371, -0.07456116592983384, 0.000219072589592394, 0.001270886972191055, -0.1083616930749109}, {0.02453279389716254, -0.005820072356487439, 0.100260287284095, 0.01277522280305745, -0.003184943445296999, 0.05814689527984152, -0.0934012278200201, -0.03017986487349484, -0.03136625380994165, 0.00988668352785117, -0.00358900410973142, -0.02017443675004764, 0.000915384582922184, -0.001460963415183106, -0.01370112443251124, 0.1130040979284457, -0.1196161771323699, -0.0005800211204222045, -0.0006153403201024954, 0.00416806428223025}, {-0.0778089244252535, -0.007055161182430869, -0.0349307504860869, -0.0811915584276571, -0.004689825871599125, -0.03726108871471753, 0.1072225647141469, -0.00917015113070944, 0.01381628985996913, -0.00123227881492089, 0.001815954515275675, 0.005708744099349901, -0.0001448985044877925, -0.001306578795561384, -0.006992743514185243, 0.1744720240732789, -0.05353628497814023, -0.0007613684227234787, -0.0003550282315997644, 0.01340106423804634}, {-0.0159527329868513, -0.007622151568160798, -0.1389875105184963, 0.1165051999914764, -0.002217810389087748, 0.01550003226513692, -0.07427664222230566, -0.003371438498619264, 0.01385754771325365, 0.004759020167383304, 0.001624078805220564, 0.02011638303109029, -0.001717827082842178, -0.0007424036708598594, -0.003978884451898934, 0.0866418927301209, -0.01280817739158123, -0.00023039242454603, 0.002309205802479111, 0.0005926106991001195}}; /* this pmb matrix decomposition due to Elisabeth Tillier */ static double pmbeigmat[20] = {0.0000001586972220,-1.8416770496147100, -1.6025046986139100,-1.5801012515121300, -1.4987794099715900,-1.3520794233801900,-1.3003469390479700,-1.2439503327631300, -1.1962574080244200,-1.1383730501367500,-1.1153278910708000,-0.4934843510654760, -0.5419014550215590,-0.9657997830826700,-0.6276075673757390,-0.6675927795018510, -0.6932641383465870,-0.8897872681859630,-0.8382698977371710,-0.8074694642446040}; static double pmbprobmat[20][20] = {{0.0771762457248147,0.0531913844998640,0.0393445076407294,0.0466756566755510, 0.0286348361997465,0.0312327748383639,0.0505410248721427,0.0767106611472993, 0.0258916271688597,0.0673140562194124,0.0965705469252199,0.0515979465932174, 0.0250628079438675,0.0503492018628350,0.0399908189418273,0.0641898881894471, 0.0517539616710987,0.0143507440546115,0.0357994592438322,0.0736218495862984}, {0.0368263046116572,-0.0006728917107827,0.0008590805287740,-0.0002764255356960, 0.0020152937187455,0.0055743720652960,0.0003213317669367,0.0000449190281568, -0.0004226254397134,0.1805040629634510,-0.0272246813586204,0.0005904606533477, -0.0183743200073889,-0.0009194625608688,0.0008173657533167,-0.0262629806302238, 0.0265738757209787,0.0002176606241904,0.0021315644838566,-0.1823229927207580}, {-0.0194800075560895,0.0012068088610652,-0.0008803318319596,-0.0016044273960017, -0.0002938633803197,-0.0535796754602196,0.0155163896648621,-0.0015006360762140, 0.0021601372013703,0.0268513218744797,-0.1085292493742730,0.0149753083138452, 0.1346457366717310,-0.0009371698759829,0.0013501708044116,0.0346352293103622, -0.0276963770242276,0.0003643142783940,0.0002074817333067,-0.0174108903914110}, {0.0557839400850153,0.0023271577185437,0.0183481103396687,0.0023339480096311, 0.0002013267015151,-0.0227406863569852,0.0098644845475047,0.0064721276774396, 0.0001389408104210,-0.0473713878768274,-0.0086984445005797,0.0026913674934634, 0.0283724052562196,0.0001063665179457,0.0027442574779383,-0.1875312134708470, 0.1279864877057640,0.0005103347834563,0.0003155113168637,0.0081451082759554}, {0.0037510125027265,0.0107095920636885,0.0147305410328404,-0.0112351252180332, -0.0001500408626446,-0.1523450933729730,0.0611532413339872,-0.0005496748939503, 0.0048714378736644,-0.0003826320053999,0.0552010244407311,0.0482555671001955, -0.0461664995115847,-0.0021165008617978,-0.0004574454232187,0.0233755883688949, -0.0035484915422384,0.0009090698422851,0.0013840637687758,-0.0073895139302231}, {-0.0111512564930024,0.1025460064723080,0.0396772456883791,-0.0298408501361294, -0.0001656742634733,-0.0079876311843289,0.0712644184507945,-0.0010780604625230, -0.0035880882043592,0.0021070399334252,0.0016716329894279,-0.1810123023850110, 0.0015141703608724,-0.0032700852781804,0.0035503782441679,0.0118634302028026, 0.0044561606458028,-0.0001576678495964,0.0023470722225751,-0.0027457045397157}, {0.1474525743949170,-0.0054432538500293,0.0853848892349828,-0.0137787746207348, -0.0008274830358513,0.0042248844582553,0.0019556229305563,-0.0164191435175148, -0.0024501858854849,0.0120908948084233,-0.0381456105972653,0.0101271614855119, -0.0061945941321859,0.0178841099895867,-0.0014577779202600,-0.0752120602555032, -0.1426985695849920,0.0002862275078983,-0.0081191734261838,0.0313401149422531}, {0.0542034611735289,-0.0078763926211829,0.0060433542506096,0.0033396210615510, 0.0013965072374079,0.0067798903832256,-0.0135291136622509,-0.0089982442731848, -0.0056744537593887,-0.0766524225176246,0.1881210263933930,-0.0065875518675173, 0.0416627569300375,-0.0953804133524747,-0.0012559228448735,0.0101622644292547, -0.0304742453119050,0.0011702318499737,0.0454733434783982,-0.1119239362388150}, {0.1069409037912470,0.0805064400880297,-0.1127352030714600,0.1001181253523260, -0.0021480427488769,-0.0332884841459003,-0.0679837575848452,-0.0043812841356657, 0.0153418716846395,-0.0079441315103188,-0.0121766182046363,-0.0381127991037620, -0.0036338726532673,0.0195324059593791,-0.0020165963699984,-0.0061222685010268, -0.0253761448771437,-0.0005246410999057,-0.0112205170502433,0.0052248485517237}, {-0.0325247648326262,0.0238753651653669,0.0203684886605797,0.0295666232678825, -0.0003946714764213,-0.0157242718469554,-0.0511737848084862,0.0084725632040180, -0.0167068828528921,0.0686962159427527,-0.0659702890616198,-0.0014289912494271, -0.0167000964093416,-0.1276689083678200,0.0036575057830967,-0.0205958145531018, 0.0000368919612829,0.0014413626622426,0.1064360941926030,0.0863372661517408}, {-0.0463777468104402,0.0394712148670596,0.1118686750747160,0.0440711686389031, -0.0026076286506751,-0.0268454015202516,-0.1464943067133240,-0.0137514051835380, -0.0094395514284145,-0.0144124844774228,0.0249103379323744,-0.0071832157138676, 0.0035592787728526,0.0415627419826693,0.0027040097365669,0.0337523666612066, 0.0316121324137152,-0.0011350177559026,-0.0349998884574440,-0.0302651879823361}, {0.0142360925194728,0.0413145623127025,0.0324976427846929,0.0580930922002398, -0.0586974207121084,0.0202001168873069,0.0492204086749069,0.1126593173463060, 0.0116620013776662,-0.0780333711712066,-0.1109786767320410,0.0407775100936731, -0.0205013161312652,-0.0653458585025237,0.0347351829703865,0.0304448983224773, 0.0068813748197884,-0.0189002309261882,-0.0334507528405279,-0.0668143558699485}, {-0.0131548829657936,0.0044244322828034,-0.0050639951827271,-0.0038668197633889, -0.1536822386530220,0.0026336969165336,0.0021585651200470,-0.0459233839062969, 0.0046854727140565,0.0393815434593599,0.0619554007991097,0.0027456299925622, 0.0117574347936383,0.0373018612990383,0.0024818527553328,-0.0133956606027299, -0.0020457128424105,0.0154178819990401,0.0246524142683911,0.0275363065682921}, {-0.1542307272455030,0.0364861558267547,-0.0090880407008181,0.0531673937889863, 0.0157585615170580,0.0029986538457297,0.0180194047699875,0.0652152443589317, 0.0266842840376180,0.0388457366405908,0.0856237634510719,0.0126955778952183, 0.0099593861698250,-0.0013941794862563,0.0294065511237513,-0.1151906949298290, -0.0852991447389655,0.0028699120202636,-0.0332087026659522,0.0006811857297899}, {0.0281300736924501,-0.0584072081898638,-0.0178386569847853,-0.0536470338171487, -0.0186881656029960,-0.0240008730656106,-0.0541064820498883,0.2217137098936020, -0.0260500001542033,0.0234505236798375,0.0311127151218573,-0.0494139126682672, 0.0057093465049849,0.0124937286655911,-0.0298322975915689,0.0006520211333102, -0.0061018680727128,-0.0007081999479528,-0.0060523759094034,0.0215845995364623}, {0.0295321046399105,-0.0088296411830544,-0.0065057049917325,-0.0053478115612781, -0.0100646496794634,-0.0015473619084872,0.0008539960632865,-0.0376381933046211, -0.0328135588935604,0.0672161874239480,0.0667626853916552,-0.0026511651464901, 0.0140451514222062,-0.0544836996133137,0.0427485157912094,0.0097455780205802, 0.0177309072915667,-0.0828759701187452,-0.0729504795471370,0.0670731961252313}, {0.0082646581043963,-0.0319918630534466,-0.0188454445200422,-0.0374976353856606, 0.0037131290686848,-0.0132507796987883,-0.0306958830735725,-0.0044119395527308, -0.0140786756619672,-0.0180512599925078,-0.0208243802903953,-0.0232202769398931, -0.0063135878270273,0.0110442171178168,0.1824538048228460,-0.0006644614422758, -0.0069909097436659,0.0255407650654681,0.0099119399501151,-0.0140911517070698}, {0.0261344441524861,-0.0714454044548650,0.0159436926233439,0.0028462736216688, -0.0044572637889080,-0.0089474834434532,-0.0177570282144517,-0.0153693244094452, 0.1160919467206400,0.0304911481385036,0.0047047513411774,-0.0456535116423972, 0.0004491494948617,-0.0767108879444462,-0.0012688533741441,0.0192445965934123, 0.0202321954782039,0.0281039933233607,-0.0590403018490048,0.0364080426546883}, {0.0115826306265004,0.1340228176509380,-0.0236200652949049,-0.1284484655137340, -0.0004742338006503,0.0127617346949511,-0.0428560878860394,0.0060030732454125, 0.0089182609926781,0.0085353834972860,0.0048464809638033,0.0709740071429510, 0.0029940462557054,-0.0483434904493132,-0.0071713680727884,-0.0036840391887209, 0.0031454003250096,0.0246243550241551,-0.0449551277644180,0.0111449232769393}, {0.0140356721886765,-0.0196518236826680,0.0030517022326582,0.0582672093364850, -0.0000973895685457,0.0021704767224292,0.0341806268602705,-0.0152035987563018, -0.0903198657739177,0.0259623214586925,0.0155832497882743,-0.0040543568451651, 0.0036477631918247,-0.0532892744763217,-0.0142569373662724,0.0104500681408622, 0.0103483945857315,0.0679534422398752,-0.0768068882938636,0.0280289727046158}} ; void init_protmats() { long l, m; eigmat = (double *) Malloc (20 * sizeof(double)); for (l = 0; l <= 19; l++) if (usejtt) eigmat[l] = jtteigmat[l]*100.0; else { if (usepmb) eigmat[l] = pmbeigmat[l]; else eigmat[l] = pameigmat[l]*100.0; } probmat = (double **) Malloc (20 * sizeof(double *)); for (l = 0; l < 20; l++) probmat[l] = (double *) Malloc (20 * sizeof(double)); for (l = 0; l <= 19; l++) for (m= 0; m <= 19; m++) if (usejtt) probmat[l][m] = jttprobmat[l][m]; else { if (usepmb) probmat[l][m] = pmbprobmat[l][m]; else probmat[l][m] = pamprobmat[l][m]; } } /* init_protmats */ void getoptions() { /* interactively set options */ long i, loopcount, loopcount2; Char ch; boolean done; boolean didchangecat, didchangercat; double probsum; fprintf(outfile, "\nAmino acid sequence\n"); fprintf(outfile, " Maximum Likelihood method with molecular "); fprintf(outfile, "clock, version %s\n\n", VERSION); putchar('\n'); auto_ = false; ctgry = false; didchangecat = false; rctgry = false; didchangercat = false; categs = 1; rcategs = 1; gama = false; global = false; hypstate = false; invar = false; jumble = false; njumble = 1; lambda = 1.0; lambda1 = 0.0; lngths = false; trout = true; usepam = false; usepmb = false; usejtt = true; usertree = false; weights = false; printdata = false; progress = true; treeprint = true; interleaved = true; loopcount = 0; do { cleerhome(); printf("\nAmino acid sequence\n"); printf(" Maximum Likelihood method with molecular clock, version %s\n\n", VERSION); printf("Settings for this run:\n"); printf(" U Search for best tree?"); if (usertree) printf(" No, use user trees in input file\n"); else printf(" Yes\n"); printf(" P JTT, PMB or PAM probability model? %s\n", usejtt ? "Jones-Taylor-Thornton" : usepmb ? "Henikoff/Tillier PMB" : "Dayhoff PAM"); if (usertree) { printf(" L Use lengths from user tree?"); if (lngths) printf(" Yes\n"); else printf(" No\n"); } printf(" C One category of substitution rates?"); if (!ctgry) printf(" Yes\n"); else printf(" %ld categories\n", categs); printf(" R Rate variation among sites?"); if (!rctgry) printf(" constant rate of change\n"); else { if (gama) printf(" Gamma distributed rates\n"); else { if (invar) printf(" Gamma+Invariant sites\n"); else printf(" user-defined HMM of rates\n"); } printf(" A Rates at adjacent sites correlated?"); if (!auto_) printf(" No, they are independent\n"); else printf(" Yes, mean block length =%6.1f\n", 1.0 / lambda); } if (!usertree) { printf(" G Global rearrangements?"); if (global) printf(" Yes\n"); else printf(" No\n"); } printf(" W Sites weighted? %s\n", (weights ? "Yes" : "No")); if (!usertree) { printf(" J Randomize input order of sequences?"); if (jumble) printf(" Yes (seed = %8ld, %3ld times)\n", inseed0, njumble); else printf(" No. Use input order\n"); } printf(" M Analyze multiple data sets?"); if (mulsets) printf(" Yes, %2ld %s\n", datasets, (justwts ? "sets of weights" : "data sets")); else printf(" No\n"); printf(" I Input sequences interleaved?"); if (interleaved) printf(" Yes\n"); else printf(" No, sequential\n"); printf(" 0 Terminal type (IBM PC, ANSI, none)?"); if (ibmpc) printf(" IBM PC\n"); if (ansi) printf(" ANSI\n"); if (!(ibmpc || ansi)) printf(" (none)\n"); printf(" 1 Print out the data at start of run"); if (printdata) printf(" Yes\n"); else printf(" No\n"); printf(" 2 Print indications of progress of run"); if (progress) printf(" Yes\n"); else printf(" No\n"); printf(" 3 Print out tree"); if (treeprint) printf(" Yes\n"); else printf(" No\n"); printf(" 4 Write out trees onto tree file?"); if (trout) printf(" Yes\n"); else printf(" No\n"); printf(" 5 Reconstruct hypothetical sequences? %s\n", (hypstate ? "Yes" : "No")); printf("\nAre these settings correct? "); printf("(type Y or the letter for one to change)\n"); #ifdef WIN32 phyFillScreenColor(); #endif scanf("%c%*[^\n]", &ch); getchar(); if (ch == '\n') ch = ' '; uppercase(&ch); done = (ch == 'Y'); if (!done) { uppercase(&ch); if (strchr("UPCRJAFWGLTMI012345", ch) != NULL){ switch (ch) { case 'C': ctgry = !ctgry; if (ctgry) { printf("\nSitewise user-assigned categories:\n\n"); initcatn(&categs); if (rate){ free(rate); } rate = (double *) Malloc( categs * sizeof(double)); didchangecat = true; initcategs(categs, rate); } break; case 'P': if (usejtt) { usejtt = false; usepmb = true; } else { if (usepmb) { usepmb = false; usepam = true; } else { usepam = false; usejtt = true; } } break; case 'R': if (!rctgry) { rctgry = true; gama = true; } else { if (gama) { gama = false; invar = true; } else { if (invar) invar = false; else rctgry = false; } } break; case 'A': auto_ = !auto_; if (auto_) { initlambda(&lambda); lambda1 = 1.0 - lambda; } break; case 'G': global = !global; break; case 'W': weights = !weights; break; case 'J': jumble = !jumble; if (jumble) initjumble(&inseed, &inseed0, seed, &njumble); else njumble = 1; break; case 'L': lngths = !lngths; break; case 'U': usertree = !usertree; break; case 'M': mulsets = !mulsets; if (mulsets) { printf("Multiple data sets or multiple weights?"); loopcount2 = 0; do { printf(" (type D or W)\n"); #ifdef WIN32 phyFillScreenColor(); #endif scanf("%c%*[^\n]", &ch2); getchar(); if (ch2 == '\n') ch2 = ' '; uppercase(&ch2); countup(&loopcount2, 10); } while ((ch2 != 'W') && (ch2 != 'D')); justwts = (ch2 == 'W'); if (justwts) justweights(&datasets); else initdatasets(&datasets); if (!jumble) { jumble = true; initjumble(&inseed, &inseed0, seed, &njumble); } } break; case 'I': interleaved = !interleaved; break; case '0': initterminal(&ibmpc, &ansi); break; case '1': printdata = !printdata; break; case '2': progress = !progress; break; case '3': treeprint = !treeprint; break; case '4': trout = !trout; break; case '5': hypstate = !hypstate; break; } } else printf("Not a possible option!\n"); } countup(&loopcount, 100); } while (!done); if (gama || invar) { loopcount = 0; do { printf( "\nCoefficient of variation of substitution rate among sites (must be positive)\n"); printf( " In gamma distribution parameters, this is 1/(square root of alpha)\n"); #ifdef WIN32 phyFillScreenColor(); #endif scanf("%lf%*[^\n]", &cv); getchar(); countup(&loopcount, 10); } while (cv <= 0.0); alpha = 1.0 / (cv * cv); } if (!rctgry) auto_ = false; if (rctgry) { printf("\nRates in HMM"); if (invar) printf(" (including one for invariant sites)"); printf(":\n"); initcatn(&rcategs); if (probcat){ free(probcat); free(rrate); } probcat = (double *) Malloc(rcategs * sizeof(double)); rrate = (double *) Malloc(rcategs * sizeof(double)); didchangercat = true; if (gama) initgammacat(rcategs, alpha, rrate, probcat); else { if (invar) { loopcount = 0; do { printf("Fraction of invariant sites?\n"); scanf("%lf%*[^\n]", &invarfrac); getchar(); countup (&loopcount, 10); } while ((invarfrac <= 0.0) || (invarfrac >= 1.0)); initgammacat(rcategs-1, alpha, rrate, probcat); for (i = 0; i < rcategs-1; i++) probcat[i] = probcat[i]*(1.0-invarfrac); probcat[rcategs-1] = invarfrac; rrate[rcategs-1] = 0.0; } else { initcategs(rcategs, rrate); initprobcat(rcategs, &probsum, probcat); } } } if (!didchangercat){ rrate = Malloc( rcategs*sizeof(double)); probcat = Malloc( rcategs*sizeof(double)); rrate[0] = 1.0; probcat[0] = 1.0; } if (!didchangecat){ rate = Malloc( categs*sizeof(double)); rate[0] = 1.0; } init_protmats(); } /* getoptions */ void makeprotfreqs() { /* calculate amino acid frequencies based on eigmat */ long i, mineig; mineig = 0; for (i = 0; i <= 19; i++) if (fabs(eigmat[i]) < fabs(eigmat[mineig])) mineig = i; memcpy(freqaa, probmat[mineig], 20 * sizeof(double)); for (i = 0; i <= 19; i++) freqaa[i] = fabs(freqaa[i]); } /* makeprotfreqs */ void reallocsites() { long i; for (i = 0; i < spp; i++) y[i] = (char *)Malloc(sites * sizeof(char)); enterorder = (long *)Malloc(spp*sizeof(long)); weight = (long *)Malloc(sites*sizeof(long)); category = (long *)Malloc(sites*sizeof(long)); alias = (long *)Malloc(sites*sizeof(long)); aliasweight = (long *)Malloc(sites*sizeof(long)); ally = (long *)Malloc(sites*sizeof(long)); location = (long *)Malloc(sites*sizeof(long)); for (i = 0; i < sites; i++) category[i] = 1; for (i = 0; i < sites; i++) weight[i] = 1; makeweights(); } void allocrest() { long i; y = (Char **)Malloc(spp*sizeof(Char *)); nayme = (naym *)Malloc(spp*sizeof(naym)); for (i = 0; i < spp; i++) y[i] = (char *)Malloc(sites * sizeof(char)); enterorder = (long *)Malloc(spp*sizeof(long)); weight = (long *)Malloc(sites*sizeof(long)); category = (long *)Malloc(sites*sizeof(long)); alias = (long *)Malloc(sites*sizeof(long)); aliasweight = (long *)Malloc(sites*sizeof(long)); ally = (long *)Malloc(sites*sizeof(long)); location = (long *)Malloc(sites*sizeof(long)); } /* allocrest */ void doinit() { /* initializes variables */ inputnumbers(&spp, &sites, &nonodes, 1); getoptions(); makeprotfreqs(); if (printdata) fprintf(outfile, "%2ld species, %3ld sites\n", spp, sites); alloctree(&curtree.nodep, nonodes, usertree); allocrest(); if (usertree) return; alloctree(&bestree.nodep, nonodes, 0); if (njumble <= 1) return; alloctree(&bestree2.nodep, nonodes, 0); } /* doinit */ void inputoptions() { long i; if (!firstset) { samenumsp(&sites, ith); reallocsites(); } if (firstset) { for (i = 0; i < sites; i++) category[i] = 1; for (i = 0; i < sites; i++) weight[i] = 1; } if (justwts || weights) inputweights(sites, weight, &weights); weightsum = 0; for (i = 0; i < sites; i++) weightsum += weight[i]; if ((ctgry && categs > 1) && (firstset || !justwts)) { inputcategs(0, sites, category, categs, "ProMLK"); if (printdata) printcategs(outfile, sites, category, "Site categories"); } if (weights && printdata) printweights(outfile, 0, sites, weight, "Sites"); fprintf(outfile, "%s model of amino acid change\n\n", (usejtt ? "Jones-Taylor-Thornton" : usepmb ? "Henikoff/Tillier PMB" : "Dayhoff PAM")); } /* inputoptions */ void input_protdata(long chars) { /* input the names and sequences for each species */ /* used by proml */ long i, j, k, l, basesread, basesnew; Char charstate; boolean allread, done; if (printdata) headings(chars, "Sequences", "---------"); basesread = 0; basesnew = 0; allread = false; while (!(allread)) { /* eat white space -- if the separator line has spaces on it*/ do { charstate = gettc(infile); } while (charstate == ' ' || charstate == '\t'); ungetc(charstate, infile); if (eoln(infile)) scan_eoln(infile); i = 1; while (i <= spp) { if ((interleaved && basesread == 0) || !interleaved) initname(i - 1); j = (interleaved) ? basesread : 0; done = false; while (!done && !eoff(infile)) { if (interleaved) done = true; while (j < chars && !(eoln(infile) || eoff(infile))) { charstate = gettc(infile); if (charstate == '\n' || charstate == '\t') charstate = ' '; if (charstate == ' ' || (charstate >= '0' && charstate <= '9')) continue; uppercase(&charstate); if ((strchr("ABCDEFGHIKLMNPQRSTVWXYZ*?-", charstate)) == NULL){ printf("ERROR: bad amino acid: %c at position %ld of species %ld\n", charstate, j, i); if (charstate == '.') { printf(" Periods (.) may not be used as gap characters.\n"); printf(" The correct gap character is (-)\n"); } exxit(-1); } j++; y[i - 1][j - 1] = charstate; } if (interleaved) continue; if (j < chars) scan_eoln(infile); else if (j == chars) done = true; } if (interleaved && i == 1) basesnew = j; scan_eoln(infile); if ((interleaved && j != basesnew) || (!interleaved && j != chars)) { printf("ERROR: SEQUENCES OUT OF ALIGNMENT AT POSITION %ld.\n", j); exxit(-1); } i++; } if (interleaved) { basesread = basesnew; allread = (basesread == chars); } else allread = (i > spp); } if (!printdata) return; for (i = 1; i <= ((chars - 1) / 60 + 1); i++) { for (j = 1; j <= spp; j++) { for (k = 0; k < nmlngth; k++) putc(nayme[j - 1][k], outfile); fprintf(outfile, " "); l = i * 60; if (l > chars) l = chars; for (k = (i - 1) * 60 + 1; k <= l; k++) { if (j > 1 && y[j - 1][k - 1] == y[0][k - 1]) charstate = '.'; else charstate = y[j - 1][k - 1]; putc(charstate, outfile); if (k % 10 == 0 && k % 60 != 0) putc(' ', outfile); } putc('\n', outfile); } putc('\n', outfile); } putc('\n', outfile); } /* input_protdata */ void makeweights() { /* make up weights vector to avoid duplicate computations */ long i; for (i = 1; i <= sites; i++) { alias[i - 1] = i; ally[i - 1] = 0; aliasweight[i - 1] = weight[i - 1]; location[i - 1] = 0; } sitesort2(sites, aliasweight); sitecombine2(sites, aliasweight); sitescrunch2(sites, 1, 2, aliasweight); for (i = 1; i <= sites; i++) { if (aliasweight[i - 1] > 0) endsite = i; } for (i = 1; i <= endsite; i++) { ally[alias[i - 1] - 1] = alias[i - 1]; location[alias[i - 1] - 1] = i; } contribution = (contribarr *) Malloc( endsite*sizeof(contribarr)); } /* makeweights */ void prot_makevalues(long categs, pointarray treenode, long endsite, long spp, sequence y, steptr alias) { /* set up fractional likelihoods at tips */ /* a version of makevalues2 found in seq.c */ /* used by proml */ long i, j, k, l; long b; for (k = 0; k < endsite; k++) { j = alias[k]; for (i = 0; i < spp; i++) { for (l = 0; l < categs; l++) { memset(treenode[i]->protx[k][l], 0, sizeof(double)*20); switch (y[i][j - 1]) { case 'A': treenode[i]->protx[k][l][0] = 1.0; break; case 'R': treenode[i]->protx[k][l][(long)arginine - (long)alanine] = 1.0; break; case 'N': treenode[i]->protx[k][l][(long)asparagine - (long)alanine] = 1.0; break; case 'D': treenode[i]->protx[k][l][(long)aspartic - (long)alanine] = 1.0; break; case 'C': treenode[i]->protx[k][l][(long)cysteine - (long)alanine] = 1.0; break; case 'Q': treenode[i]->protx[k][l][(long)glutamine - (long)alanine] = 1.0; break; case 'E': treenode[i]->protx[k][l][(long)glutamic - (long)alanine] = 1.0; break; case 'G': treenode[i]->protx[k][l][(long)glycine - (long)alanine] = 1.0; break; case 'H': treenode[i]->protx[k][l][(long)histidine - (long)alanine] = 1.0; break; case 'I': treenode[i]->protx[k][l][(long)isoleucine - (long)alanine] = 1.0; break; case 'L': treenode[i]->protx[k][l][(long)leucine - (long)alanine] = 1.0; break; case 'K': treenode[i]->protx[k][l][(long)lysine - (long)alanine] = 1.0; break; case 'M': treenode[i]->protx[k][l][(long)methionine - (long)alanine] = 1.0; break; case 'F': treenode[i]->protx[k][l][(long)phenylalanine - (long)alanine] = 1.0; break; case 'P': treenode[i]->protx[k][l][(long)proline - (long)alanine] = 1.0; break; case 'S': treenode[i]->protx[k][l][(long)serine - (long)alanine] = 1.0; break; case 'T': treenode[i]->protx[k][l][(long)threonine - (long)alanine] = 1.0; break; case 'W': treenode[i]->protx[k][l][(long)tryptophan - (long)alanine] = 1.0; break; case 'Y': treenode[i]->protx[k][l][(long)tyrosine - (long)alanine] = 1.0; break; case 'V': treenode[i]->protx[k][l][(long)valine - (long)alanine] = 1.0; break; case 'B': treenode[i]->protx[k][l][(long)asparagine - (long)alanine] = 1.0; treenode[i]->protx[k][l][(long)aspartic - (long)alanine] = 1.0; break; case 'Z': treenode[i]->protx[k][l][(long)glutamine - (long)alanine] = 1.0; treenode[i]->protx[k][l][(long)glutamic - (long)alanine] = 1.0; break; case 'X': /* unknown aa */ for (b = 0; b <= 19; b++) treenode[i]->protx[k][l][b] = 1.0; break; case '?': /* unknown aa */ for (b = 0; b <= 19; b++) treenode[i]->protx[k][l][b] = 1.0; break; case '*': /* stop codon symbol */ for (b = 0; b <= 19; b++) treenode[i]->protx[k][l][b] = 1.0; break; case '-': /* deletion event-absent data or aa */ for (b = 0; b <= 19; b++) treenode[i]->protx[k][l][b] = 1.0; break; } } } } } /* prot_makevalues */ void getinput() { long grcategs; /* reads the input data */ if (!justwts || firstset) inputoptions(); if (!justwts || firstset) input_protdata(sites); makeweights(); setuptree2(curtree); if (!usertree) { setuptree2(bestree); if (njumble > 1) setuptree2(bestree2); } grcategs = (categs > rcategs) ? categs : rcategs; prot_allocx(nonodes, grcategs, curtree.nodep, usertree); if (!usertree) { prot_allocx(nonodes, grcategs, bestree.nodep, 0); if (njumble > 1) prot_allocx(nonodes, grcategs, bestree2.nodep, 0); } prot_makevalues(rcategs, curtree.nodep, endsite, spp, y, alias); } /* getinput */ void prot_freetable(void) { long i,j,k,l; for (j = 0; j < rcategs; j++) { for (k = 0; k < categs; k++) { for (l = 0; l < 20; l++) free(ddpmatrix[j][k][l]); free(ddpmatrix[j][k]); } free(ddpmatrix[j]); } free(ddpmatrix); for (j = 0; j < rcategs; j++) { for (k = 0; k < categs; k++) { for (l = 0; l < 20; l++) free(dpmatrix[j][k][l]); free(dpmatrix[j][k]); } free(dpmatrix[j]); } free(dpmatrix); for (j = 0; j < rcategs; j++) free(tbl[j]); free(tbl); for ( i = 0 ; i < max_num_sibs ; i++ ) free_pmatrix(i); free(pmatrices); } void prot_inittable() { /* Define a lookup table. Precompute values and print them out in tables */ /* Allocate memory for the pmatrices, dpmatices and ddpmatrices */ long i, j, k, l; double sumrates; /* Allocate memory for pmatrices, the array of pointers to pmatrices */ pmatrices = (double *****) Malloc (spp * sizeof(double ****)); /* Allocate memory for the first 2 pmatrices, the matrix of conversion */ /* probabilities, but only once per run (aka not on the second jumble. */ alloc_pmatrix(0); alloc_pmatrix(1); /* Allocate memory for one dpmatrix, the first derivative matrix */ dpmatrix = (double ****) Malloc( rcategs * sizeof(double ***)); for (j = 0; j < rcategs; j++) { dpmatrix[j] = (double ***) Malloc( categs * sizeof(double **)); for (k = 0; k < categs; k++) { dpmatrix[j][k] = (double **) Malloc( 20 * sizeof(double *)); for (l = 0; l < 20; l++) dpmatrix[j][k][l] = (double *) Malloc( 20 * sizeof(double)); } } /* Allocate memory for one ddpmatrix, the second derivative matrix */ ddpmatrix = (double ****) Malloc( rcategs * sizeof(double ***)); for (j = 0; j < rcategs; j++) { ddpmatrix[j] = (double ***) Malloc( categs * sizeof(double **)); for (k = 0; k < categs; k++) { ddpmatrix[j][k] = (double **) Malloc( 20 * sizeof(double *)); for (l = 0; l < 20; l++) ddpmatrix[j][k][l] = (double *) Malloc( 20 * sizeof(double)); } } /* Allocate memory and assign values to tbl, the matrix of possible rates*/ tbl = (double **) Malloc( rcategs * sizeof(double *)); for (j = 0; j < rcategs; j++) tbl[j] = (double *) Malloc( categs * sizeof(double)); for (j = 0; j < rcategs; j++) for (k = 0; k < categs; k++) tbl[j][k] = rrate[j]*rate[k]; sumrates = 0.0; for (i = 0; i < endsite; i++) { for (j = 0; j < rcategs; j++) sumrates += aliasweight[i] * probcat[j] * tbl[j][category[alias[i] - 1] - 1]; } sumrates /= (double)sites; for (j = 0; j < rcategs; j++) for (k = 0; k < categs; k++) { tbl[j][k] /= sumrates; } if(jumb > 1) return; if (gama || invar) { fprintf(outfile, "\nDiscrete approximation to gamma distributed rates\n"); fprintf(outfile, " Coefficient of variation of rates = %f (alpha = %f)\n", cv, alpha); } if (rcategs > 1) { fprintf(outfile, "\nState in HMM Rate of change Probability\n\n"); for (i = 0; i < rcategs; i++) if (probcat[i] < 0.0001) fprintf(outfile, "%9ld%16.3f%20.6f\n", i+1, rrate[i], probcat[i]); else if (probcat[i] < 0.001) fprintf(outfile, "%9ld%16.3f%19.5f\n", i+1, rrate[i], probcat[i]); else if (probcat[i] < 0.01) fprintf(outfile, "%9ld%16.3f%18.4f\n", i+1, rrate[i], probcat[i]); else fprintf(outfile, "%9ld%16.3f%17.3f\n", i+1, rrate[i], probcat[i]); putc('\n', outfile); if (auto_) { fprintf(outfile, "Expected length of a patch of sites having the same rate = %8.3f\n", 1/lambda); putc('\n', outfile); } } if (categs > 1) { fprintf(outfile, "\nSite category Rate of change\n\n"); for (k = 0; k < categs; k++) fprintf(outfile, "%9ld%16.3f\n", k+1, rate[k]); fprintf(outfile, "\n\n"); } } /* prot_inittable */ void free_pmatrix(long sib) { long j,k,l; for (j = 0; j < rcategs; j++) { for (k = 0; k < categs; k++) { for (l = 0; l < 20; l++) free(pmatrices[sib][j][k][l]); free(pmatrices[sib][j][k]); } free(pmatrices[sib][j]); } free(pmatrices[sib]); } void alloc_pmatrix(long sib) { /* Allocate memory for a new pmatrix. Called iff num_sibs>max_num_sibs */ long j, k, l; double ****temp_matrix; temp_matrix = (double ****) Malloc (rcategs * sizeof(double ***)); for (j = 0; j < rcategs; j++) { temp_matrix[j] = (double ***) Malloc(categs * sizeof(double **)); for (k = 0; k < categs; k++) { temp_matrix[j][k] = (double **) Malloc(20 * sizeof (double *)); for (l = 0; l < 20; l++) temp_matrix[j][k][l] = (double *) Malloc(20 * sizeof(double)); } } pmatrices[sib] = temp_matrix; max_num_sibs++; } /* alloc_pmatrix */ void make_pmatrix(double **matrix, double **dmat, double **ddmat, long derivative, double lz, double rat, double *eigmat, double **probmat) { /* Computes the R matrix such that matrix[m][l] is the joint probability */ /* of m and l. */ /* Computes a P matrix such that matrix[m][l] is the conditional */ /* probability of m given l. This is accomplished by dividing all terms */ /* in the R matrix by freqaa[m], the frequency of l. */ long k, l, m; /* (l) original character state */ /* (m) final character state */ /* (k) lambda counter */ double p0, p1, p2, q; double elambdat[20], delambdat[20], ddelambdat[20]; /* exponential term for matrix */ /* and both derivative matrices */ for (k = 0; k <= 19; k++) { elambdat[k] = exp(lz * rat * eigmat[k]); if(derivative != 0) { delambdat[k] = (elambdat[k] * rat * eigmat[k]); ddelambdat[k] = (delambdat[k] * rat * eigmat[k]); } } for (m = 0; m <= 19; m++) { for (l = 0; l <= 19; l++) { p0 = 0.0; p1 = 0.0; p2 = 0.0; for (k = 0; k <= 19; k++) { q = probmat[k][m] * probmat[k][l]; p0 += (q * elambdat[k]); if(derivative !=0) { p1 += (q * delambdat[k]); p2 += (q * ddelambdat[k]); } } matrix[m][l] = p0 / freqaa[m]; if(derivative != 0) { dmat[m][l] = p1 / freqaa[m]; ddmat[m][l] = p2 / freqaa[m]; } } } } /* make_pmatrix */ void prot_nuview(node *p) { long b, i, j, k, l, m, num_sibs, sib_index; node *sib_ptr, *sib_back_ptr; psitelike prot_xx, x2; double lw, prod7; double **pmat; double maxx,correction; /* Figure out how many siblings the current node has */ /* and be sure that pmatrices is large enough */ num_sibs = count_sibs(p); for (i = 0; i < num_sibs; i++) if (pmatrices[i] == NULL) alloc_pmatrix(i); /* Recursive calls, should be called for all children */ sib_ptr = p; for (i=0 ; i < num_sibs; i++) { sib_ptr = sib_ptr->next; sib_back_ptr = sib_ptr->back; if (!(sib_back_ptr == NULL)) if (!sib_back_ptr->tip && !sib_back_ptr->initialized) prot_nuview(sib_back_ptr); } /* Make pmatrices for all possible combinations of category, rcateg */ /* and sib */ sib_ptr = p; /* return to p */ for (sib_index=0; sib_index < num_sibs; sib_index++) { sib_ptr = sib_ptr->next; sib_back_ptr = sib_ptr->back; if (sib_back_ptr != NULL) lw = fabs(p->tyme - sib_back_ptr->tyme); else lw = 0.0; for (j = 0; j < rcategs; j++) for (k = 0; k < categs; k++) make_pmatrix(pmatrices[sib_index][j][k], NULL, NULL, 0, lw, tbl[j][k], eigmat, probmat); } for (i = 0; i < endsite; i++) { correction = 0; maxx = 0; k = category[alias[i]-1] - 1; for (j = 0; j < rcategs; j++) { /* initialize to 1 all values of prot_xx */ for (m = 0; m <= 19; m++) prot_xx[m] = 1; sib_ptr = p; /* return to p */ /* loop through all sibs and calculate likelihoods for all possible*/ /* amino acid combinations */ for (sib_index=0; sib_index < num_sibs; sib_index++) { sib_ptr = sib_ptr->next; sib_back_ptr = sib_ptr->back; if (sib_back_ptr != NULL) { memcpy(x2, sib_back_ptr->protx[i][j], sizeof(psitelike)); if ( j == 0 ) correction += sib_back_ptr->underflows[i]; } else for (b = 0; b <= 19; b++) x2[b] = 1.0; pmat = pmatrices[sib_index][j][k]; for (m = 0; m <= 19; m++) { prod7 = 0; for (l = 0; l <= 19; l++) prod7 += (pmat[m][l] * x2[l]); prot_xx[m] *= prod7; if ( prot_xx[m] > maxx && sib_index == (num_sibs - 1 )) maxx = prot_xx[m]; } } /* And the final point of this whole function: */ memcpy(p->protx[i][j], prot_xx, sizeof(psitelike)); } p->underflows[i] = 0; if ( maxx < MIN_DOUBLE ) fix_protx(p,i,maxx,rcategs); p->underflows[i] += correction; } p->initialized = true; } /* prot_nuview */ void getthree(node *p, double thigh, double tlow) { /* compute likelihood at a new triple of points */ int i; double tt = p->tyme; double td = fabs(tdelta); x[0] = tt - td; x[1] = tt; x[2] = tt + td; if ( x[0] < tlow + epsilon ) { x[0] = tlow + epsilon; x[1] = ( x[0] + x[2] ) / 2; } if ( x[2] > thigh - epsilon ) { x[2] = thigh - epsilon; x[1] = ( x[0] + x[2] ) / 2; } for ( i = 0 ; i < 3 ; i++ ) { p->tyme = x[i]; prot_nuview(p); lnl[i] = prot_evaluate(p); } } /* getthree */ void makenewv(node *p) { /* improve a node time */ long it, imin, imax, i, num_sibs; double tt, tfactor, tlow, thigh, oldlike, ymin, ymax, s32, s21, yold; boolean done, already; node *s, *sdown, *sib_ptr, *sib_back_ptr; s = curtree.nodep[p->index - 1]; sdown = s->back; if (s == curtree.root) tlow = -10.0; else tlow = sdown->tyme; sib_ptr = s; num_sibs = count_sibs(p); thigh = s->next->back->tyme; for (i=0 ; i < num_sibs; i++) { sib_ptr = sib_ptr->next; sib_back_ptr = sib_ptr->back; if (sib_back_ptr->tyme < thigh) thigh = sib_back_ptr->tyme; } done = (thigh - tlow < 4.0*epsilon); it = 1; if (s != curtree.root) tdelta = (thigh - tlow) / 10.0; else tdelta = (thigh - s->tyme) / 5.0; tfactor = 1.0; if (!done) getthree(s, thigh, tlow); while (it < iterations && !done) { ymax = lnl[0]; imax = 1; for (i = 2; i <= 3; i++) { if (lnl[i - 1] > ymax) { ymax = lnl[i - 1]; imax = i; } } if (imax != 2) { ymax = x[1]; x[1] = x[imax - 1]; x[imax - 1] = ymax; ymax = lnl[1]; lnl[1] = lnl[imax - 1]; lnl[imax - 1] = ymax; } tt = x[1]; oldlike = lnl[1]; yold = tt; s32 = (lnl[2] - lnl[1]) / (x[2] - x[1]); s21 = (lnl[1] - lnl[0]) / (x[1] - x[0]); if (fabs(x[2] - x[0]) > epsilon) curv = (s32 - s21) / ((x[2] - x[0]) / 2); else curv = 0.0; slope = (s32 + s21) / 2 - curv * (x[2] - 2 * x[1] + x[0]) / 4; if (curv >= 0.0) { if (slope < 0) tdelta = -fabs(tdelta); else tdelta = fabs(tdelta); } else tdelta = -(tfactor * slope / curv); if (tt + tdelta <= tlow + epsilon) tdelta = tlow + epsilon - tt; if (tt + tdelta >= thigh - epsilon) tdelta = thigh - epsilon - tt; tt += tdelta; done = (fabs(yold - tt) < epsilon || fabs(tdelta) < epsilon); s->tyme = tt; prot_nuview(s); lnlike = prot_evaluate(s); ymin = lnl[0]; imin = 1; for (i = 2; i <= 3; i++) { if (lnl[i - 1] < ymin) { ymin = lnl[i - 1]; imin = i; } } already = (tt == x[0]) || (tt == x[1]) || (tt == x[2]); if (!already && ymin < lnlike) { x[imin - 1] = tt; lnl[imin - 1] = lnlike; } if (already || lnlike < oldlike) { tt = x[1]; tfactor /= 2; tdelta /= 2; curtree.likelihood = oldlike; lnlike = oldlike; } else tfactor = 1.0; if (!done) { sib_ptr = p; num_sibs = count_sibs(p); p->tyme = tt; for (i=0 ; i < num_sibs; i++) { sib_ptr = sib_ptr->next; sib_ptr->tyme = tt; } sib_ptr = p; prot_nuview(p); for (i=0 ; i < num_sibs; i++) { sib_ptr = sib_ptr->next; prot_nuview(sib_ptr); } } it++; } sib_ptr = p; for (i=0 ; i < num_sibs; i++) { sib_ptr = sib_ptr->next; inittrav (sib_ptr); } smoothed = smoothed && done; } /* makenewv */ void update(node *p) { node *sib_ptr, *sib_back_ptr; long i, num_sibs; /* improve time and recompute views at a node */ if (p == NULL) return; if (p->back != NULL) { if (!p->back->tip && !p->back->initialized) prot_nuview(p->back); } sib_ptr = p; num_sibs = count_sibs(p); for (i=0 ; i < num_sibs; i++) { sib_ptr = sib_ptr->next; sib_back_ptr = sib_ptr->back; if (sib_back_ptr != NULL) { if (!sib_back_ptr->tip && !sib_back_ptr->initialized) prot_nuview(sib_back_ptr); } } if ((!usertree) || (usertree && !lngths) || p->iter) { makenewv(p); return; } prot_nuview(p); sib_ptr = p; num_sibs = count_sibs(p); for (i=0 ; i < num_sibs; i++) { sib_ptr = sib_ptr->next; prot_nuview(sib_ptr); } } /* update */ void smooth(node *p) { node *sib_ptr; long i, num_sibs; if (p == NULL) return; if (p->tip) return; update(p); smoothed = false; sib_ptr = p; num_sibs = count_sibs(p); for (i=0; i < num_sibs; i++) { sib_ptr = sib_ptr->next; if (polishing || (smoothit && !smoothed)) { smooth(sib_ptr->back); p->initialized = false; sib_ptr->initialized = false; } update(p); } } /* smooth */ void promlk_add(node *below, node *newtip, node *newfork, boolean tempadd) { /* inserts the nodes newfork and its descendant, newtip, into the tree. */ long i; boolean done; node *p; below = curtree.nodep[below->index - 1]; newfork = curtree.nodep[newfork->index-1]; newtip = curtree.nodep[newtip->index-1]; if (below->back != NULL) below->back->back = newfork; newfork->back = below->back; below->back = newfork->next->next; newfork->next->next->back = below; newfork->next->back = newtip; newtip->back = newfork->next; if (newtip->tyme < below->tyme) p = newtip; else p = below; newfork->tyme = p->tyme; if (curtree.root == below) curtree.root = newfork; if (newfork->back != NULL) { if (p->tyme > newfork->back->tyme) newfork->tyme = (p->tyme + newfork->back->tyme) / 2.0; else newfork->tyme = p->tyme - epsilon; newfork->next->tyme = newfork->tyme; newfork->next->next->tyme = newfork->tyme; do { p = curtree.nodep[p->back->index - 1]; done = (p == curtree.root); if (!done) done = (curtree.nodep[p->back->index - 1]->tyme < p->tyme - epsilon); if (!done) { curtree.nodep[p->back->index - 1]->tyme = p->tyme - epsilon; curtree.nodep[p->back->index - 1]->next->tyme = p->tyme - epsilon; curtree.nodep[p->back->index - 1]->next->next->tyme = p->tyme - epsilon; } } while (!done); } else { newfork->tyme = newfork->tyme - 2*epsilon; newfork->next->tyme = newfork->tyme; newfork->next->next->tyme = newfork->tyme; } inittrav(newtip); inittrav(newtip->back); smoothed = false; i = 1; while (i < smoothings && !smoothed) { smoothed = true; smooth(newfork); smooth(newfork->back); i++; } } /* promlk_add */ void promlk_re_move(node **item, node **fork, boolean tempadd) { /* removes nodes item and its ancestor, fork, from the tree. the new descendant of fork's ancestor is made to be fork's second descendant (other than item). Also returns pointers to the deleted nodes, item and fork */ node *p, *q; long i; if ((*item)->back == NULL) { *fork = NULL; return; } *item = curtree.nodep[(*item)->index-1]; *fork = curtree.nodep[(*item)->back->index - 1]; if (curtree.root == *fork) { if (*item == (*fork)->next->back) curtree.root = (*fork)->next->next->back; else curtree.root = (*fork)->next->back; } p = (*item)->back->next->back; q = (*item)->back->next->next->back; if (p != NULL) p->back = q; if (q != NULL) q->back = p; (*fork)->back = NULL; p = (*fork)->next; while (p != *fork) { p->back = NULL; p = p->next; } (*item)->back = NULL; inittrav(p); inittrav(q); if (tempadd) return; i = 1; while (i <= smoothings) { smooth(q); if (smoothit) smooth(q->back); i++; } } /* promlk_re_move */ double prot_evaluate(node *p) { contribarr tterm; static contribarr like, nulike, clai; double sum, sum2, sumc=0, y, prod4, prodl, frexm, sumterm, lterm; double **pmat; long i, j, k, l, m, lai; node *q, *r; psitelike x1, x2; sum = 0.0; if (p == curtree.root && (count_sibs(p) == 2)) { r = p->next->back; q = p->next->next->back; y = r->tyme + q->tyme - 2 * p->tyme; if (!r->tip && !r->initialized) prot_nuview (r); if (!q->tip && !q->initialized) prot_nuview (q); } else if (p == curtree.root) { /* the next two lines copy tyme and x to p->next. Normally they are not initialized for an internal node. */ /* assumes bifurcation */ p->next->tyme = p->tyme; prot_nuview(p->next); r = p->next; q = p->next->back; y = fabs(p->next->tyme - q->tyme); } else { r = p; q = p->back; if (!r->tip && !r->initialized) prot_nuview (r); if (!q->tip && !q->initialized) prot_nuview (q); y = fabs(r->tyme - q->tyme); } for (j = 0; j < rcategs; j++) for (k = 0; k < categs; k++) make_pmatrix(pmatrices[0][j][k],NULL,NULL,0,y,tbl[j][k],eigmat,probmat); for (i = 0; i < endsite; i++) { k = category[alias[i]-1] - 1; for (j = 0; j < rcategs; j++) { memcpy(x1, r->protx[i][j], sizeof(psitelike)); memcpy(x2, q->protx[i][j], sizeof(psitelike)); prod4 = 0.0; pmat = pmatrices[0][j][k]; for (m = 0; m <= 19; m++) { prodl = 0.0; for (l = 0; l <= 19; l++) prodl += (pmat[m][l] * x2[l]); frexm = x1[m] * freqaa[m]; prod4 += (prodl * frexm); } tterm[j] = prod4; } sumterm = 0.0; for (j = 0; j < rcategs; j++) sumterm += probcat[j] * tterm[j]; if (sumterm < 0.0) sumterm = 0.00000001; /* ??? */ lterm = log(sumterm) + p->underflows[i] + q->underflows[i]; for (j = 0; j < rcategs; j++) clai[j] = tterm[j] / sumterm; memcpy(contribution[i], clai, rcategs * sizeof(double)); if (!auto_ && usertree && (which <= shimotrees)) l0gf[which - 1][i] = lterm; sum += aliasweight[i] * lterm; } if (auto_) { for (j = 0; j < rcategs; j++) like[j] = 1.0; for (i = 0; i < sites; i++) { sumc = 0.0; for (k = 0; k < rcategs; k++) sumc += probcat[k] * like[k]; sumc *= lambda; if ((ally[i] > 0) && (location[ally[i]-1] > 0)) { lai = location[ally[i] - 1]; memcpy(clai, contribution[lai - 1], rcategs*sizeof(double)); for (j = 0; j < rcategs; j++) nulike[j] = ((1.0 - lambda) * like[j] + sumc) * clai[j]; } else { for (j = 0; j < rcategs; j++) nulike[j] = ((1.0 - lambda) * like[j] + sumc); } memcpy(like, nulike, rcategs * sizeof(double)); } sum2 = 0.0; for (i = 0; i < rcategs; i++) sum2 += probcat[i] * like[i]; sum += log(sum2); } curtree.likelihood = sum; if (auto_ || !usertree) return sum; if(which <= shimotrees) l0gl[which - 1] = sum; if (which == 1) { maxwhich = 1; maxlogl = sum; return sum; } if (sum > maxlogl) { maxwhich = which; maxlogl = sum; } return sum; } /* prot_evaluate */ void tryadd(node *p, node **item, node **nufork) { /* temporarily adds one fork and one tip to the tree. if the location where they are added yields greater likelihood than other locations tested up to that time, then keeps that location as there */ long grcategs; grcategs = (categs > rcategs) ? categs : rcategs; promlk_add(p, *item, *nufork, true); like = prot_evaluate(p); if (lastsp) { if (like >= bestyet || bestyet == UNDEFINED) prot_copy_(&curtree, &bestree, nonodes, grcategs); } if (like > bestyet || bestyet == UNDEFINED) { bestyet = like; there = p; } promlk_re_move(item, nufork, true); } /* tryadd */ void addpreorder(node *p, node *item_, node *nufork_, boolean contin, boolean continagain) { /* traverses a binary tree, calling function tryadd at a node before calling tryadd at its descendants */ node *item, *nufork; item = item_; nufork = nufork_; if (p == NULL) return; tryadd(p, &item, &nufork); contin = continagain; if ((!p->tip) && contin) { addpreorder(p->next->back, item, nufork, contin, continagain); addpreorder(p->next->next->back, item, nufork, contin, continagain); } } /* addpreorder */ void restoradd(node *below, node *newtip, node *newfork, double prevtyme) { /* restore "new" tip and fork to place "below". restore tymes */ /* assumes bifurcation */ hookup(newfork, below->back); hookup(newfork->next, below); hookup(newtip, newfork->next->next); curtree.nodep[newfork->index-1] = newfork; newfork->tyme = prevtyme; /* assumes bifurcations */ newfork->next->tyme = prevtyme; newfork->next->next->tyme = prevtyme; } /* restoradd */ void tryrearr(node *p, boolean *success) { /* evaluates one rearrangement of the tree. if the new tree has greater likelihood than the old one sets success = TRUE and keeps the new tree. otherwise, restores the old tree */ node *frombelow, *whereto, *forknode; double oldlike, prevtyme; boolean wasonleft; if (p == curtree.root) return; forknode = curtree.nodep[p->back->index - 1]; if (forknode == curtree.root) return; oldlike = bestyet; prevtyme = forknode->tyme; /* the following statement presumes bifurcating tree */ if (forknode->next->back == p) { frombelow = forknode->next->next->back; wasonleft = true; } else { frombelow = forknode->next->back; wasonleft = false; } whereto = curtree.nodep[forknode->back->index - 1]; promlk_re_move(&p, &forknode, true); promlk_add(whereto, p, forknode, true); like = prot_evaluate(p); if (like <= oldlike && oldlike != UNDEFINED) { promlk_re_move(&p, &forknode, true); restoradd(frombelow, p, forknode, prevtyme); if (wasonleft && (forknode->next->next->back == p)) { hookup (forknode->next->back, forknode->next->next); hookup (forknode->next, p); } curtree.likelihood = oldlike; inittrav(forknode); inittrav(forknode->next); inittrav(forknode->next->next); } else { (*success) = true; bestyet = like; } } /* tryrearr */ void repreorder(node *p, boolean *success) { /* traverses a binary tree, calling function tryrearr at a node before calling tryrearr at its descendants */ if (p == NULL) return; tryrearr(p, success); if (p->tip) return; if (!(*success)) repreorder(p->next->back, success); if (!(*success)) repreorder(p->next->next->back, success); } /* repreorder */ void rearrange(node **r) { /* traverses the tree (preorder), finding any local rearrangement which increases the likelihood. if traversal succeeds in increasing the tree's likelihood, function rearrange runs traversal again */ boolean success; success = true; while (success) { success = false; repreorder(*r, &success); } } /* rearrange */ void nodeinit(node *p) { /* set up times at one node */ node *sib_ptr, *sib_back_ptr; long i, num_sibs; double lowertyme; sib_ptr = p; num_sibs = count_sibs(p); /* lowertyme = lowest of children's times */ lowertyme = p->next->back->tyme; for (i=0 ; i < num_sibs; i++) { sib_ptr = sib_ptr->next; sib_back_ptr = sib_ptr->back; if (sib_back_ptr->tyme < lowertyme) lowertyme = sib_back_ptr->tyme; } p->tyme = lowertyme - 0.1; sib_ptr = p; for (i=0 ; i < num_sibs; i++) { sib_ptr = sib_ptr->next; sib_back_ptr = sib_ptr->back; sib_ptr->tyme = p->tyme; sib_back_ptr->v = sib_back_ptr->tyme - p->tyme; sib_ptr->v = sib_back_ptr->v; } } /* nodeinit */ void initrav(node *p) { long i, num_sibs; node *sib_ptr, *sib_back_ptr; /* traverse to set up times throughout tree */ if (p->tip) return; sib_ptr = p; num_sibs = count_sibs(p); for (i=0 ; i < num_sibs; i++) { sib_ptr = sib_ptr->next; sib_back_ptr = sib_ptr->back; initrav(sib_back_ptr); } nodeinit(p); } /* initrav */ void travinit(node *p) { long i, num_sibs; node *sib_ptr, *sib_back_ptr; /* traverse to set up initial values */ if (p == NULL) return; if (p->tip) return; if (p->initialized) return; sib_ptr = p; num_sibs = count_sibs(p); for (i=0 ; i < num_sibs; i++) { sib_ptr = sib_ptr->next; sib_back_ptr = sib_ptr->back; travinit(sib_back_ptr); } prot_nuview(p); p->initialized = true; } /* travinit */ void travsp(node *p) { long i, num_sibs; node *sib_ptr, *sib_back_ptr; /* traverse to find tips */ if (p == curtree.root) travinit(p); if (p->tip) travinit(p->back); else { sib_ptr = p; num_sibs = count_sibs(p); for (i=0 ; i < num_sibs; i++) { sib_ptr = sib_ptr->next; sib_back_ptr = sib_ptr->back; travsp(sib_back_ptr); } } } /* travsp */ void treevaluate() { /* evaluate likelihood of tree, after iterating branch lengths */ long i, j, num_sibs; node *sib_ptr; polishing = true; smoothit = true; for (i = 0; i < spp; i++) curtree.nodep[i]->initialized = false; for (i = spp; i < nonodes; i++) { sib_ptr = curtree.nodep[i]; sib_ptr->initialized = false; num_sibs = count_sibs(sib_ptr); for (j=0 ; j < num_sibs; j++) { sib_ptr = sib_ptr->next; sib_ptr->initialized = false; } } if (!lngths) initrav(curtree.root); travsp(curtree.root); for (i = 1; i <= smoothings * 4; i++) smooth(curtree.root); prot_evaluate(curtree.root); } /* treevaluate */ void promlk_coordinates(node *p, long *tipy) { /* establishes coordinates of nodes */ node *q, *first, *last, *pp1 =NULL, *pp2 =NULL; long num_sibs, p1, p2, i; if (p->tip) { p->xcoord = 0; p->ycoord = (*tipy); p->ymin = (*tipy); p->ymax = (*tipy); (*tipy) += down; return; } q = p->next; do { promlk_coordinates(q->back, tipy); q = q->next; } while (p != q); num_sibs = count_sibs(p); p1 = (long)((num_sibs+1)/2.0); p2 = (long)((num_sibs+2)/2.0); i = 1; q = p->next; first = q->back; do { if (i == p1) pp1 = q->back; if (i == p2) pp2 = q->back; last = q->back; q = q->next; i++; } while (q != p); p->xcoord = (long)(0.5 - over * p->tyme); p->ycoord = (pp1->ycoord + pp2->ycoord) / 2; p->ymin = first->ymin; p->ymax = last->ymax; } /* promlk_coordinates */ void promlk_drawline(long i, double scale) { /* draws one row of the tree diagram by moving up tree */ node *p, *q, *r, *first =NULL, *last =NULL; long n, j; boolean extra, done; p = curtree.root; q = curtree.root; extra = false; if ((long)(p->ycoord) == i) { if (p->index - spp >= 10) fprintf(outfile, "-%2ld", p->index - spp); else fprintf(outfile, "--%ld", p->index - spp); extra = true; } else fprintf(outfile, " "); do { if (!p->tip) { r = p->next; done = false; do { if (i >= r->back->ymin && i <= r->back->ymax) { q = r->back; done = true; } r = r->next; } while (!(done || r == p)); first = p->next->back; r = p->next; while (r->next != p) r = r->next; last = r->back; } done = (p == q); n = (long)(scale * ((long)(p->xcoord) - (long)(q->xcoord)) + 0.5); if (n < 3 && !q->tip) n = 3; if (extra) { n--; extra = false; } if ((long)(q->ycoord) == i && !done) { if (p->ycoord != q->ycoord) putc('+', outfile); else putc('-', outfile); if (!q->tip) { for (j = 1; j <= n - 2; j++) putc('-', outfile); if (q->index - spp >= 10) fprintf(outfile, "%2ld", q->index - spp); else fprintf(outfile, "-%ld", q->index - spp); extra = true; } else { for (j = 1; j < n; j++) putc('-', outfile); } } else if (!p->tip) { if ((long)(last->ycoord) > i && (long)(first->ycoord) < i && i != (long)(p->ycoord)) { putc('!', outfile); for (j = 1; j < n; j++) putc(' ', outfile); } else { for (j = 1; j <= n; j++) putc(' ', outfile); } } else { for (j = 1; j <= n; j++) putc(' ', outfile); } if (p != q) p = q; } while (!done); if ((long)(p->ycoord) == i && p->tip) { for (j = 0; j < nmlngth; j++) putc(nayme[p->index - 1][j], outfile); } putc('\n', outfile); } /* promlk_drawline */ void promlk_printree() { /* prints out diagram of the tree */ long tipy; double scale; long i; node *p; if (!treeprint) return; putc('\n', outfile); tipy = 1; promlk_coordinates(curtree.root, &tipy); p = curtree.root; while (!p->tip) p = p->next->back; scale = 1.0 / (long)(p->tyme - curtree.root->tyme + 1.000); putc('\n', outfile); for (i = 1; i <= tipy - down; i++) promlk_drawline(i, scale); putc('\n', outfile); } /* promlk_printree */ void describe(node *p) { long i, num_sibs; node *sib_ptr, *sib_back_ptr; double v; if (p == curtree.root) fprintf(outfile, " root "); else fprintf(outfile, "%4ld ", p->back->index - spp); if (p->tip) { for (i = 0; i < nmlngth; i++) putc(nayme[p->index - 1][i], outfile); } else fprintf(outfile, "%4ld ", p->index - spp); if (p != curtree.root) { fprintf(outfile, "%11.5f", (p->tyme - curtree.root->tyme)); v = (p->tyme - curtree.nodep[p->back->index - 1]->tyme); fprintf(outfile, "%13.5f", v); } putc('\n', outfile); if (!p->tip) { sib_ptr = p; num_sibs = count_sibs(p); for (i=0 ; i < num_sibs; i++) { sib_ptr = sib_ptr->next; sib_back_ptr = sib_ptr->back; describe(sib_back_ptr); } } } /* describe */ void prot_reconstr(node *p, long n) { /* reconstruct and print out acid at site n+1 at node p */ long i, j, k, first, num_sibs = 0; double f, sum, xx[20]; node *q = NULL; if (p->tip) putc(y[p->index-1][n], outfile); else { num_sibs = count_sibs(p); if ((ally[n] == 0) || (location[ally[n]-1] == 0)) putc('.', outfile); else { j = location[ally[n]-1] - 1; sum = 0; for (i = 0; i <= 19; i++) { f = p->protx[j][mx-1][i]; if (!p->tip) { q = p; for (k = 0; k < num_sibs; k++) { q = q->next; f *= q->protx[j][mx-1][i]; } } f = sqrt(f); xx[i] = f * freqaa[i]; sum += xx[i]; } for (i = 0; i <= 19; i++) xx[i] /= sum; first = 0; for (i = 0; i <= 19; i++) if (xx[i] > xx[first]) first = i; if (xx[first] > 0.95) putc(aachar[first], outfile); else putc(tolower(aachar[first]), outfile); if (rctgry && rcategs > 1) mx = mp[n][mx - 1]; else mx = 1; } } } /* prot_reconstr */ void rectrav(node *p, long m, long n) { /* print out segment of reconstructed sequence for one branch */ long num_sibs, i; node *sib_ptr; putc(' ', outfile); if (p->tip) { for (i = 0; i < nmlngth; i++) putc(nayme[p->index-1][i], outfile); } else fprintf(outfile, "%4ld ", p->index - spp); fprintf(outfile, " "); mx = mx0; for (i = m; i <= n; i++) { if ((i % 10 == 0) && (i != m)) putc(' ', outfile); prot_reconstr(p, i); } putc('\n', outfile); if (!p->tip) { num_sibs = count_sibs(p); sib_ptr = p; for (i = 0; i < num_sibs; i++) { sib_ptr = sib_ptr->next; rectrav(sib_ptr->back, m, n); } } mx1 = mx; } /* rectrav */ void summarize() { long i, j, mm; double mode, sum; double like[maxcategs], nulike[maxcategs]; double **marginal; mp = (long **)Malloc(sites * sizeof(long *)); for (i = 0; i <= sites-1; ++i) mp[i] = (long *)Malloc(sizeof(long)*rcategs); fprintf(outfile, "\nLn Likelihood = %11.5f\n\n", curtree.likelihood); fprintf(outfile, " Ancestor Node Node Height Length\n"); fprintf(outfile, " -------- ---- ---- ------ ------\n"); describe(curtree.root); putc('\n', outfile); if (rctgry && rcategs > 1) { for (i = 0; i < rcategs; i++) like[i] = 1.0; for (i = sites - 1; i >= 0; i--) { sum = 0.0; for (j = 0; j < rcategs; j++) { nulike[j] = (lambda1 + lambda * probcat[j]) * like[j]; mp[i][j] = j + 1; for (k = 1; k <= rcategs; k++) { if (k != j + 1) { if (lambda * probcat[k - 1] * like[k - 1] > nulike[j]) { nulike[j] = lambda * probcat[k - 1] * like[k - 1]; mp[i][j] = k; } } } if ((ally[i] > 0) && (location[ally[i]-1] > 0)) nulike[j] *= contribution[location[ally[i] - 1] - 1][j]; sum += nulike[j]; } for (j = 0; j < rcategs; j++) nulike[j] /= sum; memcpy(like, nulike, rcategs * sizeof(double)); } mode = 0.0; mx = 1; for (i = 1; i <= rcategs; i++) { if (probcat[i - 1] * like[i - 1] > mode) { mx = i; mode = probcat[i - 1] * like[i - 1]; } } mx0 = mx; fprintf(outfile, "Combination of categories that contributes the most to the likelihood:\n\n"); for (i = 1; i <= nmlngth + 3; i++) putc(' ', outfile); for (i = 1; i <= sites; i++) { fprintf(outfile, "%ld", mx); if (i % 10 == 0) putc(' ', outfile); if (i % 60 == 0 && i != sites) { putc('\n', outfile); for (j = 1; j <= nmlngth + 3; j++) putc(' ', outfile); } mx = mp[i - 1][mx - 1]; } fprintf(outfile, "\n\n"); marginal = (double **) Malloc( sites*sizeof(double *)); for (i = 0; i < sites; i++) marginal[i] = (double *) Malloc( rcategs*sizeof(double)); for (i = 0; i < rcategs; i++) like[i] = 1.0; for (i = sites - 1; i >= 0; i--) { sum = 0.0; for (j = 0; j < rcategs; j++) { nulike[j] = (lambda1 + lambda * probcat[j]) * like[j]; for (k = 1; k <= rcategs; k++) { if (k != j + 1) nulike[j] += lambda * probcat[k - 1] * like[k - 1]; } if ((ally[i] > 0) && (location[ally[i]-1] > 0)) nulike[j] *= contribution[location[ally[i] - 1] - 1][j]; sum += nulike[j]; } for (j = 0; j < rcategs; j++) { nulike[j] /= sum; marginal[i][j] = nulike[j]; } memcpy(like, nulike, rcategs * sizeof(double)); } for (i = 0; i < rcategs; i++) like[i] = 1.0; for (i = 0; i < sites; i++) { sum = 0.0; for (j = 0; j < rcategs; j++) { nulike[j] = (lambda1 + lambda * probcat[j]) * like[j]; for (k = 1; k <= rcategs; k++) { if (k != j + 1) nulike[j] += lambda * probcat[k - 1] * like[k - 1]; } marginal[i][j] *= like[j] * probcat[j]; sum += nulike[j]; } for (j = 0; j < rcategs; j++) nulike[j] /= sum; memcpy(like, nulike, rcategs * sizeof(double)); sum = 0.0; for (j = 0; j < rcategs; j++) sum += marginal[i][j]; for (j = 0; j < rcategs; j++) marginal[i][j] /= sum; } fprintf(outfile, "Most probable category at each site if > 0.95"); fprintf(outfile, " probability (\".\" otherwise)\n\n"); for (i = 1; i <= nmlngth + 3; i++) putc(' ', outfile); for (i = 0; i < sites; i++) { sum = 0.0; for (j = 0; j < rcategs; j++) if (marginal[i][j] > sum) { sum = marginal[i][j]; mm = j; } if (sum >= 0.95) fprintf(outfile, "%ld", mm+1); else putc('.', outfile); if ((i+1) % 60 == 0) { if (i != 0) { putc('\n', outfile); for (j = 1; j <= nmlngth + 3; j++) putc(' ', outfile); } } else if ((i+1) % 10 == 0) putc(' ', outfile); } putc('\n', outfile); for (i = 0; i < sites; i++) free(marginal[i]); free(marginal); } putc('\n', outfile); putc('\n', outfile); putc('\n', outfile); if (hypstate) { fprintf(outfile, "Probable sequences at interior nodes:\n\n"); fprintf(outfile, " node "); for (i = 0; (i < 13) && (i < ((sites + (sites-1)/10 - 39) / 2)); i++) putc(' ', outfile); fprintf(outfile, "Reconstructed sequence (caps if > 0.95)\n\n"); if (!rctgry || (rcategs == 1)) mx0 = 1; for (i = 0; i < sites; i += 60) { k = i + 59; if (k >= sites) k = sites - 1; rectrav(curtree.root, i, k); putc('\n', outfile); mx0 = mx1; } } for (i = 0; i <= sites-1; ++i) free(mp[i]); free(mp); } /* summarize */ void promlk_treeout(node *p) { /* write out file with representation of final tree */ node *sib_ptr; long i, n, w, num_sibs; Char c; double x; if (p->tip) { n = 0; for (i = 1; i <= nmlngth; i++) { if (nayme[p->index - 1][i - 1] != ' ') n = i; } for (i = 0; i < n; i++) { c = nayme[p->index - 1][i]; if (c == ' ') c = '_'; putc(c, outtree); } col += n; } else { sib_ptr = p; num_sibs = count_sibs(p); putc('(', outtree); col++; for (i=0; i < (num_sibs - 1); i++) { sib_ptr = sib_ptr->next; promlk_treeout(sib_ptr->back); putc(',', outtree); col++; if (col > 55) { putc('\n', outtree); col = 0; } } sib_ptr = sib_ptr->next; promlk_treeout(sib_ptr->back); putc(')', outtree); col++; } if (p == curtree.root) { fprintf(outtree, ";\n"); return; } x = (p->tyme - curtree.nodep[p->back->index - 1]->tyme); if (x > 0.0) w = (long)(0.4342944822 * log(x)); else if (x == 0.0) w = 0; else w = (long)(0.4342944822 * log(-x)) + 1; if (w < 0) w = 0; fprintf(outtree, ":%*.5f", (int)(w + 7), x); col += w + 8; } /* promlk_treeout */ void initpromlnode(node **p, node **grbg, node *q, long len, long nodei, long *ntips, long *parens, initops whichinit, pointarray treenode, pointarray nodep, Char *str, Char *ch, FILE *intree) { /* initializes a node */ boolean minusread; double valyew, divisor; switch (whichinit) { case bottom: gnu(grbg, p); (*p)->index = nodei; (*p)->tip = false; malloc_ppheno((*p), endsite, rcategs); nodep[(*p)->index - 1] = (*p); break; case nonbottom: gnu(grbg, p); malloc_ppheno(*p, endsite, rcategs); (*p)->index = nodei; break; case tip: match_names_to_data(str, nodep, p, spp); break; case iter: (*p)->initialized = false; (*p)->v = initialv; (*p)->iter = true; if ((*p)->back != NULL) (*p)->back->iter = true; break; case length: processlength(&valyew, &divisor, ch, &minusread, intree, parens); (*p)->v = valyew / divisor; (*p)->iter = false; if ((*p)->back != NULL) { (*p)->back->v = (*p)->v; (*p)->back->iter = false; } break; case unittrwt: curtree.nodep[spp]->iter = false; break; default: /* cases hslength, hsnolength, treewt */ break; /* should never occur */ } } /* initpromlnode */ void tymetrav(node *p, double *x) { /* set up times of nodes */ node *sib_ptr, *q; long i, num_sibs; double xmax; xmax = 0.0; if (!p->tip) { sib_ptr = p; num_sibs = count_sibs(p); for (i=0; i < num_sibs; i++) { sib_ptr = sib_ptr->next; tymetrav(sib_ptr->back, x); if (xmax > (*x)) xmax = (*x); } } else (*x) = 0.0; p->tyme = xmax; if (!p->tip) { q = p; while (q->next != p) { q = q->next; q->tyme = p->tyme; } } (*x) = p->tyme - p->v; } /* tymetrav */ void free_all_protx (long nonodes, pointarray treenode) { /* used in proml */ long i, j, k; node *p; /* Zero thru spp are tips, */ for (i = 0; i < spp; i++) { for (j = 0; j < endsite; j++) free(treenode[i]->protx[j]); free(treenode[i]->protx); } /* The rest are rings (i.e. triads) */ for (i = spp; i < nonodes; i++) { if (treenode[i] != NULL) { p = treenode[i]; for (j = 1; j <= 3; j++) { for (k = 0; k < endsite; k++) free(p->protx[k]); free(p->protx); p = p->next; } } } } /* free_all_protx */ void maketree() { /* constructs a binary tree from the pointers in curtree.nodep, adds each node at location which yields highest likelihood then rearranges the tree for greatest likelihood */ long i, j; long numtrees = 0; double bestlike, gotlike, x; node *item, *nufork, *dummy, *q, *root=NULL; boolean dummy_haslengths, dummy_first, goteof; long nextnode; long grcategs; pointarray dummy_treenode=NULL; grcategs = (categs > rcategs) ? categs : rcategs; prot_inittable(); if (!usertree) { for (i = 1; i <= spp; i++) enterorder[i - 1] = i; if (jumble) randumize(seed, enterorder); curtree.root = curtree.nodep[spp]; curtree.root->back = NULL; for (i = 0; i < spp; i++) curtree.nodep[i]->back = NULL; for (i = spp; i < nonodes; i++) { q = curtree.nodep[i]; q->back = NULL; while ((q = q->next) != curtree.nodep[i]) q->back = NULL; } polishing = false; promlk_add(curtree.nodep[enterorder[0]-1], curtree.nodep[enterorder[1]-1], curtree.nodep[spp], false); if (progress) { printf("\nAdding species:\n"); writename(0, 2, enterorder); #ifdef WIN32 phyFillScreenColor(); #endif } lastsp = false; smoothit = false; for (i = 3; i <= spp; i++) { bestyet = UNDEFINED; bestree.likelihood = bestyet; there = curtree.root; item = curtree.nodep[enterorder[i - 1] - 1]; nufork = curtree.nodep[spp + i - 2]; lastsp = (i == spp); addpreorder(curtree.root, item, nufork, true, true); promlk_add(there, item, nufork, false); like = prot_evaluate(curtree.root); rearrange(&curtree.root); if (curtree.likelihood > bestree.likelihood) { prot_copy_(&curtree, &bestree, nonodes, grcategs); } if (progress) { writename(i - 1, 1, enterorder); #ifdef WIN32 phyFillScreenColor(); #endif } if (lastsp && global) { if (progress) { printf("Doing global rearrangements\n"); printf(" !"); for (j = 1; j <= nonodes; j++) if ( j % (( nonodes / 72 ) + 1 ) == 0 ) putchar('-'); printf("!\n"); } bestlike = bestyet; do { if (progress) printf(" "); gotlike = bestlike; for (j = 0; j < nonodes; j++) { bestyet = UNDEFINED; item = curtree.nodep[j]; if (item != curtree.root) { nufork = curtree.nodep[curtree.nodep[j]->back->index - 1]; promlk_re_move(&item, &nufork, false); there = curtree.root; addpreorder(curtree.root, item, nufork, true, true); promlk_add(there, item, nufork, false); } if (progress) { if ( j % (( nonodes / 72 ) + 1 ) == 0 ) putchar('.'); fflush(stdout); } } if (progress) putchar('\n'); } while (bestlike < gotlike); } } if (njumble > 1 && lastsp) { for (i = 0; i < spp; i++ ) promlk_re_move(&curtree.nodep[i], &dummy, false); if (jumb == 1 || bestree2.likelihood < bestree.likelihood) prot_copy_(&bestree, &bestree2, nonodes, grcategs); } if (jumb == njumble) { if (njumble > 1) prot_copy_(&bestree2, &curtree, nonodes, grcategs); else prot_copy_(&bestree, &curtree, nonodes, grcategs); fprintf(outfile, "\n\n"); treevaluate(); curtree.likelihood = prot_evaluate(curtree.root); promlk_printree(); summarize(); if (trout) { col = 0; promlk_treeout(curtree.root); } } } else { openfile(&intree, INTREE, "input tree file", "r", progname, intreename); numtrees = countsemic(&intree); if(numtrees > MAXSHIMOTREES) shimotrees = MAXSHIMOTREES; else shimotrees = numtrees; if (numtrees > 2) initseed(&inseed, &inseed0, seed); l0gl = (double *) Malloc(shimotrees * sizeof(double)); l0gf = (double **) Malloc(shimotrees * sizeof(double *)); for (i=0; i < shimotrees; ++i) l0gf[i] = (double *)Malloc(endsite * sizeof(double)); if (treeprint) { fprintf(outfile, "User-defined tree"); if (numtrees > 1) putc('s', outfile); fprintf(outfile, ":\n\n"); } fprintf(outfile, "\n\n"); which = 1; while (which <= numtrees) { /* These initializations required each time through the loop since multiple trees require re-initialization */ dummy_haslengths = true; nextnode = 0; dummy_first = true; goteof = false; treeread(intree, &root, dummy_treenode, &goteof, &dummy_first, curtree.nodep, &nextnode, &dummy_haslengths, &grbg, initpromlnode,false,nonodes); nonodes = nextnode; root = curtree.nodep[root->index - 1]; curtree.root = root; if (lngths) tymetrav(curtree.root, &x); if (goteof && (which <= numtrees)) { /* if we hit the end of the file prematurely */ printf ("\n"); printf ("ERROR: trees missing at end of file.\n"); printf ("\tExpected number of trees:\t\t%ld\n", numtrees); printf ("\tNumber of trees actually in file:\t%ld.\n\n", which - 1); exxit(-1); } curtree.start = curtree.nodep[0]->back; treevaluate(); promlk_printree(); summarize(); if (trout) { col = 0; promlk_treeout(curtree.root); } if(which < numtrees){ prot_freex_notip(nonodes, curtree.nodep); gdispose(curtree.root, &grbg, curtree.nodep); } which++; } FClose(intree); if (!auto_ && numtrees > 1 && weightsum > 1 ) standev2(numtrees, maxwhich, 0, endsite, maxlogl, l0gl, l0gf, aliasweight, seed); } if (usertree) { free(l0gl); for (i=0; i < shimotrees; i++) free(l0gf[i]); free(l0gf); } prot_freetable(); if (jumb < njumble) return; free(contribution); free_all_protx(nonodes2, curtree.nodep); if (!usertree) { free_all_protx(nonodes2, bestree.nodep); if (njumble > 1) free_all_protx(nonodes2, bestree2.nodep); } if (progress) { printf("\n\nOutput written to file \"%s\"\n\n", outfilename); if (trout) printf("Tree also written onto file \"%s\"\n", outtreename); putchar('\n'); } free(root); } /* maketree */ void clean_up() { /* Free and/or close stuff */ long i; free (rrate); free (probcat); free (rate); /* Seems to require freeing every time... */ for (i = 0; i < spp; i++) { free (y[i]); } free (y); free (nayme); free (enterorder); free (category); free (weight); free (alias); free (ally); free (location); free (aliasweight); free (probmat); free (eigmat); if (! (njumble <= 1)) freetree2(bestree2.nodep, nonodes2); FClose(infile); FClose(outfile); FClose(outtree); #ifdef MAC fixmacfile(outfilename); fixmacfile(outtreename); #endif } /* clean_up */ int main(int argc, Char *argv[]) { /* Protein Maximum Likelihood with molecular clock */ #ifdef MAC argc = 1; /* macsetup("Promlk", ""); */ argv[0] = "Promlk"; #endif init(argc,argv); progname = argv[0]; openfile(&infile, INFILE, "input file", "r", argv[0], infilename); openfile(&outfile, OUTFILE, "output file", "w", argv[0], outfilename); ibmpc = IBMCRT; ansi = ANSICRT; datasets = 1; mulsets = false; firstset = true; doinit(); if (trout) openfile(&outtree,OUTTREE,"output tree file","w",argv[0],outtreename); if (ctgry) openfile(&catfile,CATFILE,"categories file","r",argv[0],catfilename); if (weights || justwts) openfile(&weightfile,WEIGHTFILE,"weights file","r",argv[0],weightfilename); for (ith = 1; ith <= datasets; ith++) { if (datasets > 1) { fprintf(outfile, "Data set # %ld:\n\n", ith); if (progress) printf("\nData set # %ld:\n", ith); } getinput(); if (ith == 1) firstset = false; for (jumb = 1; jumb <= njumble; jumb++){ max_num_sibs = 0; maketree(); } } clean_up(); printf("Done.\n\n"); #ifdef WIN32 phyRestoreConsoleAttributes(); #endif return 0; } /* Protein Maximum Likelihood with molecular clock */