package beast.evolution.substitutionmodel;
import beast.core.Description;
import beast.evolution.datatype.Aminoacid;
import beast.evolution.datatype.DataType;
import beast.math.statistic.DiscreteStatistics;
@Description("WAG model of amino acid evolution by " +
"S. Whelan and N. Goldman. 2000. Bioinformatics ?.")
public class WAG extends EmpiricalSubstitutionModel {
@Override
double[][] getEmpiricalRates() {
double[][] rate = new double[20][20];
// Q matrix from Beast 1
rate[0][1] = 0.610810;
rate[0][2] = 0.569079;
rate[0][3] = 0.821500;
rate[0][4] = 1.141050;
rate[0][5] = 1.011980;
rate[0][6] = 1.756410;
rate[0][7] = 1.572160;
rate[0][8] = 0.354813;
rate[0][9] = 0.219023;
rate[0][10] = 0.443935;
rate[0][11] = 1.005440;
rate[0][12] = 0.989475;
rate[0][13] = 0.233492;
rate[0][14] = 1.594890;
rate[0][15] = 3.733380;
rate[0][16] = 2.349220;
rate[0][17] = 0.125227;
rate[0][18] = 0.268987;
rate[0][19] = 2.221870;
rate[1][2] = 0.711690;
rate[1][3] = 0.165074;
rate[1][4] = 0.585809;
rate[1][5] = 3.360330;
rate[1][6] = 0.488649;
rate[1][7] = 0.650469;
rate[1][8] = 2.362040;
rate[1][9] = 0.206722;
rate[1][10] = 0.551450;
rate[1][11] = 5.925170;
rate[1][12] = 0.758446;
rate[1][13] = 0.116821;
rate[1][14] = 0.753467;
rate[1][15] = 1.357640;
rate[1][16] = 0.613776;
rate[1][17] = 1.294610;
rate[1][18] = 0.423612;
rate[1][19] = 0.280336;
rate[2][3] = 6.013660;
rate[2][4] = 0.296524;
rate[2][5] = 1.716740;
rate[2][6] = 1.056790;
rate[2][7] = 1.253910;
rate[2][8] = 4.378930;
rate[2][9] = 0.615636;
rate[2][10] = 0.147156;
rate[2][11] = 3.334390;
rate[2][12] = 0.224747;
rate[2][13] = 0.110793;
rate[2][14] = 0.217538;
rate[2][15] = 4.394450;
rate[2][16] = 2.257930;
rate[2][17] = 0.078463;
rate[2][18] = 1.208560;
rate[2][19] = 0.221176;
rate[3][4] = 0.033379;
rate[3][5] = 0.691268;
rate[3][6] = 6.833400;
rate[3][7] = 0.961142;
rate[3][8] = 1.032910;
rate[3][9] = 0.043523;
rate[3][10] = 0.093930;
rate[3][11] = 0.533362;
rate[3][12] = 0.116813;
rate[3][13] = 0.052004;
rate[3][14] = 0.472601;
rate[3][15] = 1.192810;
rate[3][16] = 0.417372;
rate[3][17] = 0.146348;
rate[3][18] = 0.363243;
rate[3][19] = 0.169417;
rate[4][5] = 0.109261;
rate[4][6] = 0.023920;
rate[4][7] = 0.341086;
rate[4][8] = 0.275403;
rate[4][9] = 0.189890;
rate[4][10] = 0.428414;
rate[4][11] = 0.083649;
rate[4][12] = 0.437393;
rate[4][13] = 0.441300;
rate[4][14] = 0.122303;
rate[4][15] = 1.560590;
rate[4][16] = 0.570186;
rate[4][17] = 0.795736;
rate[4][18] = 0.604634;
rate[4][19] = 1.114570;
rate[5][6] = 6.048790;
rate[5][7] = 0.366510;
rate[5][8] = 4.749460;
rate[5][9] = 0.131046;
rate[5][10] = 0.964886;
rate[5][11] = 4.308310;
rate[5][12] = 1.705070;
rate[5][13] = 0.110744;
rate[5][14] = 1.036370;
rate[5][15] = 1.141210;
rate[5][16] = 0.954144;
rate[5][17] = 0.243615;
rate[5][18] = 0.252457;
rate[5][19] = 0.333890;
rate[6][7] = 0.630832;
rate[6][8] = 0.635025;
rate[6][9] = 0.141320;
rate[6][10] = 0.172579;
rate[6][11] = 2.867580;
rate[6][12] = 0.353912;
rate[6][13] = 0.092310;
rate[6][14] = 0.755791;
rate[6][15] = 0.782467;
rate[6][16] = 0.914814;
rate[6][17] = 0.172682;
rate[6][18] = 0.217549;
rate[6][19] = 0.655045;
rate[7][8] = 0.276379;
rate[7][9] = 0.034151;
rate[7][10] = 0.068651;
rate[7][11] = 0.415992;
rate[7][12] = 0.194220;
rate[7][13] = 0.055288;
rate[7][14] = 0.273149;
rate[7][15] = 1.486700;
rate[7][16] = 0.251477;
rate[7][17] = 0.374321;
rate[7][18] = 0.114187;
rate[7][19] = 0.209108;
rate[8][9] = 0.152215;
rate[8][10] = 0.555096;
rate[8][11] = 0.992083;
rate[8][12] = 0.450867;
rate[8][13] = 0.756080;
rate[8][14] = 0.771387;
rate[8][15] = 0.822459;
rate[8][16] = 0.525511;
rate[8][17] = 0.289998;
rate[8][18] = 4.290350;
rate[8][19] = 0.131869;
rate[9][10] = 3.517820;
rate[9][11] = 0.360574;
rate[9][12] = 4.714220;
rate[9][13] = 1.177640;
rate[9][14] = 0.111502;
rate[9][15] = 0.353443;
rate[9][16] = 1.615050;
rate[9][17] = 0.234326;
rate[9][18] = 0.468951;
rate[9][19] = 8.659740;
rate[10][11] = 0.287583;
rate[10][12] = 5.375250;
rate[10][13] = 2.348200;
rate[10][14] = 0.462018;
rate[10][15] = 0.382421;
rate[10][16] = 0.364222;
rate[10][17] = 0.740259;
rate[10][18] = 0.443205;
rate[10][19] = 1.997370;
rate[11][12] = 1.032220;
rate[11][13] = 0.098843;
rate[11][14] = 0.619503;
rate[11][15] = 1.073780;
rate[11][16] = 1.537920;
rate[11][17] = 0.152232;
rate[11][18] = 0.147411;
rate[11][19] = 0.342012;
rate[12][13] = 1.320870;
rate[12][14] = 0.194864;
rate[12][15] = 0.556353;
rate[12][16] = 1.681970;
rate[12][17] = 0.570369;
rate[12][18] = 0.473810;
rate[12][19] = 2.282020;
rate[13][14] = 0.179896;
rate[13][15] = 0.606814;
rate[13][16] = 0.191467;
rate[13][17] = 1.699780;
rate[13][18] = 7.154480;
rate[13][19] = 0.725096;
rate[14][15] = 1.786490;
rate[14][16] = 0.885349;
rate[14][17] = 0.156619;
rate[14][18] = 0.239607;
rate[14][19] = 0.351250;
rate[15][16] = 4.847130;
rate[15][17] = 0.578784;
rate[15][18] = 0.872519;
rate[15][19] = 0.258861;
rate[16][17] = 0.126678;
rate[16][18] = 0.325490;
rate[16][19] = 1.547670;
rate[17][18] = 2.763540;
rate[17][19] = 0.409817;
rate[18][19] = 0.347826;
// Current (May 2011) version from http://www.ebi.ac.uk/goldman/WAG/index.html
//
// 0.551571
// 0.509848 0.635346
// 0.738998 0.147304 5.429420
// 1.027040 0.528191 0.265256 0.0302949
// 0.908598 3.035500 1.543640 0.616783 0.0988179
// 1.582850 0.439157 0.947198 6.174160 0.021352 5.469470
// 1.416720 0.584665 1.125560 0.865584 0.306674 0.330052 0.567717
// 0.316954 2.137150 3.956290 0.930676 0.248972 4.294110 0.570025 0.249410
// 0.193335 0.186979 0.554236 0.039437 0.170135 0.113917 0.127395 0.0304501 0.138190
// 0.397915 0.497671 0.131528 0.0848047 0.384287 0.869489 0.154263 0.0613037 0.499462 3.170970
// 0.906265 5.351420 3.012010 0.479855 0.0740339 3.894900 2.584430 0.373558 0.890432 0.323832 0.257555
// 0.893496 0.683162 0.198221 0.103754 0.390482 1.545260 0.315124 0.174100 0.404141 4.257460 4.854020 0.934276
// 0.210494 0.102711 0.0961621 0.0467304 0.398020 0.0999208 0.0811339 0.049931 0.679371 1.059470 2.115170 0.088836 1.190630
// 1.438550 0.679489 0.195081 0.423984 0.109404 0.933372 0.682355 0.243570 0.696198 0.0999288 0.415844 0.556896 0.171329 0.161444
// 3.370790 1.224190 3.974230 1.071760 1.407660 1.028870 0.704939 1.341820 0.740169 0.319440 0.344739 0.967130 0.493905 0.545931 1.613280
// 2.121110 0.554413 2.030060 0.374866 0.512984 0.857928 0.822765 0.225833 0.473307 1.458160 0.326622 1.386980 1.516120 0.171903 0.795384 4.378020
// 0.113133 1.163920 0.0719167 0.129767 0.717070 0.215737 0.156557 0.336983 0.262569 0.212483 0.665309 0.137505 0.515706 1.529640 0.139405 0.523742 0.110864
// 0.240735 0.381533 1.086000 0.325711 0.543833 0.227710 0.196303 0.103604 3.873440 0.420170 0.398618 0.133264 0.428437 6.454280 0.216046 0.786993 0.291148 2.485390
// 2.006010 0.251849 0.196246 0.152335 1.002140 0.301281 0.588731 0.187247 0.118358 7.821300 1.800340 0.305434 2.058450 0.649892 0.314887 0.232739 1.388230 0.365369 0.314730
//
// 0.0866279 0.043972 0.0390894 0.0570451 0.0193078 0.0367281 0.0580589 0.0832518 0.0244313 0.048466 0.086209 0.0620286 0.0195027 0.0384319 0.0457631 0.0695179 0.0610127 0.0143859 0.0352742 0.0708956
//
//
// A R N D C Q E G H I L K M F P S T W Y V
// Ala Arg Asn Asp Cys Gln Glu Gly His Ile Leu Lys Met Phe Pro Ser Thr Trp Tyr Val
//
//
// Symmetrical part of the WAG rate matrix and aa frequencies,
// estimated from 3905 globular protein amino acid sequences forming 182
// protein families.
// The first part above indicates the symmetric 'exchangeability'
// parameters, where s_ij = s_ji. The s_ij above are not scaled, but the
// PAML package will perform this scaling.
// The second part gives the amino acid frequencies (pi_i)
// estimated from the 3905 sequences. The net replacement rate from i to
// j is Q_ij = s_ij*pi_j.
// Prepared by Simon Whelan and Nick Goldman, December 2000.
//
// Citation:
// Whelan, S. and N. Goldman. 2001. A general empirical model of
// protein evolution derived from multiple protein families using
// a maximum likelihood approach. Molecular Biology and
// Evolution 18:691-699.
/*
// Q matrix
rate[0][1] = 0.551571; rate[0][2] = 0.509848;
rate[0][3] = 0.738998; rate[0][4] = 1.027040;
rate[0][5] = 0.908598; rate[0][6] = 1.582850;
rate[0][7] = 1.416720; rate[0][8] = 0.316954;
rate[0][9] = 0.193335; rate[0][10] = 0.397915;
rate[0][11] = 0.906265; rate[0][12] = 0.893496;
rate[0][13] = 0.210494; rate[0][14] = 1.438550;
rate[0][15] = 3.370790; rate[0][16] = 2.121110;
rate[0][17] = 0.113133; rate[0][18] = 0.240735;
rate[0][19] = 2.006010;
rate[0][1] = 0.551571;rate[0][2] = 0.509848;
rate[0][3] = 0.738998;rate[0][4] = 1.027040;
rate[0][5] = 0.908598;rate[0][6] = 1.582850;
rate[0][7] = 1.416720;rate[0][8] = 0.316954;
rate[0][9] = 0.193335;rate[0][10] = 0.397915;
rate[0][11] = 0.906265;rate[0][12] = 0.893496;
rate[0][13] = 0.210494;rate[0][14] = 1.438550;
rate[0][15] = 3.370790;rate[0][16] = 2.121110;
rate[0][17] = 0.113133;rate[0][18] = 0.240735;
rate[0][19] = 2.006010;
rate[1][2] = 0.635346;rate[1][3] = 0.147304;
rate[1][4] = 0.528191;rate[1][5] = 3.0355;
rate[1][6] = 0.439157;rate[1][7] = 0.584665;
rate[1][8] = 2.13715;rate[1][9] = 0.186979;
rate[1][10] = 0.497671;rate[1][11] = 5.35142;
rate[1][12] = 0.683162;rate[1][13] = 0.102711;
rate[1][14] = 0.679489;rate[1][15] = 1.22419;
rate[1][16] = 0.554413;rate[1][17] = 1.16392;
rate[1][18] = 0.381533;rate[1][19] = 0.251849;
rate[2][3] = 5.42942;rate[2][4] = 0.265256;
rate[2][5] = 1.54364;rate[2][6] = 0.947198;
rate[2][7] = 1.12556;rate[2][8] = 3.95629;
rate[2][9] = 0.554236;rate[2][10] = 0.131528;
rate[2][11] = 3.01201;rate[2][12] = 0.198221;
rate[2][13] = 0.0961621;rate[2][14] = 0.195081;
rate[2][15] = 3.97423;rate[2][16] = 2.03006;
rate[2][17] = 0.0719167;rate[2][18] = 1.086;
rate[2][19] = 0.196246;
rate[3][4] = 0.0302949;rate[3][5] = 0.616783;
rate[3][6] = 6.17416;rate[3][7] = 0.865584;
rate[3][8] = 0.930676;rate[3][9] = 0.039437;
rate[3][10] = 0.0848047;rate[3][11] = 0.479855;
rate[3][12] = 0.103754;rate[3][13] = 0.0467304;
rate[3][14] = 0.423984;rate[3][15] = 1.07176;
rate[3][16] = 0.374866;rate[3][17] = 0.129767;
rate[3][18] = 0.325711;rate[3][19] = 0.152335;
rate[4][5] = 0.0988179;rate[4][6] = 0.021352;
rate[4][7] = 0.306674;rate[4][8] = 0.248972;
rate[4][9] = 0.170135;rate[4][10] = 0.384287;
rate[4][11] = 0.0740339;rate[4][12] = 0.390482;
rate[4][13] = 0.39802;rate[4][14] = 0.109404;
rate[4][15] = 1.40766;rate[4][16] = 0.512984;
rate[4][17] = 0.71707;rate[4][18] = 0.543833;
rate[4][19] = 1.00214;
rate[5][6] = 5.46947;rate[5][7] = 0.330052;
rate[5][8] = 4.29411;rate[5][9] = 0.113917;
rate[5][10] = 0.869489;rate[5][11] = 3.8949;
rate[5][12] = 1.54526;rate[5][13] = 0.0999208;
rate[5][14] = 0.933372;rate[5][15] = 1.02887;
rate[5][16] = 0.857928;rate[5][17] = 0.215737;
rate[5][18] = 0.22771;rate[5][19] = 0.301281;
rate[6][7] = 0.567717;rate[6][8] = 0.570025;
rate[6][9] = 0.127395;rate[6][10] = 0.154263;
rate[6][11] = 2.58443;rate[6][12] = 0.315124;
rate[6][13] = 0.0811339;rate[6][14] = 0.682355;
rate[6][15] = 0.704939;rate[6][16] = 0.822765;
rate[6][17] = 0.156557;rate[6][18] = 0.196303;
rate[6][19] = 0.588731;
rate[7][8] = 0.24941;rate[7][9] = 0.0304501;
rate[7][10] = 0.0613037;rate[7][11] = 0.373558;
rate[7][12] = 0.1741;rate[7][13] = 0.049931;
rate[7][14] = 0.24357;rate[7][15] = 1.34182;
rate[7][16] = 0.225833;rate[7][17] = 0.336983;
rate[7][18] = 0.103604;rate[7][19] = 0.187247;
rate[8][9] = 0.13819;rate[8][10] = 0.499462;
rate[8][11] = 0.890432;rate[8][12] = 0.404141;
rate[8][13] = 0.679371;rate[8][14] = 0.696198;
rate[8][15] = 0.740169;rate[8][16] = 0.473307;
rate[8][17] = 0.262569;rate[8][18] = 3.87344;
rate[8][19] = 0.118358;
rate[9][10] = 3.17097;rate[9][11] = 0.323832;
rate[9][12] = 4.25746;rate[9][13] = 1.05947;
rate[9][14] = 0.0999288;rate[9][15] = 0.31944;
rate[9][16] = 1.45816;rate[9][17] = 0.212483;
rate[9][18] = 0.42017;rate[9][19] = 7.8213;
rate[10][11] = 0.257555;rate[10][12] = 4.85402;
rate[10][13] = 2.11517;rate[10][14] = 0.415844;
rate[10][15] = 0.344739;rate[10][16] = 0.326622;
rate[10][17] = 0.665309;rate[10][18] = 0.398618;
rate[10][19] = 1.80034;
rate[11][12] = 0.934276;rate[11][13] = 0.088836;
rate[11][14] = 0.556896;rate[11][15] = 0.96713;
rate[11][16] = 1.38698;rate[11][17] = 0.137505;
rate[11][18] = 0.133264;rate[11][19] = 0.305434;
rate[12][13] = 1.19063;rate[12][14] = 0.171329;
rate[12][15] = 0.493905;rate[12][16] = 1.51612;
rate[12][17] = 0.515706;rate[12][18] = 0.428437;
rate[12][19] = 2.05845;
rate[13][14] = 0.161444;rate[13][15] = 0.545931;
rate[13][16] = 0.171903;rate[13][17] = 1.52964;
rate[13][18] = 6.45428;rate[13][19] = 0.649892;
rate[14][15] = 1.61328;rate[14][16] = 0.795384;
rate[14][17] = 0.139405;rate[14][18] = 0.216046;
rate[14][19] = 0.314887;
rate[15][16] = 4.37802;rate[15][17] = 0.523742;
rate[15][18] = 0.786993;rate[15][19] = 0.232739;
rate[16][17] = 0.110864;rate[16][18] = 0.291148;
rate[16][19] = 1.38823;
rate[17][18] = 2.48539;rate[17][19] = 0.365369;
rate[18][19] = 0.31473;
*/
//// WAG* matrix, a variant of the WAG matrix from the same paper
////
//// 0.589718
//// 0.514347 0.67416
//// 0.731152 0.159054 5.30821
//// 1.21324 0.568449 0.233527 0.0379056
//// 1.03344 3.02808 1.62299 0.657364 0.0999068
//// 1.55788 0.443685 1.00122 6.04299 0.0284956 5.6037
//// 1.41993 0.629768 1.12717 0.88357 0.312544 0.346823 0.588609
//// 0.317684 2.31211 3.9337 0.958529 0.341479 4.87366 0.599188 0.279542
//// 0.214596 0.187262 0.527321 0.0390513 0.198958 0.125999 0.124553 0.0310522 0.162975
//// 0.400822 0.51821 0.144354 0.0869637 0.451124 0.873266 0.154936 0.067443 0.508952 3.1554
//// 0.881639 5.74119 2.88102 0.480308 0.0719929 4.19125 2.45392 0.381514 0.854485 0.320597 0.255092
//// 0.887458 0.660816 0.198404 0.0992829 0.428648 1.64018 0.294481 0.184545 0.40117 3.94646 4.81956 0.877057
//// 0.213179 0.122792 0.0848492 0.0458258 0.485001 0.109241 0.0873936 0.0552962 0.631713 1.06458 2.10414 0.0832422 1.14516
//// 1.51861 0.711498 0.204905 0.444152 0.109081 0.913179 0.720567 0.254626 0.722123 0.111722 0.422851 0.588203 0.179858 0.165205
//// 3.52499 1.35611 3.90127 1.09965 1.35221 0.87908 0.822025 1.33618 0.876688 0.321774 0.351913 1.05314 0.554077 0.563999 1.54694
//// 2.24161 0.594177 2.06787 0.395176 0.522957 0.829315 0.889765 0.236489 0.54992 1.48876 0.351564 1.45173 1.56873 0.188237 0.802531 4.02507
//// 0.135395 1.24086 0.0746093 0.142159 0.728065 0.208163 0.176397 0.366467 0.261223 0.259584 0.706082 0.159261 0.565299 1.58681 0.135024 0.528249 0.118584
//// 0.270321 0.386714 1.05269 0.326191 0.481954 0.210494 0.209621 0.108982 4.31772 0.44009 0.427718 0.155623 0.437069 6.49269 0.212945 0.742154 0.286443 2.42261
//// 1.92496 0.282892 0.193323 0.155419 1.10899 0.32893 0.588443 0.190095 0.119749 7.48376 1.82105 0.300343 2.03324 0.653015 0.325745 0.23769 1.4088 0.396884 0.353358
////
//// 0.0866279 0.043972 0.0390894 0.0570451 0.0193078 0.0367281 0.0580589 0.0832518 0.0244313 0.048466 0.086209 0.0620286 0.0195027 0.0384319 0.0457631 0.0695179 0.0610127 0.0143859 0.0352742 0.0708956
////
////
//// A R N D C Q E G H I L K M F P S T W Y V
//// Ala Arg Asn Asp Cys Gln Glu Gly His Ile Leu Lys Met Phe Pro Ser Thr Trp Tyr Val
////
////
//// Symmetrical part of the WAG* rate matrix and aa frequencies,
//// estimated from 3905 globular protein amino acid sequences forming 182
//// protein families.
//// The first part above indicates the symmetric 'exchangeability'
//// parameters, where s_ij = s_ji. The s_ij above are not scaled, but the
//// PAML package will perform this scaling.
//// The second part gives the amino acid frequencies (pi_i)
//// estimated from the 3905 sequences. The net replacement rate from i to
//// j is Q_ij = s_ij*pi_j.
//// Prepared by Simon Whelan and Nick Goldman, December 2000.
////
//// Citation:
//// Whelan, S. and N. Goldman. 2001. A general empirical model of
//// protein evolution derived from multiple protein families using
//// a maximum likelihood approach. Molecular Biology and
//// Evolution 18:691-699.
//
// rate[0][1] = 0.589718;rate[0][2] = 0.514347;
// rate[0][3] = 0.731152;rate[0][4] = 1.21324;
// rate[0][5] = 1.03344;rate[0][6] = 1.55788;
// rate[0][7] = 1.41993;rate[0][8] = 0.317684;
// rate[0][9] = 0.214596;rate[0][10] = 0.400822;
// rate[0][11] = 0.881639;rate[0][12] = 0.887458;
// rate[0][13] = 0.213179;rate[0][14] = 1.51861;
// rate[0][15] = 3.52499;rate[0][16] = 2.24161;
// rate[0][17] = 0.135395;rate[0][18] = 0.270321;
// rate[0][19] = 1.92496;
//
// rate[1][2] = 0.67416;rate[1][3] = 0.159054;
// rate[1][4] = 0.568449;rate[1][5] = 3.02808;
// rate[1][6] = 0.443685;rate[1][7] = 0.629768;
// rate[1][8] = 2.31211;rate[1][9] = 0.187262;
// rate[1][10] = 0.51821;rate[1][11] = 5.74119;
// rate[1][12] = 0.660816;rate[1][13] = 0.122792;
// rate[1][14] = 0.711498;rate[1][15] = 1.35611;
// rate[1][16] = 0.594177;rate[1][17] = 1.24086;
// rate[1][18] = 0.386714;rate[1][19] = 0.282892;
//
//
// rate[2][3] = 5.30821;rate[2][4] = 0.233527;
// rate[2][5] = 1.62299;rate[2][6] = 1.00122;
// rate[2][7] = 1.12717;rate[2][8] = 3.9337;
// rate[2][9] = 0.527321;rate[2][10] = 0.144354;
// rate[2][11] = 2.88102;rate[2][12] = 0.198404;
// rate[2][13] = 0.0848492;rate[2][14] = 0.204905;
// rate[2][15] = 3.90127;rate[2][16] = 2.06787;
// rate[2][17] = 0.0746093;rate[2][18] = 1.05269;
// rate[2][19] = 0.193323;
//
// rate[3][4] = 0.0379056;rate[3][5] = 0.657364;
// rate[3][6] = 6.04299;rate[3][7] = 0.88357;
// rate[3][8] = 0.958529;rate[3][9] = 0.0390513;
// rate[3][10] = 0.0869637;rate[3][11] = 0.480308;
// rate[3][12] = 0.0992829;rate[3][13] = 0.0458258;
// rate[3][14] = 0.444152;rate[3][15] = 1.09965;
// rate[3][16] = 0.395176;rate[3][17] = 0.142159;
// rate[3][18] = 0.326191;rate[3][19] = 0.155419;
//
//
// rate[4][5] = 0.0999068;rate[4][6] = 0.0284956;
// rate[4][7] = 0.312544;rate[4][8] = 0.341479;
// rate[4][9] = 0.198958;rate[4][10] = 0.451124;
// rate[4][11] = 0.0719929;rate[4][12] = 0.428648;
// rate[4][13] = 0.485001;rate[4][14] = 0.109081;
// rate[4][15] = 1.35221;rate[4][16] = 0.522957;
// rate[4][17] = 0.728065;rate[4][18] = 0.481954;
// rate[4][19] = 1.10899;
//
// rate[5][6] = 5.6037;rate[5][7] = 0.346823;
// rate[5][8] = 4.87366;rate[5][9] = 0.125999;
// rate[5][10] = 0.873266;rate[5][11] = 4.19125;
// rate[5][12] = 1.64018;rate[5][13] = 0.109241;
// rate[5][14] = 0.913179;rate[5][15] = 0.87908;
// rate[5][16] = 0.829315;rate[5][17] = 0.208163;
// rate[5][18] = 0.210494;rate[5][19] = 0.32893;
//
//
// rate[6][7] = 0.588609;rate[6][8] = 0.599188;
// rate[6][9] = 0.124553;rate[6][10] = 0.154936;
// rate[6][11] = 2.45392;rate[6][12] = 0.294481;
// rate[6][13] = 0.0873936;rate[6][14] = 0.720567;
// rate[6][15] = 0.822025;rate[6][16] = 0.889765;
// rate[6][17] = 0.176397;rate[6][18] = 0.209621;
// rate[6][19] = 0.588443;
//
// rate[7][8] = 0.279542;rate[7][9] = 0.0310522;
// rate[7][10] = 0.067443;rate[7][11] = 0.381514;
// rate[7][12] = 0.184545;rate[7][13] = 0.0552962;
// rate[7][14] = 0.254626;rate[7][15] = 1.33618;
// rate[7][16] = 0.236489;rate[7][17] = 0.366467;
// rate[7][18] = 0.108982;rate[7][19] = 0.190095;
//
//
// rate[8][9] = 0.162975;rate[8][10] = 0.508952;
// rate[8][11] = 0.854485;rate[8][12] = 0.40117;
// rate[8][13] = 0.631713;rate[8][14] = 0.722123;
// rate[8][15] = 0.876688;rate[8][16] = 0.54992;
// rate[8][17] = 0.261223;rate[8][18] = 4.31772;
// rate[8][19] = 0.119749;
//
// rate[9][10] = 3.1554;rate[9][11] = 0.320597;
// rate[9][12] = 3.94646;rate[9][13] = 1.06458;
// rate[9][14] = 0.111722;rate[9][15] = 0.321774;
// rate[9][16] = 1.48876;rate[9][17] = 0.259584;
// rate[9][18] = 0.44009;rate[9][19] = 7.48376;
//
//
// rate[10][11] = 0.255092;rate[10][12] = 4.81956;
// rate[10][13] = 2.10414;rate[10][14] = 0.422851;
// rate[10][15] = 0.351913;rate[10][16] = 0.351564;
// rate[10][17] = 0.706082;rate[10][18] = 0.427718;
// rate[10][19] = 1.82105;
//
// rate[11][12] = 0.877057;rate[11][13] = 0.0832422;
// rate[11][14] = 0.588203;rate[11][15] = 1.05314;
// rate[11][16] = 1.45173;rate[11][17] = 0.159261;
// rate[11][18] = 0.155623;rate[11][19] = 0.300343;
//
//
// rate[12][13] = 1.14516;rate[12][14] = 0.179858;
// rate[12][15] = 0.554077;rate[12][16] = 1.56873;
// rate[12][17] = 0.565299;rate[12][18] = 0.437069;
// rate[12][19] = 2.03324;
//
// rate[13][14] = 0.165205;rate[13][15] = 0.563999;
// rate[13][16] = 0.188237;rate[13][17] = 1.58681;
// rate[13][18] = 6.49269;rate[13][19] = 0.653015;
//
//
// rate[14][15] = 1.54694;rate[14][16] = 0.802531;
// rate[14][17] = 0.135024;rate[14][18] = 0.212945;
// rate[14][19] = 0.325745;
//
// rate[15][16] = 4.02507;rate[15][17] = 0.528249;
// rate[15][18] = 0.742154;rate[15][19] = 0.23769;
//
//
// rate[16][17] = 0.118584;rate[16][18] = 0.286443;
// rate[16][19] = 1.4088;
//
// rate[17][18] = 2.42261;rate[17][19] = 0.396884;
//
//
// rate[18][19] = 0.353358;
//
return rate;
}
@Override
public double[] getEmpiricalFrequencies() {
double[] f = new double[20];
// frequencies from May 2011
f[0] = 0.0866279;
f[1] = 0.043972;
f[2] = 0.0390894;
f[3] = 0.0570451;
f[4] = 0.0193078;
f[5] = 0.0367281;
f[6] = 0.0580589;
f[7] = 0.0832518;
f[8] = 0.0244313;
f[9] = 0.048466;
f[10] = 0.086209;
f[11] = 0.0620286;
f[12] = 0.0195027;
f[13] = 0.0384319;
f[14] = 0.0457631;
f[15] = 0.0695179;
f[16] = 0.0610127;
f[17] = 0.0143859;
f[18] = 0.0352742;
f[19] = 0.0708956;
// frequencies from Beast 1
f[0] = 0.0866;
f[1] = 0.0440;
f[2] = 0.0391;
f[3] = 0.0570;
f[4] = 0.0193;
f[5] = 0.0367;
f[6] = 0.0581;
f[7] = 0.0833;
f[8] = 0.0244;
f[9] = 0.0485;
f[10] = 0.0862;
f[11] = 0.0620;
f[12] = 0.0195;
f[13] = 0.0384;
f[14] = 0.0458;
f[15] = 0.0695;
f[16] = 0.0610;
f[17] = 0.0144;
f[18] = 0.0353;
f[19] = 0.0709;
return f;
}
@Override
public int[] getEncodingOrder() {
Aminoacid dataType = new Aminoacid();
String codeMap = dataType.getCodeMap();
int[] codeMapNrs = new int[dataType.getStateCount()];
String encoding = "ARNDCQEGHILKMFPSTWYV";
for (int i = 0; i < dataType.getStateCount(); i++) {
codeMapNrs[i] = encoding.indexOf(codeMap.charAt(i));
}
return codeMapNrs;
}
@Override
public boolean canHandleDataType(DataType dataType) {
return dataType instanceof Aminoacid;
}
public static void main(String[] args) {
WAG wag = new WAG();
//String aminoAcids = "ARNDCQEGHILKMFPSTWYV";
boolean[] class1 = {false,true,false,false,true,true,true,false,false,true,true,
false,true,false,false,false,false,true,true,true};
int within1 = 0;
int within2 = 0;
int between = 0;
double[] w1 = new double[45];
double[] w2 = new double[45];
double[] b = new double[100];
double[][] rates = wag.getEmpiricalRates();
for (int i = 0; i < 20; i++) {
for (int j = i+1; j < 20; j++) {
if (class1[i]) {
if (class1[j]) {
w1[within1] = rates[i][j];
within1 += 1;
} else {
b[between] = rates[i][j];
between += 1;
}
} else {
if (!class1[j]) {
w2[within2] = rates[i][j];
within2 += 1;
} else {
b[between] = rates[i][j];
between += 1;
}
}
}
}
System.out.println("Within 1 mean rate = " + DiscreteStatistics.mean(w1));
System.out.println("Within 2 mean rate = " + DiscreteStatistics.mean(w2));
System.out.println("Between mean rate = " + DiscreteStatistics.mean(b));
System.out.println("Within 1 rate stdev = " + DiscreteStatistics.stdev(w1));
System.out.println("Within 2 rate stdev = " + DiscreteStatistics.stdev(w2));
System.out.println("Between rate stdev = " + DiscreteStatistics.stdev(b));
System.out.println("Within 1 rate stderr = " + DiscreteStatistics.stdev(w1)/Math.sqrt(within1));
System.out.println("Within 2 rate stderr = " + DiscreteStatistics.stdev(w2)/Math.sqrt(within2));
System.out.println("Between rate stderr = " + DiscreteStatistics.stdev(b)/Math.sqrt(between));
double sse = 0;
double meanb = DiscreteStatistics.mean(b);
for (int i = 0; i < b.length; i++) {
sse += (meanb - b[i]) * (meanb - b[i]);
}
sse /= b.length;
sse = Math.sqrt(sse);
System.out.println("Between rate sse = " + sse);
System.out.println("Within 1 count = " + within1);
System.out.println("Within 2 count = " + within2);
System.out.println("Between count = " + between);
}
} // class WAG