main.cpp

#include <iostream>
#include <vector>
#include <algorithm>
#include "newick.cpp"
#include "expectation_step.cpp"
#include "process_maf.cpp"
#include "EM_algorithm_complex.cpp"
#include "process_maf.h"
#include <gsl/gsl_eigen.h>
#include <stdlib.h>

double exchange[63*64/2] = { 0.581210,
15.563110, 0.478004,
0.632709, 21.431688, 0.529454,
0.749020, 0.077211, 0.197225, 0.159506,
0.095275, 1.032016, 0.010836, 0.389958, 15.753929,
0.134920, 0.099137, 1.110098, 0.146316, 49.318435, 21.193136,
0.093855, 0.266767, 0.000000, 0.915935, 15.814609, 32.905174, 15.876484,
2.687463, 0.062497, 0.556006, 0.147354, 0.612457, 0.000000, 0.000000, 0.021614,
0.409508, 5.014965, 0.399427, 1.883728, 0.428255, 3.005682, 0.000000, 1.242239, 0.517978,
1.739711, 0.418899, 5.309414, 0.122260, 0.072084, 0.161638, 1.648589, 0.000000, 40.487998, 1.275128,
0.257829, 0.843298, 0.212030, 3.434994, 0.000000, 0.875068, 0.313143, 2.308739, 0.331741, 46.143335, 0.863468,
0.308499, 0.034663, 0.222327, 0.000000, 1.647203, 0.530892, 0.391866, 0.522821, 0.277582, 0.000000, 0.306026, 0.000000,
0.008095, 0.697456, 0.022813, 0.000000, 0.191011, 1.071397, 0.000000, 0.411128, 0.000000, 0.549492, 0.082179, 0.111141, 11.541455,
0.104394, 0.083716, 0.272127, 0.088265, 0.432972, 0.000000, 2.943965, 0.000000, 0.054686, 0.145827, 0.492993, 0.139564, 3.147210, 0.459307,
0.047091, 0.000000, 0.000000, 0.464352, 0.007182, 0.380534, 0.247606, 0.766087, 0.011549, 0.100558, 0.000000, 0.370674, 11.003782, 27.760545, 0.568313,
0.910273, 0.334042, 0.301939, 0.133802, 0.187701, 0.000000, 0.114210, 0.134237, 0.000000, 0.153931, 0.814700, 0.110467, 0.060661, 0.000000, 0.099956, 0.073124,
0.262213, 1.650487, 0.352155, 0.268885, 0.080474, 0.184282, 0.000000, 0.053395, 0.000000, 0.846812, 0.655845, 0.000000, 0.000252, 0.000000, 0.059180, 0.236031, 0.855325,
0.486331, 0.000000, 1.122447, 0.487166, 0.117094, 0.260023, 0.338839, 0.000000, 0.698496, 0.226905, 0.000000, 0.181440, 0.120316, 0.178966, 0.134563, 0.000000, 32.534553, 2.495203,
0.286895, 0.148017, 0.255525, 1.196662, 0.045141, 0.052327, 0.120422, 0.130288, 0.000000, 0.000000, 0.000000, 0.482767, 0.000000, 0.209644, 0.031002, 0.000000, 1.046003, 58.410734, 0.961689,
0.066199, 0.016481, 0.044313, 0.000000, 0.728036, 0.072765, 0.000000, 0.024678, 0.046264, 0.007519, 0.063470, 0.000000, 0.000000, 0.103493, 0.077039, 0.081201, 0.738042, 0.028871, 0.191031, 0.000000,
0.040067, 0.206127, 0.015669, 0.000000, 0.180330, 0.583841, 0.077471, 0.279440, 0.000000, 0.000000, 0.068522, 0.490914, 0.601179, 0.462959, 0.000000, 0.000000, 0.067797, 1.040454, 0.031001, 0.251090, 20.159234,
0.052950, 0.000000, 0.060574, 0.061778, 0.000000, 0.000000, 2.194430, 0.178304, 0.077773, 0.000000, 0.000000, 0.011175, 0.000000, 0.000000, 0.138892, 0.125348, 0.147478, 0.025787, 1.829750, 0.075639, 46.887092, 39.633304,
0.022091, 0.000000, 0.033558, 0.132681, 0.097633, 0.262964, 0.051213, 0.455749, 0.000000, 0.408639, 0.000000, 0.058674, 0.556299, 0.000000, 0.000000, 0.152429, 0.000039, 0.230734, 0.095356, 0.741125, 22.763060, 49.492507, 18.490177,
0.944127, 0.000000, 0.843122, 0.552476, 0.047695, 1.070144, 0.000000, 0.000000, 17.187125, 0.757537, 1.015590, 0.000000, 0.258196, 0.000000, 0.105490, 0.085306, 7.314215, 0.000000, 0.000000, 0.732925, 1.291685, 0.623259, 0.000000, 0.579306,
0.229394, 1.915409, 0.789940, 0.000000, 0.135078, 0.727031, 0.003531, 0.000000, 5.794902, 3.482350, 1.218529, 0.310979, 0.039164, 0.000000, 0.053583, 0.441751, 0.591085, 14.418730, 0.859429, 0.466483, 0.057566, 2.640641, 0.002621, 0.037862, 20.618005,
1.517375, 0.797452, 0.000000, 0.000000, 0.000000, 0.000000, 0.752398, 1.551200, 2.248566, 0.000000, 44.061056, 0.899375, 0.152647, 0.451626, 0.029082, 0.000000, 0.000000, 0.051836, 17.017765, 0.514944, 0.000000, 0.000000, 7.166142, 0.431489, 150.361161, 55.270133,
0.009896, 0.000000, 0.000000, 0.639135, 0.000562, 0.000000, 0.000000, 0.141571, 4.984891, 0.000000, 5.396914, 0.829176, 0.000000, 0.149664, 0.004566, 0.000000, 0.103796, 0.000000, 0.066265, 4.116663, 0.000000, 0.000000, 0.000000, 0.690723, 20.304468, 98.868462, 26.296639,
0.076185, 0.000000, 0.000000, 0.286589, 0.224284, 0.060457, 0.021870, 0.000000, 0.039925, 0.000000, 0.000000, 0.162746, 1.568783, 0.417423, 0.739826, 0.137057, 0.410596, 0.000000, 0.518875, 0.234408, 1.110975, 0.000000, 0.000000, 0.439714, 1.081426, 0.806999, 0.523769, 0.000000,
0.000000, 0.000000, 0.067671, 0.249910, 0.555469, 0.000000, 0.000000, 0.231677, 0.000000, 0.525589, 0.297108, 0.091708, 1.736653, 4.899255, 0.412733, 0.000000, 0.000000, 1.510406, 0.717409, 0.428116, 1.475873, 2.332245, 0.000000, 0.000000, 0.000000, 4.634631, 0.442689, 0.004883, 12.733571,
0.059779, 0.000000, 0.042869, 0.533397, 0.183715, 0.072118, 0.274229, 0.000000, 0.000000, 0.000000, 0.161290, 0.291983, 1.034426, 0.722227, 1.564296, 0.000000, 0.000000, 0.000000, 2.225157, 0.160754, 0.523704, 0.000000, 0.837382, 1.178229, 0.000000, 2.758903, 2.687397, 0.000000, 34.614559, 15.129294,
0.049996, 0.317931, 0.000000, 0.000000, 0.356306, 0.000000, 0.486823, 0.000000, 0.000000, 0.306962, 0.000000, 0.300807, 1.827523, 0.000000, 0.213560, 2.509984, 0.250805, 0.405624, 0.000000, 0.713324, 0.385827, 0.000000, 2.337610, 0.290645, 0.690732, 0.000000, 0.058561, 1.456656, 11.227375, 39.457012, 12.904200,
0.987910, 0.207730, 0.000000, 0.271599, 0.277918, 0.164910, 0.067738, 0.000000, 0.114197, 0.213350, 0.000000, 0.249641, 0.124566, 0.146009, 0.012716, 0.000000, 1.249151, 0.462433, 0.000000, 0.000000, 0.119174, 0.083300, 0.032816, 0.000000, 0.449821, 0.291481, 0.000000, 0.052229, 0.018576, 0.061246, 0.030042, 0.000000,
0.155840, 3.158621, 0.140173, 0.000000, 0.009189, 0.518859, 0.113705, 0.000000, 0.000000, 1.806903, 0.010675, 0.000000, 0.007381, 0.081087, 0.009648, 0.000000, 0.404504, 1.063950, 0.000000, 0.000000, 0.000000, 0.539935, 0.305632, 0.000000, 0.216981, 0.000000, 0.000000, 0.116171, 0.216011, 0.600828, 0.788736, 0.000000, 1.292750,
0.000000, 0.336075, 1.441360, 0.303222, 0.162858, 0.000000, 0.448497, 0.228369, 0.000000, 0.374496, 0.641901, 0.133087, 0.109522, 0.000000, 0.103542, 0.208657, 0.000000, 0.000000, 3.747700, 0.547344, 0.050998, 0.000000, 0.252905, 0.118388, 0.000000, 0.000000, 0.892253, 0.153856, 0.111432, 0.000000, 0.019578, 0.040944, 20.688427, 2.273031,
0.119965, 0.000000, 0.102133, 2.261594, 0.056324, 0.000000, 0.037167, 0.311712, 0.000000, 0.000000, 0.000000, 0.853791, 0.002873, 0.007741, 0.024273, 0.049773, 0.104988, 0.000000, 0.425768, 0.487591, 0.035713, 0.000000, 0.000000, 0.184667, 0.000000, 0.407509, 0.206089, 0.000000, 0.000000, 0.000000, 0.000000, 0.263377, 1.133080, 23.651072, 1.384351,
0.221583, 0.106956, 0.100740, 0.083688, 4.006297, 0.407657, 0.000000, 0.449088, 0.194676, 0.096544, 0.060718, 0.000000, 0.521904, 0.027122, 0.325806, 0.000000, 0.218561, 0.000000, 0.108089, 0.000000, 0.903898, 0.324377, 0.000000, 0.318105, 0.000000, 0.000809, 0.000000, 0.049001, 0.171605, 1.054405, 0.023873, 0.337838, 0.624720, 0.018192, 0.224498, 0.038672,
0.000000, 0.107133, 0.132738, 0.220992, 0.090073, 3.924319, 0.178088, 0.000000, 0.098026, 1.165314, 0.000000, 0.195615, 0.014375, 0.031834, 0.018956, 0.423232, 0.000000, 0.239858, 0.148558, 0.000000, 0.298245, 0.307821, 0.000000, 0.560624, 0.000000, 0.043917, 1.143561, 0.042065, 0.660493, 0.000000, 1.175586, 1.225266, 0.098629, 0.572943, 0.000000, 0.055702, 16.830036,
0.125299, 0.090062, 0.290758, 0.100231, 0.000000, 0.475077, 8.873984, 0.389810, 0.006809, 0.000000, 0.591616, 0.298536, 0.449019, 0.000000, 0.814238, 0.060188, 0.090764, 0.015783, 0.288974, 0.000000, 0.000000, 0.294032, 3.462958, 0.167758, 0.000000, 0.000000, 1.087950, 0.072984, 0.034194, 0.000000, 0.335123, 0.838866, 0.129048, 0.102170, 1.427760, 0.097373, 54.926871, 23.740868,
0.079728, 0.193276, 0.015683, 0.183531, 0.106547, 0.000000, 0.065435, 2.317427, 0.000000, 0.245919, 0.176927, 0.809508, 0.000000, 0.465597, 0.023674, 0.160960, 0.076264, 0.000000, 0.000000, 0.106577, 0.040440, 0.540362, 0.465913, 0.338494, 0.481099, 0.223491, 0.000000, 0.030650, 0.000000, 2.319882, 0.000000, 0.000000, 0.000000, 0.069241, 0.132037, 0.388986, 16.496601, 33.395296, 14.732489,
0.480661, 0.303715, 0.083047, 0.196624, 0.488275, 0.000000, 0.000000, 0.824541, 1.550415, 0.000000, 0.000000, 0.000000, 0.000000, 0.355596, 0.082294, 0.000000, 0.283191, 0.096306, 0.000000, 0.005266, 0.000000, 0.376161, 0.000000, 0.407138, 1.999530, 0.482876, 0.000000, 0.000000, 0.031756, 0.138489, 0.003763, 0.000000, 1.945803, 0.005497, 0.623839, 0.180562, 1.791773, 0.626305, 0.000000, 0.000000,
0.017167, 0.000000, 0.190885, 1.152310, 0.012930, 0.864321, 0.100982, 0.000000, 0.421233, 7.881538, 0.000000, 0.559986, 0.037409, 0.000000, 0.019720, 0.216379, 0.000000, 1.000745, 0.016604, 0.000000, 0.290029, 0.000000, 0.000000, 0.517306, 0.652268, 0.410752, 0.000000, 0.482156, 0.008488, 0.546516, 0.315243, 0.000000, 0.000000, 1.764232, 0.458357, 0.610885, 0.219036, 2.282234, 0.347131, 0.000000, 21.233897,
0.095359, 0.448136, 0.390840, 0.000000, 0.000000, 1.688249, 0.985432, 0.000000, 0.000000, 0.000000, 3.532209, 0.000000, 0.000000, 0.000000, 0.153060, 0.285654, 0.000000, 0.000000, 0.431492, 0.183953, 0.000000, 0.770997, 0.225006, 0.208561, 0.000000, 0.020922, 3.887695, 0.000000, 0.014578, 0.000000, 0.029792, 0.130489, 0.310164, 0.222473, 2.579923, 0.000000, 0.006457, 0.000000, 3.987643, 0.765572, 58.548571, 19.497275,
0.027388, 0.310643, 0.000000, 0.095622, 0.060460, 0.000000, 0.006577, 0.203578, 0.000000, 0.000000, 0.240123, 2.403076, 0.188588, 0.071722, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.175975, 0.121440, 0.041025, 0.389425, 0.000000, 0.000000, 0.376751, 0.165425, 0.218142, 0.000000, 0.000000, 0.000000, 0.071363, 0.000000, 0.204550, 0.000000, 0.580710, 0.129068, 0.000000, 0.039864, 0.729593, 14.539872, 31.246112, 20.652534,
0.207410, 0.000000, 0.033833, 0.141152, 1.797225, 0.358117, 0.429069, 0.000000, 0.273201, 0.154393, 0.000000, 0.000000, 8.600236, 0.827320, 0.758821, 0.292609, 0.165898, 0.062135, 0.000000, 0.000000, 0.000000, 0.175716, 0.000000, 0.573440, 0.070858, 0.000000, 0.124088, 0.000000, 1.680901, 1.547508, 0.323719, 0.621556, 0.516150, 0.000000, 0.076676, 0.014377, 4.021073, 0.000000, 0.000000, 0.303957, 0.142697, 0.171286, 0.909555, 0.000000,
0.000000, 0.593534, 0.108717, 0.000000, 0.000000, 1.297648, 0.000000, 0.313229, 0.000000, 0.800027, 0.000000, 0.000000, 0.541970, 5.837924, 0.000000, 1.181059, 0.000000, 0.000000, 0.254773, 0.192650, 0.000000, 0.000000, 1.170849, 0.000000, 0.159523, 0.000000, 0.000000, 0.246581, 0.884189, 2.232994, 1.683233, 0.000000, 0.214697, 0.229985, 0.000000, 0.108241, 0.214938, 2.331261, 0.096801, 0.207652, 0.000000, 0.240941, 0.412811, 0.009712, 14.297408,
0.090486, 0.165946, 0.068338, 0.000000, 0.792196, 0.000000, 2.444380, 0.000000, 0.000000, 0.000000, 0.423463, 0.070000, 2.438833, 0.000000, 2.577749, 0.568614, 0.022860, 0.000000, 0.204827, 0.016758, 0.000000, 0.596930, 0.327958, 0.214064, 0.010888, 0.111986, 0.000000, 0.000000, 0.069180, 0.293073, 1.723185, 1.120367, 0.088987, 0.000000, 0.697477, 0.000000, 0.000000, 0.285306, 6.784114, 0.000000, 0.547546, 0.000000, 0.000000, 0.036953, 40.639291, 21.481324,
0.034290, 0.000000, 0.013382, 0.307393, 0.000000, 0.292293, 0.000000, 0.969271, 0.000000, 0.000000, 0.031346, 0.442067, 0.483278, 1.356857, 0.000000, 4.146653, 0.066357, 0.199970, 0.000000, 0.000000, 0.378033, 0.000000, 0.000000, 0.000000, 0.000000, 0.545784, 0.127831, 0.000000, 0.000000, 0.000000, 0.000000, 1.214617, 0.000000, 0.110261, 0.331540, 0.177800, 0.175719, 0.228170, 0.237924, 1.373525, 0.185112, 0.088704, 0.000000, 0.102844, 14.663966, 33.108782, 16.892095,
0.886628, 0.646899, 0.515932, 0.247176, 0.083186, 0.000000, 0.000000, 0.000000, 0.063306, 0.554063, 0.000000, 0.003250, 0.650818, 0.349594, 0.027112, 0.000000, 6.981227, 0.000000, 0.000000, 0.137770, 0.766745, 0.000000, 0.000000, 0.000000, 5.302647, 0.000000, 0.029797, 0.000000, 0.000000, 0.057130, 0.000000, 0.015138, 1.441951, 0.303182, 0.000000, 0.000000, 0.239458, 0.210965, 0.163587, 0.000000, 0.606056, 0.211722, 0.097119, 0.000000, 0.489901, 0.397413, 0.236327, 0.000000,
0.084336, 0.498741, 0.000000, 0.000000, 0.034532, 0.120610, 0.021158, 0.034849, 0.015844, 0.379120, 0.000000, 0.000000, 0.051024, 0.000000, 0.031396, 0.059221, 0.084293, 3.513888, 0.152216, 0.000000, 0.000000, 0.061043, 0.040112, 0.000000, 0.000000, 1.410304, 0.129854, 0.000000, 0.000000, 0.028949, 0.001850, 0.214896, 0.001822, 0.336081, 0.072220, 0.000000, 0.043623, 0.059941, 0.032237, 0.045128, 0.000000, 0.210611, 0.000000, 0.000000, 0.013692, 0.000000, 0.057039, 0.078064, 0.587740,
0.960921, 0.250810, 0.474287, 0.692960, 0.000000, 0.000000, 0.034693, 0.000000, 0.000000, 0.172164, 1.269599, 0.548721, 0.200360, 0.000000, 0.111537, 0.230672, 0.000000, 0.965801, 14.234842, 0.140443, 0.000000, 0.000000, 1.841766, 0.000000, 2.459291, 1.785422, 10.091001, 0.000000, 1.583833, 0.000000, 1.218662, 0.000000, 0.000000, 0.000000, 2.888922, 0.138947, 0.161138, 0.000000, 0.276271, 0.136482, 0.499417, 0.000000, 0.699368, 0.000000, 0.500526, 0.000000, 0.329900, 0.186302, 486.412169, 2.684468,
0.000000, 0.000000, 0.034368, 0.492827, 0.025758, 0.018903, 0.055829, 0.093677, 0.000000, 0.000000, 0.040647, 0.238821, 0.003243, 0.071543, 0.032339, 0.039721, 0.109106, 0.000000, 0.116193, 2.380201, 0.000000, 0.047225, 0.000000, 0.059261, 0.273288, 0.000000, 0.034783, 0.446094, 0.356646, 0.247371, 0.052196, 0.000000, 0.041847, 0.000000, 0.002442, 0.244748, 0.027902, 0.066161, 0.058467, 0.040985, 0.044615, 0.000000, 0.045401, 0.060115, 0.043970, 0.103019, 0.024720, 0.015502, 2.245816, 37.140983, 1.783312,
0.174660, 0.105661, 0.120582, 0.150751, 2.654171, 0.316747, 0.239080, 0.220325, 0.030904, 0.557200, 0.000000, 0.226196, 0.000000, 0.116552, 0.162299, 0.000000, 0.245453, 0.088307, 0.191545, 0.031558, 3.285822, 0.000000, 0.000000, 0.055485, 1.232935, 0.297499, 0.000000, 0.000000, 0.724705, 0.488099, 0.381352, 0.647326, 0.138114, 0.000000, 0.102226, 0.071339, 2.546651, 0.148625, 0.456219, 0.299769, 0.000000, 0.272685, 0.000000, 0.158932, 0.120549, 0.232446, 0.000000, 0.239668, 2.633689, 0.000000, 0.838252, 0.102140,
0.076290, 0.222235, 0.091728, 0.126827, 0.435155, 2.568649, 0.112151, 0.216897, 0.000000, 0.922295, 0.247360, 1.059574, 0.196550, 0.000000, 0.027697, 0.204131, 0.000000, 0.065786, 0.279558, 0.234001, 0.143670, 6.273089, 0.000000, 0.000000, 0.000000, 0.180281, 0.000000, 0.143620, 0.000000, 4.435553, 0.000000, 0.000000, 0.075017, 0.480922, 0.000000, 0.000000, 0.503215, 4.046132, 0.234913, 0.000000, 0.396891, 0.000000, 0.387623, 0.194181, 0.527405, 0.000000, 0.039654, 0.220256, 0.229124, 0.420790, 0.000000, 0.014856, 12.802339,
0.112591, 0.170737, 0.152960, 0.118820, 0.409862, 0.023975, 4.040695, 0.434860, 0.000000, 0.137047, 0.210743, 0.563215, 0.000000, 0.000000, 0.280330, 0.111136, 0.076451, 0.042725, 0.328045, 0.139056, 0.000000, 0.203348, 8.116849, 0.000000, 0.000000, 0.000000, 2.227129, 0.000000, 0.000000, 0.728749, 0.000000, 0.844170, 0.087893, 0.165974, 0.176141, 0.000000, 0.654636, 0.540798, 4.066529, 0.016221, 0.000000, 0.000000, 0.240944, 0.206063, 0.000000, 0.306564, 0.039927, 0.134277, 0.412406, 0.104584, 3.983463, 0.000000, 42.416576, 19.141427,
0.075938, 0.161257, 0.080334, 0.179790, 0.221106, 0.308403, 0.449379, 1.633883, 0.032948, 1.017172, 0.000000, 0.879498, 0.304811, 0.212739, 0.000000, 0.000000, 0.088776, 0.352034, 0.000000, 0.200678, 0.000000, 0.000000, 0.930661, 3.330835, 0.000000, 0.284515, 0.000000, 0.186103, 0.478359, 0.000000, 1.095536, 0.593863, 0.016012, 0.000000, 0.106686, 0.234584, 0.486522, 0.000000, 0.586732, 2.194112, 0.208192, 0.648705, 0.175287, 0.000000, 0.067774, 0.263861, 0.353851, 0.000000, 0.032511, 0.005846, 0.303390, 0.421322, 13.025586, 29.255618, 13.712447,
0.423497, 0.150552, 0.000000, 0.361415, 0.000000, 0.000000, 0.000000, 0.000000, 1.671338, 0.662016, 1.770083, 0.747868, 0.478512, 0.000000, 0.121256, 0.159343, 0.000000, 0.006066, 4.940797, 0.000000, 0.594209, 0.000000, 1.341997, 0.046355, 13.055068, 2.249021, 2.828416, 1.182330, 1.540525, 0.427894, 0.759690, 0.000000, 0.216479, 0.000000, 0.982736, 0.000000, 0.452749, 0.000000, 0.191647, 0.028696, 3.153440, 0.002316, 1.034944, 0.120343, 0.247372, 0.000000, 0.000000, 0.353709, 423.207339, 0.558197, 209.336485, 0.000000, 3.322661, 0.224165, 1.624617, 0.668449,
0.167846, 0.593887, 0.000000, 0.000000, 0.033474, 0.379160, 0.111165, 0.000000, 0.000000, 2.734169, 0.414991, 0.264326, 0.081819, 0.626054, 0.040304, 0.000000, 0.046298, 0.000000, 0.221687, 0.810484, 0.000000, 0.194011, 0.218109, 0.082335, 0.193604, 7.605495, 0.000000, 0.000000, 0.000000, 0.025548, 0.000000, 0.587455, 0.000000, 0.433385, 0.889536, 0.000000, 0.198213, 0.440281, 0.046882, 0.086429, 0.000000, 1.155222, 0.348591, 0.065344, 0.000000, 0.817450, 0.208139, 0.000000, 0.000000, 2.954758, 1.420488, 0.000000, 0.202105, 1.914725, 0.429550, 0.273775, 6.909055,
0.015444, 0.000965, 0.029316, 0.002769, 0.026030, 0.017431, 0.058191, 0.004839, 0.000000, 0.041495, 0.216275, 0.016229, 0.048090, 0.021558, 0.059608, 0.015619, 0.004252, 0.085617, 0.099629, 0.025858, 0.017249, 0.008373, 0.072950, 0.000000, 0.118416, 0.050120, 3.265562, 0.000000, 0.043810, 0.081384, 0.135791, 0.003276, 0.000900, 0.007199, 0.027115, 0.008059, 0.022075, 0.003495, 0.073267, 0.000000, 0.046231, 0.011252, 0.409717, 0.000000, 0.041895, 0.004284, 0.048473, 0.009314, 0.000000, 0.180511, 3.575284, 0.193517, 0.058504, 0.000000, 0.423116, 0.000000, 7.963153, 0.516820,
0.000000, 0.000000, 0.077798, 0.272073, 0.058593, 0.001869, 0.085667, 0.175545, 0.000000, 0.000000, 0.000000, 1.126982, 0.000000, 0.000000, 0.045010, 0.251098, 0.090385, 0.658032, 0.000000, 0.000000, 0.025499, 0.066118, 0.000000, 0.121227, 0.063344, 0.000000, 0.361418, 1.744951, 0.161152, 0.545267, 0.188062, 0.086206, 0.219463, 0.000000, 0.000000, 0.172792, 0.110211, 0.116180, 0.165967, 0.286618, 0.185170, 0.184486, 0.000000, 0.319763, 0.170665, 0.000000, 0.034537, 0.363201, 0.330882, 0.000000, 0.000000, 1.537824, 0.079165, 0.126285, 0.054285, 0.958372, 4.001877, 98.188238, 0.241872,
0.057473, 0.262884, 0.005060, 0.000000, 0.085291, 0.031113, 0.000000, 0.130545, 0.039827, 0.124990, 0.051355, 0.000000, 1.143533, 0.000000, 0.560616, 0.438764, 0.181176, 0.000000, 0.000000, 0.021100, 0.000000, 1.297757, 0.851549, 0.000000, 0.301899, 0.000000, 0.000000, 0.346638, 29.700775, 7.007735, 10.452363, 5.270753, 0.000000, 0.000000, 0.108007, 0.235866, 0.079923, 0.000000, 0.000000, 0.327583, 0.059397, 0.000000, 0.033672, 0.046509, 0.978196, 0.000000, 0.175446, 0.511700, 1.202017, 0.000000, 0.905227, 0.218817, 1.207795, 1.004379, 0.139326, 0.000000, 1.949425, 0.162242, 0.064332, 0.107138,
0.000000, 0.000000, 0.077123, 0.105725, 0.000000, 0.228748, 0.004437, 0.000000, 0.051380, 0.000000, 0.000000, 0.046525, 0.000000, 0.693290, 0.140976, 0.000000, 0.165638, 0.000000, 0.000000, 0.213329, 0.000000, 0.448879, 0.179594, 0.037935, 0.213165, 0.000000, 0.000000, 0.159334, 1.148017, 5.554200, 0.667162, 0.000000, 0.000000, 0.000000, 0.020986, 0.103900, 0.000000, 0.528857, 0.057361, 0.000000, 0.028700, 0.072115, 0.044128, 0.051650, 0.000000, 1.043061, 0.027150, 0.000000, 1.266404, 1.181732, 0.000000, 0.972926, 0.000000, 1.059611, 0.000000, 0.050903, 0.000000, 0.436223, 0.153078, 0.183917, 0.677572,
0.000000, 0.136851, 0.071342, 0.000000, 0.000000, 0.110783, 0.136471, 0.081232, 0.033349, 0.130757, 0.000000, 0.000000, 0.528313, 0.070096, 1.137683, 0.070841, 0.085511, 0.266065, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.066063, 0.072512, 0.116559, 0.229946, 0.000000, 6.158381, 5.774936, 25.210182, 7.048552, 0.043426, 0.057664, 0.000000, 0.000000, 0.000000, 0.000000, 0.056143, 0.163533, 0.020904, 0.076903, 0.051841, 0.000000, 0.000000, 0.000000, 0.826724, 0.231042, 0.977674, 0.367270, 0.260490, 0.000000, 0.000000, 0.127730, 1.868945, 0.022373, 0.159692, 0.269767, 0.254645, 0.034510, 19.970821, 0.000000,
0.052141, 0.102223, 0.000000, 0.000000, 0.001868, 0.000000, 0.010455, 0.155385, 0.016960, 0.075331, 0.061180, 0.000000, 0.155947, 0.000000, 0.174958, 0.499488, 0.000000, 0.299579, 0.335970, 0.000000, 0.037930, 0.000000, 0.000000, 0.296135, 0.000000, 0.487602, 0.211727, 0.000000, 0.000000, 0.000000, 0.000000, 3.126326, 0.004730, 0.166950, 0.000000, 0.000000, 0.030902, 0.000000, 0.000000, 0.237758, 0.033598, 0.062251, 0.024403, 0.052678, 0.122858, 0.000000, 0.044946, 0.582508, 0.000000, 1.033681, 1.250157, 1.407673, 0.000000, 0.000000, 0.036481, 0.736912, 0.997027, 0.382341, 0.145731, 0.278208, 0.221501, 27.533240, 1.044282 };

int main() {

    auto start_time = std::chrono::steady_clock::now();

    double coding_freq[64] = {0.042558, 0.024513, 0.030138, 0.036293, 0.017995, 0.012479, 0.008171, 0.020108, 0.020893, 0.010055, 0.009439, 0.014649, 0.018597, 0.016947, 0.020770, 0.030135, 0.026648, 0.007776, 0.012271, 0.013887, 0.017633, 0.006954, 0.005409, 0.013514, 0.003213, 0.002707, 0.001910, 0.006327, 0.013506, 0.005776, 0.010752, 0.012711, 0.045143, 0.020107, 0.019335, 0.037538, 0.016269, 0.012098, 0.006217, 0.020106, 0.011230, 0.009807, 0.006088, 0.022349, 0.012260, 0.011228, 0.010804, 0.021456, 0.001050, 0.014490, 0.000519, 0.019166, 0.019118, 0.014150, 0.008827, 0.023470, 0.000684, 0.005011, 0.010397, 0.008147, 0.026476, 0.018304, 0.026556, 0.026866 };

    std::vector<CDS_info> all_CDS_info = get_CDS_info("/Users/sukhwanpark/Downloads/sacCer3.ncbiRefSeq.gtf");

    gsl_matrix *qmatrix = gsl_matrix_alloc(64, 64);
    gsl_vector *x = gsl_vector_alloc(2016);//freed

    for (int i = 0; i < 2016 ; i++) {
        if (exchange[i] != 0) {
            gsl_vector_set(x, i, exchange[i]);
        } else {
            gsl_vector_set(x, i, 0.001);
        }
    }

    for (size_t row = 1; row < 64; row++) {
        for (size_t col = 0; col < 64; col++) {
            if (row > col) {
                gsl_matrix_set(qmatrix, row, col, gsl_vector_get(x, (row - 1) * (row) / 2 + col));
                gsl_matrix_set(qmatrix, col, row, gsl_matrix_get(qmatrix, row, col));
                gsl_matrix_set(qmatrix, row, col, gsl_matrix_get(qmatrix, row, col) * coding_freq[col]);
                gsl_matrix_set(qmatrix, col, row, gsl_matrix_get(qmatrix, col, row) * coding_freq[row]);
            }
        }
    }
    double sum = 0.0;
    for (int row = 0; row < 64; row++) {
        sum = 0.0;
        for (int col = 0; col < 64; col++) {
            if (row != col) {
                sum += gsl_matrix_get(qmatrix, row, col);
            }
        }
        gsl_matrix_set(qmatrix, row, row, -sum);
    }
    sum = 0.0;
    for (size_t dia = 0; dia < 64; dia++) {
        sum -= gsl_matrix_get(qmatrix, dia, dia) * coding_freq[dia];
    }
    gsl_matrix_scale(qmatrix, 1 / sum);

    newick_start* start_point = new newick_start;
    newick_graph* end_point = new newick_graph;
    int species_num;
    int order_max_value = 0;
    int *newick_order_max;
    newick_order_max = &order_max_value;

    process_newick("/Users/sukhwanpark/Downloads/7yeast.nh", start_point, end_point, species_num, newick_order_max);

    std::vector<std::vector<aligned_codon>> aligned_codon_set = make_msa_to_aligned_coding_codon("/Users/sukhwanpark/Downloads/yeast_maf/yeast_CDS_ss.maf", all_CDS_info, species_num);

    auto middle_time = std::chrono::steady_clock::now();
    std::cout << std::chrono::duration_cast<std::chrono::milliseconds>(middle_time - start_time).count() / 1000.0 << std::endl;

    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
    std::shuffle(aligned_codon_set.begin(), aligned_codon_set.end(), std::default_random_engine(seed));

    //test newick
    /*int order = *newick_order_max;
    std::vector<newick_graph*> iterator = start_point->next;
    size_t size;
    while (order >= 0) {
        size = iterator.size();
        for (int num = 0; num < size; num++) {
            if (iterator[num]->order == order) {
                // Todo: Do something
                if (std::find(iterator.begin(), iterator.end(), iterator[num]->next) == iterator.end()) {
                    if (iterator[num]->next != NULL) {
                        iterator.emplace_back(iterator[num]->next);
                    }
                }
                iterator.erase(iterator.begin() + num);
                num--;
                size--;
            }
        }
        order--;
    }*/

    gsl_matrix *eigenvector = gsl_matrix_alloc(64, 64);//freed
    gsl_matrix *eigenvec_inverse = gsl_matrix_alloc(64, 64);//freed
    gsl_vector *eigenvalue = gsl_vector_alloc(64);//freed
    gsl_vector_complex *eigenvalue_com = gsl_vector_complex_alloc(64);//freed
    gsl_matrix_complex *eigenvector_com = gsl_matrix_complex_alloc(64, 64);//freed
    gsl_matrix_complex *eigenvec_inverse_com = gsl_matrix_complex_alloc(64, 64);//freed

    gsl_matrix *total_count = gsl_matrix_alloc(64, 64);//freed
    gsl_matrix *eigencount = gsl_matrix_alloc(64, 64);//freed
    gsl_matrix *eigencount_col = gsl_matrix_alloc(64, 64);//freed
    gsl_matrix_complex *eigencount_com = gsl_matrix_complex_alloc(64, 64);//freed
    gsl_matrix_complex *eigencount_col_com = gsl_matrix_complex_alloc(64, 64);//freed

    normalizing_qmatrix(qmatrix, coding_freq, qmatrix);

    double omega[64];

    double likelihood = 0.0;
    double likelihood_old = 1.0;
    int count = 0;
    bool eigen_is_complex = false;
    gsl_matrix *qmatrix_temp = gsl_matrix_alloc(64, 64);//freed

    while (abs(likelihood - likelihood_old) > 1.0e-8) {
        eigen_is_complex = false;
        likelihood_old = likelihood;
        likelihood = 0.0;
        int order = *newick_order_max;
        std::vector<newick_graph*> iterator = start_point->next;
        size_t size;
        while (order >= 0) {
            size = iterator.size();
            for (int num = 0; num < size; num++) {
                if (iterator[num]->order == order) {
                    iterator[num]->zero_entries.clear();
                    if (std::find(iterator.begin(), iterator.end(), iterator[num]->next) == iterator.end()) {
                        if (iterator[num]->next != NULL) {
                            iterator.emplace_back(iterator[num]->next);
                        }
                    }
                    iterator.erase(iterator.begin() + num);
                    num--;
                    size--;
                }
            }
            order--;
        }
        std::cout << "Iteration start : " << count << std::endl;

        gsl_matrix_memcpy(qmatrix_temp, qmatrix);

        get_eigenvector_and_inverse(qmatrix_temp, eigenvector, eigenvec_inverse, eigenvalue, eigenvector_com, eigenvalue_com, eigen_is_complex);

        if (eigen_is_complex == false) {
            set_matrices(start_point, eigenvector, eigenvec_inverse, eigenvalue, newick_order_max);
        } else {
            get_inverse_of_eigenvec_com(eigenvector_com, eigenvec_inverse_com);

            for (int row = 0; row < 64; row++) {
                for (int col = 0; col < 64; col++) {
                    gsl_matrix_complex_set(eigenvector_com, row, col,
                                           gsl_complex_mul(gsl_matrix_complex_get(eigenvector_com, row, col), gsl_vector_complex_get(eigenvalue_com, col)));
                }
            }

            set_matrices_complex(start_point, eigenvector_com, eigenvec_inverse_com, eigenvalue_com, newick_order_max);
        }

        ///////////////////////////
        //start of expectation step
        gsl_matrix_set_all(eigencount, 0.0);
        gsl_complex default_element = gsl_complex_rect(0.0, 0.0);
        gsl_matrix_complex_set_all(eigencount_com, default_element);

        double element, denominator, sum;

        //Todo : put the amount of the aligned codons to train
        for (int num_codon = 0; num_codon < 20000; num_codon++) {
            std::vector<newick_graph *> next_iterator = start_point->next;
            for (size_t num_t = 0; num_t < next_iterator.size(); num_t++) {
                for (char base = 0; base < 64; base++) {
                    next_iterator[num_t]->base[base] = false;
                }
            }

            denominator = 0.0;

            denominator = conduct_felsenstein(start_point, aligned_codon_set[num_codon], eigenvec_inverse, eigenvec_inverse_com, newick_order_max, end_point, coding_freq, eigen_is_complex);
            if (denominator == 0) {
                std::cout << "denominator is 0" << std::endl;
                exit(-1);
            }
            likelihood += denominator;
            denominator = 1 / denominator;

            calculate_upper(end_point, coding_freq, eigenvector, eigenvector_com, eigen_is_complex);

            int order = *newick_order_max;
            std::vector<newick_graph*> iterator = start_point->next;
            size_t size;
            gsl_matrix_set_all(eigencount_col, 0.0);
            gsl_matrix_complex_set_all(eigencount_col_com, default_element);
            if (eigen_is_complex == false) {
                while (order >= 0) {
                    size = iterator.size();
                    for (int num = 0; num < size; num++) {
                        if (iterator[num]->order == order) {
                            for (int row = 0; row < 64; row++) {
                                for (int col = 0; col < 64; col++) {
                                    gsl_matrix_set(eigencount_col, row, col, gsl_matrix_get(eigencount_col, row, col) + next_iterator[num]->down_projection[row] *
                                                                             next_iterator[num]->up_projection[col] * gsl_matrix_get(next_iterator[num]->expon_eigen, row, col));
                                }
                            }
                            if (std::find(iterator.begin(), iterator.end(), iterator[num]->next) == iterator.end()) {
                                if (iterator[num]->next != NULL) {
                                    iterator.emplace_back(iterator[num]->next);
                                }
                            }
                            iterator.erase(iterator.begin() + num);
                            num--;
                            size--;
                        }
                    }
                    order--;
                }
            } else {
                while (order >= 0) {
                    size = iterator.size();
                    for (int num = 0; num < size; num++) {
                        if (iterator[num]->order == order) {
                            for (int row = 0; row < 64; row++) {
                                for (int col = 0; col < 64; col++) {
                                    gsl_matrix_complex_set(eigencount_col_com, row, col, gsl_complex_add(
                                            gsl_matrix_complex_get(eigencount_col_com, row, col),
                                            gsl_complex_mul(gsl_complex_mul(iterator[num]->down_projection_com[row],gsl_matrix_complex_get(iterator[num]->expon_eigen_com, row, col)),
                                                            iterator[num]->up_projection_com[col])));
                                    gsl_complex subtract = gsl_complex_rect(0.0, 0.0);
                                    for (int zero = 0; zero < iterator[num]->zero_entries.size(); zero++) {
                                        int b = iterator[num]->zero_entries[zero] % 64;
                                        int a = (iterator[num]->zero_entries[zero] - b) / 64;
                                        gsl_complex_add(subtract, gsl_complex_mul_real(gsl_complex_mul(gsl_complex_mul(gsl_matrix_complex_get(iterator[num]->expon_eigen_com, row, col),
                                                                                                                       gsl_matrix_complex_get(eigenvec_inverse_com, col, b)),
                                                                                                       gsl_matrix_complex_get(eigenvector_com, a, row)),
                                                                                       iterator[num]->felsenstein[b] * iterator[num]->upper[a]));
                                    }
                                    gsl_matrix_complex_set(eigencount_col_com, row, col, gsl_complex_sub(gsl_matrix_complex_get(eigencount_col_com, row, col), subtract));
                                }
                            }
                            if (std::find(iterator.begin(), iterator.end(), iterator[num]->next) == iterator.end()) {
                                if (iterator[num]->next != NULL) {
                                    iterator.emplace_back(iterator[num]->next);
                                }
                            }
                            iterator.erase(iterator.begin() + num);
                            num--;
                            size--;
                        }
                    }
                    order--;
                }
            }
            if (eigen_is_complex == false) {
                gsl_matrix_scale(eigencount_col, denominator);
                gsl_matrix_add(eigencount, eigencount_col);
            } else {
                gsl_complex denominator_to_complex = gsl_complex_rect(denominator, 0.0);
                gsl_matrix_complex_scale(eigencount_col_com, denominator_to_complex);
                gsl_matrix_complex_add(eigencount_com, eigencount_col_com);
            }
        }

        for (int row = 0; row < 64; row++) {
            for (int col = 0; col < 64; col++) {
                if (eigen_is_complex == false) {
                    sum = 0.0;
                    for (int k = 0; k < 64; k++) {
                        element = 0.0;
                        for (int l = 0; l < 64; l++) {
                            element += gsl_matrix_get(eigenvector, col, l) * gsl_matrix_get(eigencount, k, l);
                            //sum += (gsl_matrix_get(eigenvector, col, l) * gsl_matrix_get(eigenvec_inverse, k, row) +
                                    //gsl_matrix_get(eigenvector, row, l) * gsl_matrix_get(eigenvec_inverse, k, col)) *
                                    //gsl_matrix_get(eigencount, k, l);
                        }
                        sum += gsl_matrix_get(eigenvec_inverse, k, row) * element;
                    }
                    if (((row != col && gsl_matrix_get(qmatrix, row, col) * sum >= 0.0) ||
                        (row == col && gsl_matrix_get(qmatrix, row, col) * sum <= 0)) /*&& row >= col*/) {
                        gsl_matrix_set(total_count, row, col, gsl_matrix_get(qmatrix, row, col) * sum);
                        //gsl_matrix_set(total_count, row, col, gsl_matrix_get(qmatrix, row, col) * sum / 2);
                    //} else if (row < col) {
                        //continue;
                    } else {
                        std::cout << "count is weird : (" << row << ", " << col << ", "
                                  << gsl_matrix_get(qmatrix, row, col) * sum << ")" << std::endl;
                        exit(-1);
                    }
                } else {
                    sum = 0.0;
                    for (int k = 0; k < 64; k++) {
                        default_element = gsl_complex_rect(0.0, 0.0);
                        for (int l = 0; l < 64; l++) {
                            default_element = gsl_complex_add(default_element, gsl_complex_mul(gsl_matrix_complex_get(eigenvector_com, col, l), gsl_matrix_complex_get(eigencount_com, k, l)));
                        }
                        sum += gsl_complex_mul(gsl_matrix_complex_get(eigenvec_inverse_com, k, row), default_element).dat[0];
                    }
                    if ((row != col && gsl_matrix_get(qmatrix, row, col) * sum >= 0.0) ||
                        (row == col && gsl_matrix_get(qmatrix, row, col) * sum <= 0)) {
                        gsl_matrix_set(total_count, row, col, gsl_matrix_get(qmatrix, row, col) * sum);
                    } else {
                        std::cout << "count is weird : (" << row << ", " << col << ", "
                                  << gsl_matrix_get(qmatrix, row, col) * sum << ")" << std::endl;
                        exit(-1);
                    }
                }
            }
        }
        for (int i = 0; i < 64; i++) {
            omega[i] = gsl_matrix_get(total_count, i, i) / gsl_matrix_get(qmatrix, i, i);
            //std::cout << "omega" << i << " : " << omega[i] << std::endl;
            if (abs(omega[i]) == 1.0e-14) {
                std::cout << "omega being 0, stop the process" << std::endl;
                exit(-1);
            }
        }

        //trying to make reversible
        for (int row = 0; row < 64; row++) {
            for (int col = 0; col < 64; col++) {
                //if (row > col) {
                    //Todo : testing without dividing the omega
                    gsl_matrix_set(qmatrix, row, col, gsl_matrix_get(total_count, row, col) / omega[row]);
                    //gsl_matrix_set(qmatrix, col, row, gsl_matrix_get(qmatrix, row, col) * coding_freq[row] / coding_freq[col]);
                //}
            }
        }

        for (int row = 0; row < 64; row++) {
            sum = 0.0;
            for (int col = 0; col < 64; col++) {
                if (row != col) {
                    sum += gsl_matrix_get(qmatrix, row, col);
                }
            }
            gsl_matrix_set(qmatrix, row, row, -sum);
        }

        /*double sum_codon_freq = 1.0;
        codon_freq[0] = 1.0;
        for (int i = 1; i < 64; i++) {
            codon_freq[i] = gsl_matrix_get(qmatrix, 0, i) / gsl_matrix_get(qmatrix, i, 0);
            sum_codon_freq += codon_freq[i];
        }
        for (int i = 0; i < 64; i++) {
            codon_freq[i] /= sum_codon_freq;
        }*/

        //normalize the matrix
        sum = 0.0;
        for (size_t dia = 0; dia < 64; dia++) {
            sum -= gsl_matrix_get(qmatrix, dia, dia) * coding_freq[dia];
        }
        gsl_matrix_scale(qmatrix, 1 / sum);

        for (size_t row = 0; row < 64; row++) {
            for (size_t col = 0; col < 64; col++) {
                std::cout << gsl_matrix_get(qmatrix, row, col) << ", ";
                if (row != col && gsl_matrix_get(qmatrix, row, col) < 0) {
                    std::cout << "nondiagonal being negative" << std::endl;
                    exit(-1);
                } else if (row == col && gsl_matrix_get(qmatrix, row, col) > 0) {
                    std::cout << "diagonal being positive" << std::endl;
                    exit(-1);
                }
            }
            std::cout << std::endl;
        }
        /////////////////////////////////////////////////////////////////////////
        std::cout << "likelihood : " << likelihood << std::endl;
        std::cout << "Difference : " << abs(likelihood - likelihood_old) << std::endl;
        count++;
    }

    gsl_matrix_free(eigenvector);
    gsl_matrix_free(eigenvec_inverse);
    gsl_vector_free(eigenvalue);
    gsl_matrix_free(total_count);
    gsl_matrix_free(eigencount);
    gsl_matrix_free(eigencount_col);
    gsl_vector_free(x);
    gsl_matrix_complex_free(eigenvector_com);
    gsl_vector_complex_free(eigenvalue_com);
    gsl_matrix_complex_free(eigenvec_inverse_com);
    gsl_matrix_free(qmatrix_temp);
    gsl_matrix_complex_free(eigencount_com);
    gsl_matrix_complex_free(eigencount_col_com);

    gsl_vector *scale = gsl_vector_alloc(64);
    for (int i = 0; i < 64; i++) {
        gsl_vector_set(scale, i, 1 / coding_freq[i]);
    }

    gsl_matrix_scale_columns(qmatrix, scale);
    gsl_vector_free(scale);

    for (size_t row = 0; row < 64; row++) {
        for (size_t col = 0; col < 64; col++) {
            std::cout << gsl_matrix_get(qmatrix, row, col) << ", ";
        }
        std::cout << std::endl;
    }
    auto end_time = std::chrono::steady_clock::now();
    std::cout << std::chrono::duration_cast<std::chrono::milliseconds>(end_time - start_time).count() / 1000.0 << std::endl;
    return 0;
}