diff --git a/model/modelmarkov.cpp b/model/modelmarkov.cpp index 6064b9d4b..36208859d 100644 --- a/model/modelmarkov.cpp +++ b/model/modelmarkov.cpp @@ -29,6 +29,10 @@ using namespace Eigen; /** number of squaring for scaling-squaring technique */ const int TimeSquare = 10; +/** Cassius : threshold for two eigenvalues to be too similar*/ +const double thr_eigenvalues = 1e-3; + + //----- declaration of some helper functions -----/ int matexp (double Q[], double t, int n, int TimeSquare, double space[]); int computeStateFreqFromQMatrix (double Q[], double pi[], int n, double space[]); @@ -115,6 +119,18 @@ void ModelMarkov::setReversible(bool reversible) { if (!ceval) ceval = aligned_alloc >(num_states); + + /*Cassius: I initialize the real counterparts, because I use them when the eigenvalues are real*/ + if (!eigenvalues) + eigenvalues = aligned_alloc(num_states); + if (!eigenvalues) + eigenvalues_imag = aligned_alloc(num_states); + if (!eigenvectors) + eigenvectors = aligned_alloc(num_states*num_states); + if (!inv_eigenvectors) + inv_eigenvectors = aligned_alloc(num_states*num_states); + /*Cassius: End of my edits*/ + if (!cevec) cevec = aligned_alloc >(num_states*num_states); if (!cinv_evec) @@ -389,9 +405,13 @@ void ModelMarkov::report_state_freqs(ostream& out, double *custom_state_freq) { void ModelMarkov::computeTransMatrix(double time, double *trans_matrix, int mixture) { if (!is_reversible) { - if (phylo_tree->params->matrix_exp_technique == MET_EIGEN_DECOMPOSITION) { + + // Cassius: If eigenvalues are different enough, use --eigen mehtod (closed formulae) + // This thr_eigenvalues is also used in decomposeRateMatrix() + // For consistency, the threshold must be the same in both instances + if (phylo_tree->params->matrix_exp_technique == MET_EIGEN_DECOMPOSITION and minimum_difference_eval > thr_eigenvalues and num_states == 4) { computeTransMatrixEigen(time, trans_matrix); - } else if (phylo_tree->params->matrix_exp_technique == MET_SCALING_SQUARING) { + } else if (phylo_tree->params->matrix_exp_technique == MET_SCALING_SQUARING or minimum_difference_eval <= thr_eigenvalues or num_states != 4) { // scaling and squaring technique int statesqr = num_states*num_states; memcpy(trans_matrix, rate_matrix, statesqr*sizeof(double)); @@ -405,46 +425,47 @@ void ModelMarkov::computeTransMatrix(double time, double *trans_matrix, int mixt // trans_matrix[i] *= time; // double space[NCODE*NCODE*3] = {0}; // matexp2(trans_matrix, num_states, 7, 5, space); - } - - /* compute P(t) */ - double evol_time = time / total_num_subst; - double exptime[num_states]; - int i, j, k; - - for (i = 0; i < num_states; i++) - exptime[i] = exp(evol_time * eigenvalues[i]); + } else { + /* compute P(t) */ + double evol_time = time / total_num_subst; + double exptime[num_states]; + int i, j, k; - int row_offset; - for (i = 0, row_offset = 0; i < num_states; i++, row_offset+=num_states) { - double *trans_row = trans_matrix + row_offset; - for (j = i+1; j < num_states; j ++) { - // compute upper triangle entries - double *trans_entry = trans_row + j; + for (i = 0; i < num_states; i++) + exptime[i] = exp(evol_time * eigenvalues[i]); + + int row_offset; + for (i = 0, row_offset = 0; i < num_states; i++, row_offset+=num_states) { + double *trans_row = trans_matrix + row_offset; + for (j = i+1; j < num_states; j ++) { + // compute upper triangle entries + double *trans_entry = trans_row + j; // double *coeff_entry = eigen_coeff + ((row_offset+j)*num_states); - *trans_entry = 0.0; - for (k = 0; k < num_states; k ++) { - *trans_entry += eigenvectors[i*num_states+k] * inv_eigenvectors[k*num_states+j] * exptime[k]; - } - if (*trans_entry < 0.0) { - *trans_entry = 0.0; - } - // update lower triangle entries - trans_matrix[j*num_states+i] = (state_freq[i]/state_freq[j]) * (*trans_entry); - } - trans_row[i] = 0.0; // initialize diagonal entry - // taking the sum of row - double sum = 0.0; - for (j = 0; j < num_states; j++) - sum += trans_row[j]; - trans_row[i] = 1.0 - sum; // update diagonal entry - } + *trans_entry = 0.0; + for (k = 0; k < num_states; k ++) { + *trans_entry += eigenvectors[i*num_states+k] * inv_eigenvectors[k*num_states+j] * exptime[k]; + } + if (*trans_entry < 0.0) { + *trans_entry = 0.0; + } + // update lower triangle entries + trans_matrix[j*num_states+i] = (state_freq[i]/state_freq[j]) * (*trans_entry); + } + trans_row[i] = 0.0; // initialize diagonal entry + // taking the sum of row + double sum = 0.0; + for (j = 0; j < num_states; j++) + sum += trans_row[j]; + trans_row[i] = 1.0 - sum; // update diagonal entry + } // delete [] exptime; + } } double ModelMarkov::computeTrans(double time, int state1, int state2) { if (is_reversible) { + double evol_time = time / total_num_subst; int i; double trans_prob = 0.0; @@ -867,11 +888,8 @@ void ModelMarkov::decomposeRateMatrix(){ if (!is_reversible) { double sum; - //double m[num_states]; - double *space = new double[num_states*(num_states+1)]; - for (i = 0; i < num_states; i++) - state_freq[i] = 1.0/num_states; + auto *space = new double[num_states*(num_states+1)]; for (i = 0, k = 0; i < num_states; i++) { rate_matrix[i*num_states+i] = 0.0; @@ -882,6 +900,7 @@ void ModelMarkov::decomposeRateMatrix(){ } rate_matrix[i*num_states+i] = -row_sum; } + computeStateFreqFromQMatrix(rate_matrix, state_freq, num_states, space); @@ -900,56 +919,91 @@ void ModelMarkov::decomposeRateMatrix(){ } delete [] space; - if (phylo_tree->params->matrix_exp_technique == MET_EIGEN_DECOMPOSITION) { - eigensystem_nonrev(rate_matrix, state_freq, eigenvalues, eigenvalues_imag, eigenvectors, inv_eigenvectors, num_states); + + // Begin of Cassius' edits + + // The closed formulae only apply for 4 states, for nonreversible and different-than-4 states there is no implementation + if (num_states == 4) { + // Cassius: We use the closed formula for the eigenvalues + eigenvalues_closed_formula_nonrev(rate_matrix, ceval, eigenvalues, eigenvalues_imag); + all_eigenvalues_real = true; + for (int i = 0; i < num_states; ++i) + if (abs(eigenvalues_imag[i]) > 1e-6) all_eigenvalues_real = false; + + minimum_difference_eval = 1000; + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < i; ++j) { + minimum_difference_eval = min (minimum_difference_eval, abs(ceval[i] - ceval[j])); + } + } + + //Cassius: When two eigenvalues are equal, we should forget eigenvectors and go for the squaring technique in computeTransMatrix + if (phylo_tree->params->matrix_exp_technique == MET_EIGEN_DECOMPOSITION and minimum_difference_eval > thr_eigenvalues) { + if (all_eigenvalues_real) { + eigensystem_closed_formula_nonrev_real(rate_matrix, state_freq, eigenvalues, eigenvectors, inv_eigenvectors); + } else { + eigensystem_closed_formula_nonrev(rate_matrix, state_freq, ceval, cevec, cinv_evec); + } + } else { + /*Our effort was useless, we will use approximate the exponential using the squaring technique*/ + //ASSERT(0 && "this line should not be reached"); + return; + } + } else { + // Cassius: I commented the following line because this function is from the previous implementation and doesn't work. + // eigensystem_nonrev(rate_matrix, state_freq, eigenvalues, eigenvalues_imag, eigenvectors, inv_eigenvectors, num_states); + return; } + + /*end of Cassius edits*/ + } else if (num_params == -1) { // reversible model - // manual compute eigenvalues/vectors for F81-style model - eigenvalues[0] = 0.0; - double mu = 0.0; - for (i = 0; i < num_states; i++) - mu += state_freq[i]*state_freq[i]; - mu = total_num_subst/(1.0 - mu); - - // compute eigenvalues - for (i = 1; i < num_states; i++) - eigenvalues[i] = -mu; - -// double *f = new double[num_states]; -// for (i = 0; i < num_states; i++) f[i] = sqrt(state_freq[i]); - // compute eigenvectors - memset(eigenvectors, 0, num_states*num_states*sizeof(double)); - memset(inv_eigenvectors, 0, num_states*num_states*sizeof(double)); - eigenvectors[0] = 1.0; - for (i = 1; i < num_states; i++) - eigenvectors[i] = -1.0; -// eigenvectors[i] = f[i]/f[num_states-1]; - for (i = 1; i < num_states; i++) { - eigenvectors[i*num_states] = 1.0; - eigenvectors[i*num_states+i] = state_freq[0]/state_freq[i]; - } + // manual compute eigenvalues/vectors for F81-style model + eigenvalues[0] = 0.0; + double mu = 0.0; + for (i = 0; i < num_states; i++) + mu += state_freq[i]*state_freq[i]; + mu = total_num_subst/(1.0 - mu); + + // compute eigenvalues + for (i = 1; i < num_states; i++) + eigenvalues[i] = -mu; + + // double *f = new double[num_states]; + // for (i = 0; i < num_states; i++) f[i] = sqrt(state_freq[i]); + // compute eigenvectors + memset(eigenvectors, 0, num_states*num_states*sizeof(double)); + memset(inv_eigenvectors, 0, num_states*num_states*sizeof(double)); + eigenvectors[0] = 1.0; + for (i = 1; i < num_states; i++) + eigenvectors[i] = -1.0; + // eigenvectors[i] = f[i]/f[num_states-1]; + for (i = 1; i < num_states; i++) { + eigenvectors[i*num_states] = 1.0; + eigenvectors[i*num_states+i] = state_freq[0]/state_freq[i]; + } - for (i = 0; i < num_states; i++) - for (j = 0; j < num_states; j++) - inv_eigenvectors[i*num_states+j] = state_freq[j]*eigenvectors[j*num_states+i]; - writeInfo(cout); - // sanity check - double *q = new double[num_states*num_states]; - getQMatrix(q); - double zero; - for (j = 0; j < num_states; j++) { - for (i = 0, zero = 0.0; i < num_states; i++) { - for (k = 0; k < num_states; k++) zero += q[i*num_states+k] * eigenvectors[k*num_states+j]; - zero -= eigenvalues[j] * eigenvectors[i*num_states+j]; - if (fabs(zero) > 1.0e-5) { - cout << "\nERROR: Eigenvector doesn't satisfy eigenvalue equation! (gap=" << fabs(zero) << ")" << endl; - abort(); - } - } - } - delete [] q; -/* } else if (Params::getInstance().matrix_exp_technique == MET_EIGEN3LIB_DECOMPOSITION) { + for (i = 0; i < num_states; i++) + for (j = 0; j < num_states; j++) + inv_eigenvectors[i*num_states+j] = state_freq[j]*eigenvectors[j*num_states+i]; + writeInfo(cout); + // sanity check + double *q = new double[num_states*num_states]; + getQMatrix(q); + double zero; + for (j = 0; j < num_states; j++) { + for (i = 0, zero = 0.0; i < num_states; i++) { + for (k = 0; k < num_states; k++) zero += q[i*num_states+k] * eigenvectors[k*num_states+j]; + zero -= eigenvalues[j] * eigenvectors[i*num_states+j]; + if (fabs(zero) > 1.0e-5) { + cout << "\nERROR: Eigenvector doesn't satisfy eigenvalue equation! (gap=" << fabs(zero) << ")" << endl; + abort(); + } + } + } + delete [] q; + /* } else if (Params::getInstance().matrix_exp_technique == MET_EIGEN3LIB_DECOMPOSITION) { // Using Eigen3 libary for general time-reversible model Eigen::MatrixXd Q; Q.resize(num_states, num_states); @@ -961,37 +1015,83 @@ void ModelMarkov::decomposeRateMatrix(){ } else { // full matrix } -*/ - } else { + */ + } else { // general reversible model - double **rate_matrix = new double*[num_states]; + //the rates are a symmetrix matrix, don't confuse with the rate matrix + double **matrix_rates = new double *[num_states]; - for (i = 0; i < num_states; i++) - rate_matrix[i] = new double[num_states]; + for (i = 0; i < num_states; i++) + matrix_rates[i] = new double[num_states]; if (half_matrix) { for (i = 0, k = 0; i < num_states; i++) { - rate_matrix[i][i] = 0.0; - for (j = i+1; j < num_states; j++, k++) { - rate_matrix[i][j] = (state_freq[i] <= ZERO_FREQ || state_freq[j] <= ZERO_FREQ) ? 0 : rates[k]; - rate_matrix[j][i] = rate_matrix[i][j]; + matrix_rates[i][i] = 0.0; + for (j = i + 1; j < num_states; j++, k++) { + matrix_rates[i][j] = (state_freq[i] <= ZERO_FREQ || state_freq[j] <= ZERO_FREQ) ? 0 : rates[k]; + matrix_rates[j][i] = matrix_rates[i][j]; } } } else { // full matrix for (i = 0; i < num_states; i++) { - memcpy(rate_matrix[i], &rates[i*num_states], num_states*sizeof(double)); - rate_matrix[i][i] = 0.0; + memcpy(matrix_rates[i], &rates[i * num_states], num_states * sizeof(double)); + matrix_rates[i][i] = 0.0; } } - /* eigensystem of 1 PAM rate matrix */ - eigensystem_sym(rate_matrix, state_freq, eigenvalues, eigenvectors, inv_eigenvectors, num_states); - //eigensystem(rate_matrix, state_freq, eigenvalues, eigenvectors, inv_eigenvectors, num_states); - for (i = num_states-1; i >= 0; i--) - delete [] rate_matrix[i]; - delete [] rate_matrix; - } + + /*Cassius' part begins*/ + + if (phylo_tree->params->matrix_exp_technique == MET_EIGEN_DECOMPOSITION and num_states == 4) { + + if (!rate_matrix) + rate_matrix = aligned_alloc(num_states*num_states); + //compute the actual rate matrix + for (i = 0; i < num_states; ++i) { + for (j = 0; j < num_states; ++j) rate_matrix[num_states*i+j] = matrix_rates[i][j]*state_freq[j]; + } + + //rows sum zero + for (i = 0; i < num_states; ++i) { + for (j = i+1; j < num_states; ++j) rate_matrix[num_states*i + i] -= rate_matrix[num_states*i + j]; + for (j = 0; j < i; ++j) rate_matrix[num_states*i + i] -= rate_matrix[num_states*i + j]; + } + //average mutation equal 1 + double average_mut = 0; + for (i = 0; i < num_states; ++i) average_mut -= rate_matrix[num_states*i+i]*state_freq[i]; + for (i = 0; i < num_states; ++i) { + for (j = 0; j < num_states; ++j) rate_matrix[num_states*i+j] /= average_mut; + } + + // Compute eigenvalues. It holds that eigenvalues[3] == 0 + eigenvalues_closed_formula_rev(rate_matrix, eigenvalues); + + + //only if eigenvalues are different we can apply the closed formulas + minimum_difference_eval = 1000; + for (int i = 0; i < num_states; ++i) { + for (int j = 0; j < i; ++j) { + minimum_difference_eval = min (minimum_difference_eval, abs(eigenvalues[i] - eigenvalues[j])); + } + } + + if (minimum_difference_eval > thr_eigenvalues) + eigensystem_closed_formula_rev(rate_matrix, state_freq, eigenvalues, eigenvectors, inv_eigenvectors); + else + eigensystem_sym(matrix_rates, state_freq, eigenvalues, eigenvectors, inv_eigenvectors, num_states); + + // Cassius' part ends. + + } else { + /* eigensystem of 1 PAM rate matrix */ + eigensystem_sym(matrix_rates, state_freq, eigenvalues, eigenvectors, inv_eigenvectors, num_states); + //eigensystem(matrix_rates, state_freq, eigenvalues, eigenvectors, inv_eigenvectors, num_states); + } + for (i = num_states - 1; i >= 0; i--) + delete[] matrix_rates[i]; + delete[] matrix_rates; + } } void ModelMarkov::readRates(istream &in) throw(const char*, string) { @@ -1300,56 +1400,79 @@ void ModelMarkov::update_eigen_pointers(double *eval, double *evec, double *inv_ } void ModelMarkov::computeTransMatrixEigen(double time, double *trans_matrix) { - /* compute P(t) */ - double evol_time = time / total_num_subst; - int nstates_2 = num_states*num_states; - double *exptime = new double[nstates_2]; - int i, j, k; - memset(exptime, 0, sizeof(double)*nstates_2); - for (i = 0; i < num_states; i++) - if (eigenvalues_imag[i] == 0.0) { - exptime[i*num_states+i] = exp(evol_time * eigenvalues[i]); - } else { - ASSERT(i < num_states-1 && eigenvalues_imag[i+1] != 0.0 && eigenvalues_imag[i] > 0.0); - complex exp_eval(eigenvalues[i] * evol_time, eigenvalues_imag[i] * evol_time); - exp_eval = exp(exp_eval); - exptime[i*num_states+i] = exp_eval.real(); - exptime[i*num_states+i+1] = exp_eval.imag(); - i++; - exptime[i*num_states+i] = exp_eval.real(); - exptime[i*num_states+i-1] = -exp_eval.imag(); - } + if (all_eigenvalues_real) { + /* compute P(t) assuming everything is real */ + double evol_time = time / total_num_subst; + int nstates_2 = num_states*num_states; + auto *exptime = new double[num_states]; + auto *trans_aux = new double[nstates_2]; + int i, j, k; - // compute V * exp(L t) - for (i = 0; i < num_states; i++) - for (j = 0; j < num_states; j++) { - double val = 0; - for (k = 0; k < num_states; k++) - val += eigenvectors[i*num_states+k] * exptime[k*num_states+j]; - trans_matrix[i*num_states+j] = val; + + memset(exptime, 0, sizeof(double)*num_states); + for (i = 0; i < num_states; i++) + exptime[i] = exp(eigenvalues[i]*evol_time); + + memset(trans_aux, 0, sizeof(double)*nstates_2); + // Compute V * exp(L t) * V^{-1} + for (i = 0; i < num_states; i++) + for (j = 0; j < num_states; j++) { + for (k = 0; k < num_states; ++k) + trans_aux[num_states*i + j] += eigenvectors[num_states*i + k] * exptime[k]*inv_eigenvectors[num_states*k + j]; + } + + for (i = 0; i < num_states; i++) { + double row_sum = 0.0; + for (j = 0; j < num_states; j++) { + trans_matrix[num_states * i + j] = trans_aux[num_states * i + j]; + row_sum += trans_matrix[num_states * i + j]; + ASSERT(trans_matrix[num_states * i + j] >= -0.001); + } + ASSERT(fabs(row_sum - 1.0) < 1e-4); } - memcpy(exptime, trans_matrix, sizeof(double)*nstates_2); - // then compute V * exp(L t) * V^{-1} - for (i = 0; i < num_states; i++) { - double row_sum = 0.0; - for (j = 0; j < num_states; j++) { - double val = 0; - for (k = 0; k < num_states; k++) - val += exptime[i*num_states+k] * inv_eigenvectors[k*num_states+j]; - // make sure that trans_matrix are non-negative - ASSERT(val >= -0.001); - val = fabs(val); - trans_matrix[i*num_states+j] = val; - row_sum += val; + delete [] trans_aux; + delete [] exptime; + + } else { + /* compute P(t) assuming everything is complex */ + + double evol_time = time / total_num_subst; + int nstates_2 = num_states*num_states; + auto *exptime = new complex[num_states]; + auto *trans_aux = new complex[nstates_2]; + int i, j, k; + + + memset(exptime, 0, sizeof(std::complex)*num_states); + for (i = 0; i < num_states; i++) + exptime[i] = exp(ceval[i]*evol_time); + + memset(trans_aux, 0, sizeof(std::complex)*nstates_2); + // Compute V * exp(L t) * V^{-1} + for (i = 0; i < num_states; i++) + for (j = 0; j < num_states; j++) { + for (k = 0; k < num_states; ++k) + trans_aux[num_states*i + j] += cevec[num_states*i + k] * exptime[k]*cinv_evec[num_states*k + j]; + } + + for (i = 0; i < num_states; i++) { + double row_sum = 0.0; + for (j = 0; j < num_states; j++) { + trans_matrix[num_states * i + j] = trans_aux[num_states * i + j].real(); + row_sum += trans_matrix[num_states * i + j]; + ASSERT(abs(trans_aux[num_states * i + j].imag()) < 1e-6); + ASSERT(trans_matrix[num_states * i + j] >= -0.001); + } + ASSERT(fabs(row_sum - 1.0) < 1e-4); } - ASSERT(fabs(row_sum-1.0) < 1e-4); - } - delete [] exptime; + delete [] trans_aux; + delete [] exptime; + } } @@ -1490,3 +1613,5 @@ int matexp (double Q[], double t, int n, int TimeSquare, double space[]) for (i=0;i values if necessary compute the transition probability matrix using (complex) eigenvalues @param time time between two events @param trans_matrix (OUT) the transition matrix between all pairs of states. @@ -484,6 +484,17 @@ class ModelMarkov : public ModelSubst, public EigenDecomposition */ std::complex *ceval, *cevec, *cinv_evec; + /** + Cassius: equals true if all eigenvalues are real + */ + bool all_eigenvalues_real; + + /** + Cassius: minimum difference between two eigenvalues +*/ + double minimum_difference_eval; + + }; #endif diff --git a/utils/eigendecomposition.cpp b/utils/eigendecomposition.cpp index 699e17f8a..2f8a45f1e 100644 --- a/utils/eigendecomposition.cpp +++ b/utils/eigendecomposition.cpp @@ -1300,6 +1300,7 @@ void EigenDecomposition::checkevector(double *evec, double *ivec, int nn) { } } } + if (error) { cout << "\nWARNING: Inversion of eigenvector matrix not perfect!\n"; } @@ -1308,3 +1309,610 @@ void EigenDecomposition::checkevector(double *evec, double *ivec, int nn) { delete [] matx[i]; delete [] matx; } /* checkevector */ + + + +/*Cassius: From now on, these are my functions to implement the closed formulae of the decomposition*/ + +/* check inversion with complex entries*/ +void EigenDecomposition::checkevector_complex(complex *evec, complex *ivec, int nn) { + int i, j, ia, ib, ic, error; + complex **matx = (complex**) new complex [nn]; + complex sum; + + for (i = 0; i < nn; i++) + matx[i] = new complex[nn]; + + /* multiply matrix of eigenvectors and its inverse */ + for (ia = 0; ia < nn; ia++) { + for (ic = 0; ic < nn; ic++) { + sum = 0.0; + for (ib = 0; ib < nn; ib++) sum += evec[ia*nn+ib] * ivec[ib*nn+ic]; + matx[ia][ic] = sum; + } + } + /* check whether the unitary matrix is obtained */ + error = 0; + for (i = 0; i < nn; i++) { + for (j = 0; j < nn; j++) { + if (i == j) { + if (abs(matx[i][j] - 1.0) > 1.0e-5) + error = 1; + } else { + if (abs(matx[i][j]) > 1.0e-5) + error = 1; + } + } + } + + if (error) { + cout << "\nWARNING: Inversion of eigenvector matrix not perfect!\n"; + } + + for (i = nn-1; i >= 0; i--) + delete [] matx[i]; + delete [] matx; +} + + +/* eigensystem */ + +/* solve the cubic equation*/ + +void EigenDecomposition::solve_cubic_equation(double a, double b, double c, double d, double* root_real, double* root_imag) { + const double PI = 4.0 * atan( 1.0 ); + + // Reduced equation: X^3 - 3pX - 2q = 0, where X = x-b/(3a) + double p = ( b * b - 3.0 * a * c ) / ( 9.0 * a * a ); + double q = ( 9.0 * a * b * c - 27.0 * a * a * d - 2.0 * b * b * b ) / ( 54.0 * a * a * a ); + double offset = b / ( 3.0 * a ); + + // Discriminant + double discriminant = p * p * p - q * q; + + if (abs(discriminant) < 1e-12) + discriminant = 0; + if (abs(q) < 1e-18) + q = 0; + + if ( discriminant > 0) // set X = 2 sqrt(p) cos(theta) and compare 4 cos^3(theta)-3 cos(theta) = cos(3 theta) + { + double theta = acos( q / ( p * sqrt( p ) ) ); + double r = 2.0 * sqrt( p ); + for ( int n = 0; n < 3; n++ ) { + root_real[n] = r * cos( ( theta + 2.0 * n * PI ) / 3.0 ) - offset; + root_imag[n] = 0; + } + } + else + { + double gamma1 = cbrt( q + sqrt( -discriminant ) ); + double gamma2 = cbrt( q - sqrt( -discriminant ) ); + + root_real[0] = gamma1 + gamma2 - offset; + root_imag[0] = 0; + + double re = -0.5 * ( gamma1 + gamma2 ) - offset; + double im = ( gamma1 - gamma2 ) * sqrt( 3.0 ) / 2.0; + if ( discriminant == 0 ) // Equal roots + { + root_real[1] = root_real[2] = re; + root_imag[1] = root_imag[2] = 0; + } + else + { + root_real[1] = root_real[2] = re; + root_imag[1] = im; + root_imag[2] = -im; + } + } +} + +/* For NONREVERSIBLE matrices WITH 4 STATES, find the eigenvalues using the closed formulae */ +void EigenDecomposition::eigenvalues_closed_formula_nonrev( + double *rate_matrix, complex *ceval, double *eval, double *ieval) +{ + auto *root_real = new double[3]; + auto *root_imag = new double[3]; + + double c0 = 0; + double c1 = 0; + double c2 = 0; + for (int i = 0; i < 4; i++){ + for (int j = 0; j < i; j++) { + c2 += rate_matrix[4*i + j]; + } + for (int j = i+1; j < 4; j++) { + c2 += rate_matrix[4*i + j]; + } + } + + //The following are the parameters as output by Mathematica using CForm. I extract common factor whenever possible, to minimize computation. + + c1 = rate_matrix[11]*(rate_matrix[12] + rate_matrix[13] + rate_matrix[1] + rate_matrix[2] + rate_matrix[3] + rate_matrix[4] + rate_matrix[6] + rate_matrix[7]) + (rate_matrix[4] + rate_matrix[6])*rate_matrix[8] + rate_matrix[4]*rate_matrix[9] + + rate_matrix[2]*(rate_matrix[4] + rate_matrix[6] + rate_matrix[7] + rate_matrix[9]) + rate_matrix[7]*(rate_matrix[8] + rate_matrix[9]) + + rate_matrix[13]*(rate_matrix[1] + rate_matrix[2] + rate_matrix[3] + rate_matrix[4] + rate_matrix[6] + rate_matrix[8] + rate_matrix[9]) + rate_matrix[1]*(rate_matrix[6] + rate_matrix[7] + rate_matrix[8] + rate_matrix[9]) + + rate_matrix[3]*(rate_matrix[4] + rate_matrix[6] + rate_matrix[7] + rate_matrix[8] + rate_matrix[9]) + rate_matrix[12]*(rate_matrix[1] + rate_matrix[2] + rate_matrix[4] + rate_matrix[6] + rate_matrix[7] + rate_matrix[8] + rate_matrix[9]) + + rate_matrix[14]*(rate_matrix[1] + rate_matrix[2] + rate_matrix[3] + rate_matrix[4] + rate_matrix[6] + rate_matrix[7] + rate_matrix[8] + rate_matrix[9]); + + c0 = rate_matrix[11]*(rate_matrix[13]*(rate_matrix[1] + rate_matrix[2] + rate_matrix[3] + rate_matrix[4]) + rate_matrix[1]*(rate_matrix[6] + rate_matrix[7]) + + rate_matrix[2]*(rate_matrix[4] + rate_matrix[6] + rate_matrix[7]) + rate_matrix[3]*(rate_matrix[4] + rate_matrix[6] + rate_matrix[7]) + + rate_matrix[12]*(rate_matrix[1] + rate_matrix[4] + rate_matrix[6] + rate_matrix[7])) + rate_matrix[2]*rate_matrix[7]*rate_matrix[9] + + rate_matrix[1]*rate_matrix[7]*(rate_matrix[8] + rate_matrix[9]) + rate_matrix[3]*((rate_matrix[4] + rate_matrix[6])*rate_matrix[8] + rate_matrix[4]*rate_matrix[9] + + rate_matrix[7]*(rate_matrix[8] + rate_matrix[9])) + rate_matrix[12]*((rate_matrix[4] + rate_matrix[6])*rate_matrix[8] + rate_matrix[4]*rate_matrix[9] + + rate_matrix[2]*(rate_matrix[4] + rate_matrix[6] + rate_matrix[7] + rate_matrix[9]) + rate_matrix[7]*(rate_matrix[8] + rate_matrix[9]) + rate_matrix[1]*(rate_matrix[6] + rate_matrix[8] + rate_matrix[9])) + + rate_matrix[13]*((rate_matrix[4] + rate_matrix[6])*rate_matrix[8] + rate_matrix[4]*rate_matrix[9] + rate_matrix[2]*(rate_matrix[4] + + rate_matrix[6] + rate_matrix[9]) + rate_matrix[1]*(rate_matrix[6] + rate_matrix[8] + rate_matrix[9]) + rate_matrix[3]*(rate_matrix[6] + rate_matrix[8] + rate_matrix[9])) + + rate_matrix[14]*((rate_matrix[4] + rate_matrix[6])*rate_matrix[8] + rate_matrix[7]*rate_matrix[8] + rate_matrix[4]*rate_matrix[9] + rate_matrix[2]*(rate_matrix[4] + rate_matrix[6] + + rate_matrix[7] + rate_matrix[9]) + rate_matrix[3]*(rate_matrix[4] + rate_matrix[6] + rate_matrix[7] + rate_matrix[9]) + rate_matrix[1]*(rate_matrix[6] + rate_matrix[7] + rate_matrix[8] + rate_matrix[9])); + + solve_cubic_equation(1, c2, c1, c0, root_real, root_imag); + + + + ceval[3] = (0.0, 0.0); + eval[3] = 0.0; + ieval[3] = 0.0; + for (int i = 0; i < 3; ++i) { + complex aux(root_real[i], root_imag[i]); + ceval[i] = aux; + eval[i] = root_real[i]; + ieval[i] = root_imag[i]; + } + + delete [] root_real; + delete [] root_imag; + +} + +/* Find the eigenvalues of a REVERSIBLE rate matrix with DIMENSION 4 using closed formulae*/ +void EigenDecomposition::eigenvalues_closed_formula_rev( + double *rate_matrix, double *eval) +{ + + auto *root_real = new double[3]; + auto *root_imag = new double[3]; + + double c0; + double c1; + double c2 = 0; + for (int i = 0; i < 4; i++){ + for (int j = 0; j < i; j++) { + c2 += rate_matrix[4*i + j]; + } + for (int j = i+1; j < 4; j++) { + c2 += rate_matrix[4*i + j]; + } + } + + //The following are the parameters as output by Mathematica using CForm. I extract common factor whenever possible, to minimize computation. + + c1 = rate_matrix[11]*(rate_matrix[12] + rate_matrix[13] + rate_matrix[1] + rate_matrix[2] + rate_matrix[3] + rate_matrix[4] + rate_matrix[6] + rate_matrix[7]) + (rate_matrix[4] + rate_matrix[6])*rate_matrix[8] + rate_matrix[4]*rate_matrix[9] + + rate_matrix[2]*(rate_matrix[4] + rate_matrix[6] + rate_matrix[7] + rate_matrix[9]) + rate_matrix[7]*(rate_matrix[8] + rate_matrix[9]) + + rate_matrix[13]*(rate_matrix[1] + rate_matrix[2] + rate_matrix[3] + rate_matrix[4] + rate_matrix[6] + rate_matrix[8] + rate_matrix[9]) + rate_matrix[1]*(rate_matrix[6] + rate_matrix[7] + rate_matrix[8] + rate_matrix[9]) + + rate_matrix[3]*(rate_matrix[4] + rate_matrix[6] + rate_matrix[7] + rate_matrix[8] + rate_matrix[9]) + rate_matrix[12]*(rate_matrix[1] + rate_matrix[2] + rate_matrix[4] + rate_matrix[6] + rate_matrix[7] + rate_matrix[8] + rate_matrix[9]) + + rate_matrix[14]*(rate_matrix[1] + rate_matrix[2] + rate_matrix[3] + rate_matrix[4] + rate_matrix[6] + rate_matrix[7] + rate_matrix[8] + rate_matrix[9]); + + c0 = rate_matrix[11]*(rate_matrix[13]*(rate_matrix[1] + rate_matrix[2] + rate_matrix[3] + rate_matrix[4]) + rate_matrix[1]*(rate_matrix[6] + rate_matrix[7]) + + rate_matrix[2]*(rate_matrix[4] + rate_matrix[6] + rate_matrix[7]) + rate_matrix[3]*(rate_matrix[4] + rate_matrix[6] + rate_matrix[7]) + + rate_matrix[12]*(rate_matrix[1] + rate_matrix[4] + rate_matrix[6] + rate_matrix[7])) + rate_matrix[2]*rate_matrix[7]*rate_matrix[9] + + rate_matrix[1]*rate_matrix[7]*(rate_matrix[8] + rate_matrix[9]) + rate_matrix[3]*((rate_matrix[4] + rate_matrix[6])*rate_matrix[8] + rate_matrix[4]*rate_matrix[9] + + rate_matrix[7]*(rate_matrix[8] + rate_matrix[9])) + rate_matrix[12]*((rate_matrix[4] + rate_matrix[6])*rate_matrix[8] + rate_matrix[4]*rate_matrix[9] + + rate_matrix[2]*(rate_matrix[4] + rate_matrix[6] + rate_matrix[7] + rate_matrix[9]) + rate_matrix[7]*(rate_matrix[8] + rate_matrix[9]) + rate_matrix[1]*(rate_matrix[6] + + rate_matrix[8] + rate_matrix[9])) + rate_matrix[13]*((rate_matrix[4] + rate_matrix[6])*rate_matrix[8] + rate_matrix[4]*rate_matrix[9] + rate_matrix[2]*(rate_matrix[4] + + rate_matrix[6] + rate_matrix[9]) + rate_matrix[1]*(rate_matrix[6] + rate_matrix[8] + rate_matrix[9]) + rate_matrix[3]*(rate_matrix[6] + rate_matrix[8] + rate_matrix[9])) + + rate_matrix[14]*((rate_matrix[4] + rate_matrix[6])*rate_matrix[8] + rate_matrix[7]*rate_matrix[8] + rate_matrix[4]*rate_matrix[9] + rate_matrix[2]*(rate_matrix[4] + rate_matrix[6] + + rate_matrix[7] + rate_matrix[9]) + rate_matrix[3]*(rate_matrix[4] + rate_matrix[6] + rate_matrix[7] + rate_matrix[9]) + rate_matrix[1]*(rate_matrix[6] + rate_matrix[7] + rate_matrix[8] + rate_matrix[9])); + + solve_cubic_equation(1, c2, c1, c0, root_real, root_imag); + + eval[3] = 0.0; + for (int i = 0; i < 3; ++i) + eval[i] = root_real[i]; + + delete [] root_real; + delete [] root_imag; +} + +/* Closed formula for the eigenvectors of a (nonreversible) rate matrix with DIMENSION 4 and DIFFERENT EIGENVALUES */ + +void EigenDecomposition::eigensystem_closed_formula_nonrev(double *rate_matrix, double *state_freq, complex *ceval, + complex *cevec, complex *cinv_evec) { + + int i,j,k; + + double square_rate_mat[16]; + for (i = 0; i < 4; ++i) { + for (j = 0; j < 4; ++j) { + square_rate_mat[4*i + j] = 0.; + for (k = 0; k < 4; ++k) square_rate_mat[4*i + j] += rate_matrix[4*i + k]*rate_matrix[4*k + j]; + } + } + + + // Declare the right eigenvectors using close formulae + + for (i = 0; i < 4; ++i) + cevec[4*i + 3] = 1.0; + + for (k = 0; k < 3; ++k) { + int pre = (k + 2) % 3; + int post = (k + 1) % 3; + for (i = 0; i < 4; ++i) { + cevec[4*i + k] = ceval[pre]*ceval[post]*( bool(i==0) - state_freq[0]) - (ceval[pre] + ceval[post])*rate_matrix[4*i] + square_rate_mat[4*i]; + } + } + + // Declare the left eigenvectors using close formulae + + for (i = 0; i < 4; ++i) + cinv_evec[4*3 + i] = state_freq[i]; + + for (k = 0; k < 3; ++k) { + int pre = (k + 2) % 3; + int post = (k + 1) % 3; + for (i = 0; i < 4; ++i) { + cinv_evec[4*k + i] = ceval[pre]*ceval[post]*( bool(i==0) - state_freq[i]) - (ceval[pre] + ceval[post])*rate_matrix[i] + square_rate_mat[i]; + } + } + + // Rescale the eigenvectors so that inv_evec is the inverse of evec. This is arbitrary up to a constant factor. + + complex root_scalar_products[4]; + for (i = 0; i < 4; ++i) { + root_scalar_products[i] = 0.; + for (j = 0; j < 4; ++j) root_scalar_products[i] += cevec[4*j + i] *cinv_evec[4*i + j]; + if (abs(root_scalar_products[i]) < 1e-6) { + cout << "Warning: The module of one eigenvector was arbitrarily set too small." << endl; + cout << "Here comes the right vector" << endl; + for (j = 0; j < 4; ++j) { + cout << cevec[4*j + i] << " "; + } + cout << endl; + cout << "Here comes the left vector" << endl; + for (j = 0; j < 4; ++j) { + cout << cinv_evec[4*i + j] << " "; + } + cout << endl; + + } + root_scalar_products[i] = sqrt(root_scalar_products[i]); + } + + for (i = 0; i < 4; ++i) + for (j = 0; j < 4; ++j) cevec[4 * i + j] /= root_scalar_products[j]; + + + for (i = 0; i < 4; ++i) + for (j = 0; j < 4; ++j) cinv_evec[4 * i + j] /= root_scalar_products[i]; + + + checkevector_complex(cevec, cinv_evec, 4); /* check whether inversion was OK */ + +} + +/* Closed formula for the eigenvectors of a rate matrix with DIMENSION 4 and DIFFERENT REAL eigenvalues */ + +void EigenDecomposition::eigensystem_closed_formula_nonrev_real(double *rate_matrix, double *state_freq, double *eval, double *evec, double *inv_evec) { + + // Set true to output all interesting parts of the process + bool verbose = false; + + int i,j,k; + + // Some threshold for very small numbers + double small = 1e-6; + + auto *square_rate_mat = new double [16]; + memset(square_rate_mat, 0, sizeof(double)*16); + + for (i = 0; i < 4; ++i) { + for (j = 0; j < 4; ++j) { + for (k = 0; k < 4; ++k) square_rate_mat[4*i + j] += rate_matrix[4*i + k]*rate_matrix[4*k + j]; + } + } + + vector got_right(3, false); + vector got_left(3, false); + + // Declare the right eigenvectors using close formulae + + + for (i = 0; i < 4; ++i) + evec[4*i + 3] = 1.0; + + for (k = 0; k < 3; ++k) { + int pre = (k + 2) % 3; + int post = (k + 1) % 3; + // the loop for "candidate" iterates looking for a nonzero eigenvector + for (int candidate = 0; candidate < 4 and not got_right[k]; ++ candidate) { + for (i = 0; i < 4; ++i) { + evec[4 * i + k] = eval[pre] * eval[post] * (bool(i == candidate) - state_freq[candidate]) - (eval[pre] + eval[post]) * rate_matrix[4 * i + candidate] + square_rate_mat[4 * i + candidate]; + if (abs(evec[4 * i + k]) > small) got_right[k] = true; + } + } + } + + + + // Declare the left eigenvectors using close formulae + + for (i = 0; i < 4; ++i) + inv_evec[4*3 + i] = state_freq[i]; + + + for (k = 0; k < 3; ++k) { + int pre = (k + 2) % 3; + int post = (k + 1) % 3; + // the loop for "candidate" iterates looking for a nonzero eigenvector + for (int candidate = 0; candidate < 4 and not got_left[k]; ++ candidate) { + for (i = 0; i < 4; ++i) { + inv_evec[4*k + i] = eval[pre]*eval[post]*( bool(i == candidate) - state_freq[i]) - (eval[pre] + eval[post])*rate_matrix[4* candidate + i] + square_rate_mat[4*candidate + i]; + if (abs(inv_evec[4*k + i]) > small) got_left[k] = true; + } + } + } + + + // Rescale the eigenvectors so that inv_evec is the inverse of evec. This is arbitrary up to a constant factor. + + auto *scalar_products = new double[4]; + memset(scalar_products, 0, sizeof(double)*4); + + for (i = 0; i < 4; ++i) { + for (j = 0; j < 4; ++j) scalar_products[i] += evec[4*j + i] *inv_evec[4*i + j]; + + if (abs(scalar_products[i]) < 1e-6) { + cout << "Warning: The module of one eigenvector was arbitrarily set too small." << endl; + + // The following "if" is for debugging + if (verbose) { + cout << "Here comes the right vector" << endl; + for (j = 0; j < 4; ++j) { + cout << evec[4*j + i] << " "; + } + cout << endl; + cout << "Here comes the left vector" << endl; + for (j = 0; j < 4; ++j) { + cout << inv_evec[4*i + j] << " "; + } + cout << endl; + } + } + } + + + + for (i = 0; i < 4; ++i) + for (j = 0; j < 4; ++j) evec[4 * i + j] /= scalar_products[j]; + + + checkevector(evec, inv_evec, 4); /* check whether inversion was OK */ + + delete [] square_rate_mat; + delete [] scalar_products; + + // The following is just for debugging + if (verbose) { + + cout << "Roots of scalar products " << endl; + for (int i = 0; i < 4; ++i) cout << scalar_products[i] << " "; + cout << endl; + + + + vector scalars(16); + + for (i = 0; i < 4; ++i){ + for (j = 0; j < 4; ++j) { + scalars[4*i + j] = 0; + for (k = 0; k < 4; ++k) { + scalars[4*i+j] += evec[4*i+k]*evec[4*j+k]; + } + } + } + + cout << "Scalar cross-products " << endl; + for (i = 0; i < 4; ++i) { + for (j = 0; j < 4; ++j) { + cout << scalars[4*i + j]/sqrt(scalars[4*i+i] * scalars[4*j + j]) << " "; + } + cout << endl; + } + + cout << "Right eigenvectors " << endl; + for (i = 0; i < 4; ++i) { + for (j = 0; j < 4; ++j) { + cout << evec[4*i + j] << " "; + } + cout << endl; + } + } + + +} + + +/* Closed formula for the eigenvectors of a REVERSIBLE rate matrix with different eigenvalues */ + +void EigenDecomposition::eigensystem_closed_formula_rev(double *rate_matrix, double *state_freq, double *eval, + double *evec, double *inv_evec) { + + + // Set true to output all interesting parts of the process + bool verbose = false; + + int i,j,k; + + // Some threshold for very small numbers + double small = 1e-6; + + // These bools indicate if we found a nonzero eigenvector; + // only three suffice, because eigenvalue 0 has trivial associated eigenvectors + + vector got_right(3, false); + vector got_left(3, false); + + + auto *square_rate_mat = new double[16]; + memset(square_rate_mat, 0, sizeof(double)*16); + for (i = 0; i < 4; ++i) { + for (j = 0; j < 4; ++j) { + square_rate_mat[4*i + j] = 0.; + for (k = 0; k < 4; ++k) square_rate_mat[4*i + j] += rate_matrix[4*i + k]*rate_matrix[4*k + j]; + } + } + + // Declare the right eigenvectors using close formulae + + + for (i = 0; i < 4; ++i) + evec[4*i + 3] = 1.0; + + for (k = 0; k < 3; ++k) { + int pre = (k + 2) % 3; + int post = (k + 1) % 3; + // the loop for "candidate" iterates looking for a nonzero eigenvector + for (int candidate = 0; candidate < 4 and not got_right[k]; ++ candidate) { + for (i = 0; i < 4; ++i) { + evec[4 * i + k] = eval[pre] * eval[post] * (bool(i == candidate) - state_freq[candidate]) - (eval[pre] + eval[post]) * rate_matrix[4 * i + candidate] + square_rate_mat[4 * i + candidate]; + if (abs(evec[4 * i + k]) > small) got_right[k] = true; + } + } + } + + + + // Declare the left eigenvectors using close formulae + + for (i = 0; i < 4; ++i) + inv_evec[4*3 + i] = state_freq[i]; + + + for (k = 0; k < 3; ++k) { + int pre = (k + 2) % 3; + int post = (k + 1) % 3; + // the loop for "candidate" iterates looking for a nonzero eigenvector + for (int candidate = 0; candidate < 4 and not got_left[k]; ++ candidate) { + for (i = 0; i < 4; ++i) { + inv_evec[4*k + i] = eval[pre]*eval[post]*( bool(i == candidate) - state_freq[i]) - (eval[pre] + eval[post])*rate_matrix[4* candidate + i] + square_rate_mat[4*candidate + i]; + if (abs(inv_evec[4*k + i]) > small) got_left[k] = true; + } + } + } + + + // Rescale the eigenvectors so that inv_evec is the inverse of evec, as also that inv_evec[4*i + j] = evec[4*j + i] * state_freq[j]. + + double small_scalar = 4*small*small; + double scalar_products[4]; + for (i = 0; i < 4; ++i) { + scalar_products[i] = 0.; + for (j = 0; j < 4; ++j) scalar_products[i] += evec[4*j + i] *inv_evec[4*i + j]; + if (abs(scalar_products[i]) < small_scalar) { + cout << "Warning: The module of one eigenvector was arbitrarily set very small." << endl; + + if (verbose) { + cout << "This is the right eigenvector" << endl; + for (j = 0; j < 4; ++j) cout << evec[4*j + i] << " "; + cout << endl; + + cout << "This is the left eigenvector" << endl; + for (j = 0; j < 4; ++j) cout << inv_evec[4*i + j] << " "; + cout << endl; + } + } + } + + double quotients[4]; + for (i = 0; i < 4; ++i) { + quotients[i] = 0; + for (j = 0; j < 4 and quotients[i] == 0; ++j) { + if (abs(evec[4*j + i]) > small) { + quotients[i] = inv_evec[4*i + j]/(evec[4*j + i]*state_freq[j]); + } + } + if (abs(quotients[i]) < small) cout << "Warning: The module of left and right eigenvectors were arbitrarily set very different." << endl; + } + + double scaling_right[4]; + for (i = 0; i < 4; ++i) scaling_right[i] = sqrt(quotients[i]/scalar_products[i]); + + for (i = 0; i < 4; ++i) { + for (j = 0; j < 4; ++j) evec[4 * i + j] *= scaling_right[j]; + } + + double scaling_left[4]; + for (i = 0; i < 4; ++i) scaling_left[i] = sqrt(1/(quotients[i]*scalar_products[i])); + + for (i = 0; i < 4; ++i) { + for (j = 0; j < 4; ++j) inv_evec[4 * i + j] *= scaling_left[i]; + } + + checkevector(evec, inv_evec, 4); /* check whether inversion was OK */ + delete [] square_rate_mat; + + if (verbose) { + double prod[16]; + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < 4; ++j) { + prod[4*i + j] = 0.; + for (int k = 0; k < 4; ++k) prod[4*i + j] += evec[4*i + k] *inv_evec[4*k + j]; + } + } + + double aux[16]; + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < 4; ++j) { + aux[4*i + j] = evec[4*i + j]*eval[j]; + } + } + double recons[16]; + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < 4; ++j) { + recons[4*i + j] = 0.; + for (int k = 0; k < 4; ++k) recons[4*i + j] += aux[4*i + k]*inv_evec[4*k+j]; + } + } + + cout << endl << "I give you the rate matrix: " << endl; + for (int k = 0; k < 4; ++k) { + for (int j = 0; j < 4; ++j) { + cout << rate_matrix[4*k + j] << " "; + } + cout << endl; + } + + cout << endl << "I give you the eigenvalues: " << endl; + for (int k = 0; k < 4; ++k) { + cout << eval[k] << " "; + cout << endl; + } + + cout << endl << "These are my right eigenvectors: " << endl; + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < 4; ++j) cout << evec[4*i + j] << " "; + cout << endl; + } + + cout << endl << "These are my left eigenvectors : " << endl; + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < 4; ++j) cout << inv_evec[4*i + j] << " "; + cout << endl; + } + + cout << endl << "I give you the product between left and right eigenvectors: " << endl; + for (int k = 0; k < 4; ++k) { + for (int j = 0; j < 4; ++j) { + cout << real(prod[4*j + k]) << " "; + } + cout << endl; + } + + cout << endl << "This is the differnce between my aimed reconstruction and the rate matrix: " << endl; + for (int i = 0; i < 4; ++i) { + for (int j = 0; j < 4; ++j) cout << rate_matrix[4*i + j] - recons[4*i + j] << " "; + cout << endl; + } + } +} diff --git a/utils/eigendecomposition.h b/utils/eigendecomposition.h index 314acaa0f..d38dc37a1 100644 --- a/utils/eigendecomposition.h +++ b/utils/eigendecomposition.h @@ -20,6 +20,8 @@ #ifndef EIGENDECOMPOSITION_H #define EIGENDECOMPOSITION_H +#include + //const double ZERO_FREQ = 0.000001; const double ZERO_FREQ = 1e-10; @@ -186,6 +188,77 @@ class EigenDecomposition{ void checkevector(double *evec, double *ivec, int nn); + + /********************************************************* +* Cassius' functions to implement the closed formulas for the eigensystem. +*********************************************************/ + + + /** check that inversion was ok using complex numbers*/ + + void checkevector_complex(std::complex *evec, std::complex *civec, int nn); + + /** + Formula to solve the cubic equation + @param (IN) Solve ax^3 + bx^2 + cx + d = 0 + @param (OUT) root_real real part of eigenvalues + @param (OUT) root_imag imaginary part of eigenvalues + */ + void solve_cubic_equation(double a, double b, double c, double d, double* root_real, double* root_imag); + + /** + Eigenvalues for general 4 DIMENSIONAL non-symmetric matrix using the closed formula of the cubic + @param rate_matrix rate matrix + @param eval (OUT) real part of eigenvalues + @param eval_imag (OUT) imaginary part of eigenvalues + */ + void eigenvalues_closed_formula_nonrev(double *rate_matrix, std::complex *ceval, double *eigenvalues, double *eigenvalues_imag); + + /** + Eigenvalues for general 4 DIMENSIONAL reversible matrix using the closed formula of the cubic + @param rate_matrix rate matrix with all the 16 entries + @param eval (OUT) real part of eigenvalues + @param eval_imag (OUT) imaginary part of eigenvalues + */ + void eigenvalues_closed_formula_rev(double *rate_matrix, double *eval); + + + /** + EigenSystem for 4 DIMENSIONAL general non-symmetric matrix with state frequencies + @param rate_matrix rate matrix with all the 16 entries + @param eval (IN) real part of eigenvalues (computed in eigenvalues_closed_formuka_nonrev() ) + @param eval_imag (IN) imaginary part of eigenvalues (computed in eigenvalues_closed_formuka_nonrev() ) + @param evec (OUT) eigenvectors + @param inv_evec (OUT) inverse matrix of eigenvectors + @param num_state (IN) number of states + */ + void eigensystem_closed_formula_nonrev(double *rate_matrix, double *state_freq, + std::complex *ceval, std::complex *cevec, std::complex *inv_cevec); + + /** + EigenSystem for 4 DIMENSIONAL general non-symmetric matrix with state frequencies + @param rate_matrix rate matrix with all the 16 entries + @param eval (IN) real part of eigenvalues (computed in eigenvalues_closed_formuka_nonrev() ) + @param eval_imag (IN) imaginary part of eigenvalues (computed in eigenvalues_closed_formuka_nonrev() ) + @param evec (OUT) eigenvectors + @param inv_evec (OUT) inverse matrix of eigenvectors + @param num_state (IN) number of states + */ + void eigensystem_closed_formula_nonrev_real(double *rate_matrix, double *state_freq, double *eval, double *evec, double *inv_evec); + + + /** + EigenSystem for 4 DIMENSIONAL general reversible matrix with state frequencies + @param rate_matrix rate matrix with all the 16 entries + @param eval (IN) real eigenvalues (computed in eigenvalues_closed_formula_rev() ) + @param evec (OUT) eigenvectors + @param inv_evec (OUT) inverse matrix of eigenvectors + @param num_state (IN) number of states + */ + void eigensystem_closed_formula_rev(double *rate_matrix, double *state_freq, + double *eval, double *evec, double *inv_evec); + + }; #endif