Mac binaries

[jabaws.git] / website / archive / binaries / mac / src / muscle / msadistkimura.cpp

diff --git a/website/archive/binaries/mac/src/muscle/msadistkimura.cpp b/website/archive/binaries/mac/src/muscle/msadistkimura.cpp
new file mode 100644 (file)
index 0000000..6903c83
--- /dev/null
+++ b/website/archive/binaries/mac/src/muscle/msadistkimura.cpp
@@ -0,0 +1,88 @@
+#include "muscle.h"\r
+#include "msa.h"\r
+#include <math.h>\r
+\r
+// "Standard" NJ distance: the Kimura measure.\r
+// This is defined to be:\r
+//\r
+//             log_e(1 - p - p*p/5)\r
+//\r
+// where p is the fraction of residues that differ, i.e.:\r
+//\r
+//             p = (1 - fractional_conservation)\r
+//\r
+// This measure is infinite for p = 0.8541 and is considered\r
+// unreliable for p >= 0.75 (according to the ClustalW docs).\r
+// ClustalW uses a table lookup for values > 0.75.\r
+// The following table was copied from the ClustalW file dayhoff.h.\r
+\r
+static int dayhoff_pams[]={\r
+  195,   /* 75.0% observed d; 195 PAMs estimated = 195% estimated d */\r
+  196,   /* 75.1% observed d; 196 PAMs estimated */\r
+                  197,    198,    199,    200,    200,    201,    202,  203,    \r
+  204,    205,    206,    207,    208,    209,    209,    210,    211,  212,    \r
+  213,    214,    215,    216,    217,    218,    219,    220,    221,  222,    \r
+  223,    224,    226,    227,    228,    229,    230,    231,    232,  233,    \r
+  234,    236,    237,    238,    239,    240,    241,    243,    244,  245,    \r
+  246,    248,    249,    250,    /* 250 PAMs = 80.3% observed d */          \r
+                                  252,    253,    254,    255,    257,  258,    \r
+  260,    261,    262,    264,    265,    267,    268,    270,    271,  273,    \r
+  274,    276,    277,    279,    281,    282,    284,    285,    287,  289,    \r
+  291,    292,    294,    296,    298,    299,    301,    303,    305,  307,    \r
+  309,    311,    313,    315,    317,    319,    321,    323,    325,  328,    \r
+  330,    332,    335,    337,    339,    342,    344,    347,    349,  352,    \r
+  354,    357,    360,    362,    365,    368,    371,    374,    377,  380,    \r
+  383,    386,    389,    393,    396,    399,    403,    407,    410,  414,    \r
+  418,    422,    426,    430,    434,    438,    442,    447,    451,  456,    \r
+  461,    466,    471,    476,    482,    487,    493,    498,    504,  511,    \r
+  517,    524,    531,    538,    545,    553,    560,    569,    577,  586,    \r
+  595,    605,    615,    626,    637,    649,    661,    675,    688,  703,    \r
+  719,    736,    754,    775,    796,    819,    845,    874,    907,  945,\r
+         /* 92.9% observed; 945 PAMs */    \r
+  988    /* 93.0% observed; 988 PAMs */\r
+};\r
+static int iTableEntries = sizeof(dayhoff_pams)/sizeof(dayhoff_pams[0]);\r
+\r
+double KimuraDist(double dPctId)\r
+       {\r
+       double p = 1 - dPctId;\r
+// Typical case: use Kimura's empirical formula\r
+       if (p < 0.75)\r
+               return -log(1 - p - (p*p)/5);\r
+\r
+// Per ClustalW, return 10.0 for anything over 93%\r
+       if (p > 0.93)\r
+               return 10.0;\r
+\r
+// If p >= 0.75, use table lookup\r
+       assert(p <= 1 && p >= 0.75);\r
+// Thanks for Michael Hoel for pointing out a bug\r
+// in the table index calculation in versions <= 3.52.\r
+       int iTableIndex = (int) ((p - 0.75)*1000 + 0.5);\r
+       if (iTableIndex < 0 || iTableIndex >= iTableEntries)\r
+               Quit("Internal error in MSADistKimura::ComputeDist");\r
+\r
+       return dayhoff_pams[iTableIndex] / 100.0;\r
+       }\r
+\r
+//double MSADistKimura::ComputeDist(const MSA &msa, unsigned uSeqIndex1,\r
+//  unsigned uSeqIndex2)\r
+//     {\r
+//     double dPctId = msa.GetPctIdentityPair(uSeqIndex1, uSeqIndex2);\r
+//     return KimuraDist(dPctId);\r
+//     }\r
+\r
+double KimuraDistToPctId(double dKimuraDist)\r
+       {\r
+// Solve quadratic equation\r
+       const double a = 0.2;\r
+       const double b = 1;\r
+       const double c = 1.0 - exp(-dKimuraDist);\r
+       const double p = (-b + sqrt(b*b + 4*a*c))/(2*a);\r
+       return 1 - p;\r
+       }\r
+\r
+double PctIdToHeightKimura(double dPctId)\r
+       {\r
+       return KimuraDist(dPctId);\r
+       }\r