Fix submodel indexing logic in ModelMarkov functions, fix likelihood calculations in AlignmentPairwise

alignment/alignmentpairwise.h

@@ -112,10 +112,22 @@ class AlignmentPairwise : public Alignment, public Optimization
 	*/
     virtual ~AlignmentPairwise();
 
+private:
+    /**
+     *  Implementation of both computeFunction() and computeFuncDerv()
+     *  @param value A distance between the two sequences
+     *  @param[out] lh The log-likelihood at this distance
+     *  @param[out] df The log-likelihood first derivative at this distance
+     *  @param[out] ddf The log-likelihood second derivative at this distance
+     */
+    template <bool COMPUTE_DERV>
+    void likelihoodKernelFunction(double value, double &lh, double &df, double &ddf);
+
+public:
     size_t pairCount;
     size_t derivativeCalculationCount;
     size_t costCalculationCount;
-    
+
 protected:
 	PhyloTree* tree;          //multi-species alignment tree from which sequences
                               //to be aligned are to be drawn
@@ -125,13 +137,12 @@ class AlignmentPairwise : public Alignment, public Optimization
                               //size is num_states_squared times 1 (or by the number
                               //of categories).
     int        trans_size;    //number of elements (rows x columns) in transition matrices
-    double*    trans_mat;     //used in computeFunction(),
-    double*    sum_trans_mat; //used in computeFunction()
+    double*    trans_mat;     //used in computeFunction() and computeFuncDerv()
     double*    trans_derv1;   //used in computeFuncDerv()
     double*    trans_derv2;   //used in computeFuncDerv()
+    double*    sum_trans;     //used in computeFunction() and computeFuncDerv()
     double*    sum_derv1;     //used in computeFuncDerv()
     double*    sum_derv2;     //used in computeFuncDerv()
-    double*    sum_trans;     //used in computeFuncDerv()
 
     int        seq_id1;
     int        seq_id2;

model/modelcodon.cpp

@@ -832,7 +832,7 @@ void ModelCodon::decomposeRateMatrix() {
     ModelMarkov::decomposeRateMatrix();
 }
 
-void ModelCodon::getQMatrix(double *q_mat, int mixture) {
+void ModelCodon::getQMatrix(double *q_mat, int) {
 
     double **rate_matrix = (double**) new double[num_states];
     int i, j;

model/modelcodon.h

@@ -100,7 +100,7 @@ class ModelCodon: public ModelMarkov {
      * compute Q matrix
      * @param q_mat (OUT) Q matrix, assuming of size num_states * num_states
      */
-    virtual void getQMatrix(double *q_mat, int mixture = 0);
+    virtual void getQMatrix(double *q_mat, int model_id = -1);
 
 	/**
 	 * read codon model from a stream, modying rates and state_freq accordingly

model/modelfactory.cpp

@@ -1771,17 +1771,18 @@ void ModelFactory::stopStoringTransMatrix() {
 }
 
 
-double ModelFactory::computeTrans(double time, int state1, int state2) {
-    return model->computeTrans(time, state1, state2);
+double ModelFactory::computeTrans(double time, int state1, int state2, int model_id) {
+    return model->computeTrans(time, state1, state2, model_id);
 }
 
-double ModelFactory::computeTrans(double time, int state1, int state2, double &derv1, double &derv2) {
-    return model->computeTrans(time, state1, state2, derv1, derv2);
+double ModelFactory::computeTrans(double time, int state1, int state2,
+                                  double &derv1, double &derv2, int model_id) {
+    return model->computeTrans(time, state1, state2, derv1, derv2, model_id);
 }
 
-void ModelFactory::computeTransMatrix(double time, double *trans_matrix, int mixture, int selected_row) {
+void ModelFactory::computeTransMatrix(double time, double *trans_matrix, int model_id, int selected_row) {
     if (!store_trans_matrix || !is_storing || model->isSiteSpecificModel()) {
-        model->computeTransMatrix(time, trans_matrix, mixture, selected_row);
+        model->computeTransMatrix(time, trans_matrix, model_id, selected_row);
         return;
     }
     int mat_size = model->num_states * model->num_states;
@@ -1790,20 +1791,19 @@ void ModelFactory::computeTransMatrix(double time, double *trans_matrix, int mix
         // allocate memory for 3 matricies
         double *trans_entry = new double[mat_size * 3];
         trans_entry[mat_size] = trans_entry[mat_size+1] = 0.0;
-        model->computeTransMatrix(time, trans_entry, mixture, selected_row);
+        model->computeTransMatrix(time, trans_entry, model_id, selected_row);
         ass_it = insert(value_type(round(time * 1e6), trans_entry)).first;
     } else {
         //if (verbose_mode >= VB_MAX)
             //cout << "ModelFactory bingo" << endl;
     }
-
     memcpy(trans_matrix, ass_it->second, mat_size * sizeof(double));
 }
 
 void ModelFactory::computeTransDerv(double time, double *trans_matrix,
-    double *trans_derv1, double *trans_derv2, int mixture) {
+                                    double *trans_derv1, double *trans_derv2, int model_id) {
     if (!store_trans_matrix || !is_storing || model->isSiteSpecificModel()) {
-        model->computeTransDerv(time, trans_matrix, trans_derv1, trans_derv2, mixture);
+        model->computeTransDerv(time, trans_matrix, trans_derv1, trans_derv2, model_id);
         return;
     }
     int mat_size = model->num_states * model->num_states;
@@ -1812,11 +1812,11 @@ void ModelFactory::computeTransDerv(double time, double *trans_matrix,
         // allocate memory for 3 matricies
         double *trans_entry = new double[mat_size * 3];
         trans_entry[mat_size] = trans_entry[mat_size+1] = 0.0;
-        model->computeTransDerv(time, trans_entry, trans_entry+mat_size, trans_entry+(mat_size*2), mixture);
+        model->computeTransDerv(time, trans_entry, trans_entry+mat_size, trans_entry+(mat_size*2), model_id);
         ass_it = insert(value_type(round(time * 1e6), trans_entry)).first;
     } else if (ass_it->second[mat_size] == 0.0 && ass_it->second[mat_size+1] == 0.0) {
         double *trans_entry = ass_it->second;
-        model->computeTransDerv(time, trans_entry, trans_entry+mat_size, trans_entry+(mat_size*2), mixture);
+        model->computeTransDerv(time, trans_entry, trans_entry+mat_size, trans_entry+(mat_size*2), model_id);
     }
     memcpy(trans_matrix, ass_it->second, mat_size * sizeof(double));
     memcpy(trans_derv1, ass_it->second + mat_size, mat_size * sizeof(double));

model/modelfactory.h

@@ -139,47 +139,29 @@ class ModelFactory : public unordered_map<int, double*>, public Optimization, pu
 	*/
 	void stopStoringTransMatrix();
 
-	/**
-		Wrapper for computing the transition probability matrix from the model. It use ModelFactory
-		that stores matrix computed before for effiency purpose.
-		@param time time between two events
-        @param mixture (optional) class for mixture model
-        @param selected_row (optional) only compute the entries of one selected row. By default, compute all rows
-		@param trans_matrix (OUT) the transition matrix between all pairs of states.
-			Assume trans_matrix has size of num_states * num_states.
-	*/
-	void computeTransMatrix(double time, double *trans_matrix, int mixture = 0, int selected_row = -1);
+    /**
+     *  Wrapper for the computeTransMatrix() function of the model.
+     *  ModelFactory stores the matrix computed before for efficiency purposes
+     */
+    void computeTransMatrix(double time, double *trans_matrix, int model_id = -1, int selected_row = -1);
 
-	/**
-		Wrapper for computing the transition probability between two states.
-		@param time time between two events
-		@param state1 first state
-		@param state2 second state
-	*/
-	double computeTrans(double time, int state1, int state2);
+    /**
+     *  Wrapper for the computeTransDerv() function of the model.
+     *  ModelFactory stores the matrix computed before for efficiency purposes
+     */
+    void computeTransDerv(double time, double *trans_matrix,
+                          double *trans_derv1, double *trans_derv2, int model_id = -1);
 
-	/**
-		Wrapper for computing the transition probability between two states
-		@param time time between two events
-		@param state1 first state
-		@param state2 second state
-		@param derv1 (OUT) 1st derivative
-		@param derv2 (OUT) 2nd derivative
-	*/
-	virtual double computeTrans(double time, int state1, int state2, double &derv1, double &derv2);
+    /**
+     *  Wrapper for the computeTrans() function of the model
+     */
+    double computeTrans(double time, int state1, int state2, int model_id = -1);
 
-	/**
-		Wrapper for computing the transition probability matrix and the derivative 1 and 2 from the model.
-		It use ModelFactory that stores matrix computed before for effiency purpose.
-		@param time time between two events
-        @param mixture (optional) class for mixture model
-		@param trans_matrix (OUT) the transition matrix between all pairs of states. 
-			Assume trans_matrix has size of num_states * num_states.
-		@param trans_derv1 (OUT) the 1st derivative matrix between all pairs of states. 
-		@param trans_derv2 (OUT) the 2nd derivative matrix between all pairs of states. 
-	*/
-	void computeTransDerv(double time, double *trans_matrix, 
-		double *trans_derv1, double *trans_derv2, int mixture = 0);
+    /**
+     *  Wrapper for the computeTrans() function of the model
+     */
+    double computeTrans(double time, int state1, int state2,
+                        double &derv1, double &derv2, int model_id = -1);
 
 	/**
 		 destructor

model/modelliemarkov.cpp

@@ -2269,9 +2269,11 @@ void ModelLieMarkov::decomposeRateMatrixClosedForm() {
     }
 }
 
+// OBSOLETE
+/*
 void ModelLieMarkov::computeTransMatrix(double time, double *trans_matrix, int mixture, int selected_row) {
     return ModelMarkov::computeTransMatrix(time, trans_matrix, mixture, selected_row);
-    /*
+
   MatrixExpTechnique technique = phylo_tree->params->matrix_exp_technique;
   if (technique == MET_SCALING_SQUARING || nondiagonalizable ) {
         Matrix4d A = Map<Matrix4d>(rate_matrix);
@@ -2331,6 +2333,5 @@ void ModelLieMarkov::computeTransMatrix(double time, double *trans_matrix, int m
 
     } else
         ModelMarkov::computeTransMatrix(time, trans_matrix);
-     */
 }
-
+*/

model/modelliemarkov.h

@@ -95,15 +95,6 @@ class ModelLieMarkov: public ModelMarkov {
     /** decompose rate matrix using Eigen library */
     virtual void decomposeRateMatrixEigen3lib();
 
-	/**
-		compute the transition probability matrix.
-		@param time time between two events
-        @param mixture (optional) class for mixture model
-        @param selected_row (optional) only compute the entries of one selected row. By default, compute all rows
-		@param trans_matrix (OUT) the transition matrix between all pairs of states.
-			Assume trans_matrix has size of num_states * num_states.
-	*/
-	virtual void computeTransMatrix(double time, double *trans_matrix, int mixture = 0, int selected_row = -1);
 	// overrides Optimization::restartParameters
 	bool restartParameters(double guess[], int ndim, double lower[], double upper[], bool bound_check[], int iteration);