Fix Coala IO

Author: lgueguen

src/Bpp/Phyl/Io/BppOFrequencySetFormat.cpp

@@ -530,12 +530,12 @@ std::unique_ptr<FrequencySetInterface> BppOFrequencySetFormat::readFrequencySet(
     if (dynamic_cast<Coala*>(model.get()) == 0)
       throw Exception("MVAprotein frequency set needs a Coala model");
 
-    auto coala = unique_ptr<Coala>(dynamic_cast<Coala*>(model.release()));
+    auto coala = shared_ptr<Coala>(dynamic_cast<Coala*>(model.release()));
     // map<string, string> unparsedParameterValuesNested(nestedReader.getUnparsedArguments());
 
     auto mvaFS = make_unique<MvaFrequencySet>(dynamic_pointer_cast<const ProteicAlphabet>(alphabet));
     // mvaFS->setParamValues(args);
-    mvaFS->initSet(*coala);
+    mvaFS->initSet(coala);
 
     pFS = std::move(mvaFS);
   }
@@ -724,6 +724,24 @@ void BppOFrequencySetFormat::writeFrequencySet(
   catch (bad_cast&)
   {}
 
+  // MV
+
+  try
+  {
+    auto pmvaFS = dynamic_cast<const MvaFrequencySet&>(freqset);
+    if (comma)
+      out << ", ";
+    out << "model=";
+
+    BppOSubstitutionModelFormat bIO(BppOSubstitutionModelFormat::ALL, true, true, true, false, 0);
+
+    bIO.write(*(pmvaFS.getModel()), out, globalAliases, writtenNames);
+    comma = true;
+  }
+  catch (bad_cast&)
+  {}
+
+
 
   // FullPerAA
   try

src/Bpp/Phyl/Io/BppOSubstitutionModelFormat.cpp

@@ -856,9 +856,15 @@ unique_ptr<SubstitutionModelInterface> BppOSubstitutionModelFormat::readSubstitu
         if (TextTools::isEmpty(nbrOfParametersPerBranch))
           throw Exception("'nbrAxes' argument missing to define the number of axis of the Correspondence Analysis.");
         // Now we create the Coala model
-        model = make_unique<Coala>(alpha, *nestedModel, TextTools::to<unsigned int>(nbrOfParametersPerBranch));
-        if (nData)
-          model->setFreqFromData(*mData.at(nData));
+        auto coala = make_unique<Coala>(alpha, *nestedModel, TextTools::to<unsigned int>(nbrOfParametersPerBranch));
+        if (!nData)
+          throw Exception("'data' argument missing to build the Correspondence Analysis.");
+        else
+        {
+          coala->setNData(nData);
+          coala->setFreqFromData(*mData.at(nData));
+        }
+        model = unique_ptr<SubstitutionModelInterface>(dynamic_cast<SubstitutionModelInterface*>(coala.release()));
       }
       else
         throw Exception("Model '" + modelName + "' is unknown, or does not fit proteic alphabet.");
@@ -1698,6 +1704,8 @@ void BppOSubstitutionModelFormat::write(const BranchModelInterface& model,
   {
     const Coala& coalaModel = dynamic_cast<const Coala&>(model);
     out << "exch=" << coalaModel.getExch() << ", nbrAxes=" << coalaModel.getNbrOfAxes();
+    if (coalaModel.getNData() > 0)
+      out << ", data=" << coalaModel.getNData();
     comma = true;
   }
   catch (bad_cast&)

src/Bpp/Phyl/Legacy/App/PhylogeneticsApplicationTools.cpp

@@ -543,7 +543,7 @@ void LegacyPhylogeneticsApplicationTools::setSubstitutionModelSet(
     {
       auto tmp2 = dynamic_pointer_cast<Coala>(tmp);
       if (tmp2)
-        dynamic_pointer_cast<MvaFrequencySet>(rootFrequencies)->initSet(*tmp2);
+        dynamic_pointer_cast<MvaFrequencySet>(rootFrequencies)->initSet(tmp2);
       else
         throw Exception("The MVAprotein frequencies set at the root can only be used if a Coala model is used on branches.");
     }

src/Bpp/Phyl/Model/FrequencySet/MvaFrequencySet.cpp

@@ -14,20 +14,21 @@ MvaFrequencySet::MvaFrequencySet(shared_ptr<const ProteicAlphabet> alpha) :
       make_shared<CanonicalStateMap>(alpha, false),
       "MVA.",
       "MVAprotein"),
+  model_(),
   tPpalAxes_(),
   rowCoords_(),
   nbrOfAxes_(0),
-  model_(),
   columnWeights_(),
   paramValues_()
 {}
 
-void MvaFrequencySet::initSet(const CoalaCore& coala)
+void MvaFrequencySet::initSet(std::shared_ptr<const Coala> coala)
 {
-  setNbrOfAxes(coala.getNbrOfAxes());
-  setTransposeMatrixOfPpalAxes(coala.getTppalAxesMatrix());
-  setMatrixOfRowCoords(coala.getRowCoordinates());
-  setVectorOfColumnWeights(coala.getColumnWeights());
+  model_=coala;
+  setNbrOfAxes(coala->getNbrOfAxes());
+  setTransposeMatrixOfPpalAxes(coala->getTppalAxesMatrix());
+  setMatrixOfRowCoords(coala->getRowCoordinates());
+  setVectorOfColumnWeights(coala->getColumnWeights());
   defineParameters();
   updateFrequencies();
 }

src/Bpp/Phyl/Model/FrequencySet/MvaFrequencySet.h

@@ -22,6 +22,9 @@ class MvaFrequencySet :
   public virtual ProteinFrequencySetInterface,
   public AbstractFrequencySet
 {
+private:
+  std::shared_ptr<const Coala> model_;
+
 public:
   /**
    * @brief Constructor
@@ -48,7 +51,6 @@ class MvaFrequencySet :
   RowMatrix<double> tPpalAxes_;
   RowMatrix<double> rowCoords_;
   size_t nbrOfAxes_;
-  std::string model_;
   std::vector<double> columnWeights_;
   std::map<std::string, std::string> paramValues_;
 
@@ -58,10 +60,20 @@ class MvaFrequencySet :
     return std::dynamic_pointer_cast<const ProteicAlphabet>(getAlphabet());
   }
 
+  const Coala& model() const
+  {
+    return *model_;
+  }
+
+  std::shared_ptr<const Coala> getModel() const
+  {
+    return model_;
+  }
+
   void setTransposeMatrixOfPpalAxes(const RowMatrix<double>& matrix) { tPpalAxes_ = matrix; }
   void setMatrixOfRowCoords(const RowMatrix<double>& matrix) { rowCoords_ = matrix; }
   void setNbrOfAxes(const size_t& nAxes) { nbrOfAxes_ = nAxes; }
-  void setModelName(const std::string& modelName) { model_ = modelName; }
+//  void setModelName(const std::string& modelName) { model_ = modelName; }
   void setVectorOfColumnWeights(const std::vector<double>& cw) { columnWeights_ = cw; }
   void setParamValues(std::map<std::string, std::string>& valuesSettings) {paramValues_ = valuesSettings;}
 
@@ -71,7 +83,7 @@ class MvaFrequencySet :
   void fireParameterChanged(const ParameterList& parameters) override;
   void updateFrequencies();
 
-  void initSet(const CoalaCore& coala);
+  void initSet(std::shared_ptr<const Coala> coala);
 
   void computeReversePCA(const std::vector<double>& positions, std::vector<double>& tmpFreqs);
   void computeCoordsFirstSpacePCA(std::vector<double>& tmpFreqs, std::vector<double>& freqs);

src/Bpp/Phyl/Model/Protein/CoalaCore.h

@@ -38,6 +38,7 @@ class CoalaCore
   std::vector<double> colWeights_;
   ParameterList paramValues_;
 
+  size_t nData_; // number of the used data (for output consistency)
 public:
   CoalaCore(size_t nbAxes = 0);
 
@@ -51,6 +52,16 @@ class CoalaCore
   const RowMatrix<double>& getRowCoordinates() const { return R_; }
   const std::vector<double>& getColumnWeights() const { return colWeights_; }
 
+  size_t getNData() const
+  {
+    return nData_;
+  }
+
+  void setNData(size_t n)
+  {
+    nData_=n;
+  }
+  
 protected:
   /**
    * @brief Comoute COA from data and return the axis position parameters