Add tests for rate variation

Basic tests for rate variation.

Author: SR-Martin

model/model_dna_rate_variation.cpp

@@ -535,4 +535,25 @@ void ModelDNARateVariation::updateCountsAndWaitingTimesAcrossRoot( PositionType
             waitingTimes[i][childState] += weight * distToRoot;
         }      
     }
+}
+
+void ModelDNARateVariation::setAllMatricesToDefault() {
+    for(int i = 0; i < genomeSize; i++) {
+        for(int stateA = 0; stateA < num_states_; stateA++) {
+            RealNumType rowSum = 0;
+            for(int stateB = 0; stateB < num_states_; stateB++) {
+                mutationMatrices[i * matSize + (stateB + row_index[stateA])] = mutation_mat[stateB + row_index[stateA]];
+                transposedMutationMatrices[i * matSize + (stateB + row_index[stateA])] = transposed_mut_mat[stateB + row_index[stateA]];
+                freqiFreqjQijs[i * matSize + (stateB + row_index[stateA])] = freqi_freqj_qij[stateB + row_index[stateA]];
+            }
+
+            diagonalMutationMatrices[i * num_states_ + stateA] = diagonal_mut_mat[stateA];
+
+            // pre-compute matrix to speedup
+            const RealNumType* transposed_mut_mat_row = getTransposedMutationMatrixRow(stateA, i);
+            RealNumType* freqjTransposedijsRow = freqjTransposedijs + (i * matSize) + row_index[stateA];
+            setVecByProduct<4>(freqjTransposedijsRow, root_freqs, transposed_mut_mat_row);
+        
+        }
+    }  
 }

model/model_dna_rate_variation.h

@@ -26,6 +26,7 @@ class ModelDNARateVariation : public ModelDNA {
     }; 
 
     virtual inline const cmaple::RealNumType *const getMutationMatrixRow(StateType row, PositionType i) const override {
+        assert(i < genomeSize);
         return mutationMatrices + (i * matSize) + row_index[row];
     }; 
 
@@ -70,6 +71,8 @@ class ModelDNARateVariation : public ModelDNA {
    */
   virtual bool updateMutationMatEmpirical() override;
 
+  void setAllMatricesToDefault();
+
   void printMatrix(const RealNumType* matrix, std::ostream* outStream);
   void printCountsAndWaitingTimes(const RealNumType* counts, const RealNumType* waitingTImes, std::ostream* outStream);
 

unittest/CMakeLists.txt

@@ -10,6 +10,7 @@ add_executable(
   model_test.cpp
   alignment_test.cpp
   matrix_test.cpp
+  ratevariation_test.cpp
   sequence_test.cpp
   seqregion_test.cpp
   mutation_test.cpp