refactor: Simplify isotopologues handling

src/base/serialiser.h

@@ -7,6 +7,7 @@
 
 #include "templates/keyedVector.h"
 #include "templates/orderedMap.h"
+#include "templates/resolvableKeyedVector.h"
 #include <map>
 #include <vector>
 
@@ -111,7 +112,15 @@ template <typename... Contexts> class Serialisable
         for (const auto &[resolvable, value] : keyedVector)
             group[std::string(resolvable.name())] = value;
         return group;
-    };
+    }
+    template <typename KeyClass, typename ValueClass, typename Lambda>
+    static SerialisedValue fromVectorToTable(const ResolvableKeyedVector<KeyClass, ValueClass> &keyedVector, Lambda getInner)
+    {
+        SerialisedValue group;
+        for (const auto &[resolvable, value] : keyedVector)
+            group[std::string(resolvable.name())] = getInner(value);
+        return group;
+    }
     // A helper function to add elements of a vector to a node under the named heading
     template <typename T, typename Lambda>
     static void fromVectorToTable(const std::vector<T> &vector, std::string name, SerialisedValue &node, Lambda getName)

src/classes/CMakeLists.txt

@@ -34,7 +34,6 @@ add_library(
   histogramSet.cpp
   isotopeMix.cpp
   isotopologue.cpp
-  isotopologues.cpp
   isotopologueSet.cpp
   kVector.cpp
   localMolecule.cpp
@@ -96,7 +95,6 @@ add_library(
   interactionPotential.h
   isotopeMix.h
   isotopologue.h
-  isotopologues.h
   isotopologueSet.h
   kVector.h
   localMolecule.h

src/classes/isotopeMix.cpp

@@ -3,7 +3,7 @@
 
 #include "classes/isotopeMix.h"
 #include "classes/atomType.h"
-#include "classes/isotopologues.h"
+#include "classes/isotopologueSet.h"
 #include "classes/species.h"
 #include "data/elements.h"
 #include "data/isotopes.h"
@@ -43,27 +43,27 @@ double IsotopeMix::population(const AtomType *atomType) const
 double IsotopeMix::fraction(const AtomType *atomType) const { return population(atomType) / totalPopulation_; }
 
 // Create mix from Isotopologues
-void IsotopeMix::create(const std::vector<Isotopologues> &isotopologues,
+void IsotopeMix::create(const std::map<const Species *, double> &speciesPopulations, const IsotopologueSet &isotopologues,
                         const std::vector<std::shared_ptr<AtomType>> &exchangeableTypes)
 {
     mix_.clear();
     totalPopulation_ = 0.0;
     exchangeables_.clear();
 
-    // Loop over Isotopologues and add to the mix
-    for (auto &topos : isotopologues)
+    // Loop over species / populations
+    for (auto &[species, speciesPopulation] : speciesPopulations)
     {
-        // Get normalised Isotopologue populations and species atom type populations
-        auto normalised = topos.normalised();
+        // Get the normalised populations for this species
+        auto topes = isotopologues.normalisedIsotopologues(species);
 
         // Loop over the Isotopologues in the mixture
-        for (const auto &[iso, weight] : normalised)
+        for (const auto &[iso, weight] : topes)
         {
             // Loop over Atoms in the Species, searching for the AtomType/Isotope entry in the isotopes list of the
             // Isotopologue
-            for (const auto &[atomType, atomTypePopulation] : topos.species()->atomTypePopulations())
+            for (const auto &[atomType, atomTypePopulation] : species->atomTypePopulations())
             {
-                auto population = atomTypePopulation * weight * topos.speciesPopulation();
+                auto population = speciesPopulation * atomTypePopulation * weight;
 
                 auto &isotopes = mix_[atomType];
                 if (isotopes.contains(iso->atomTypeIsotope(atomType)))

src/classes/isotopeMix.h

@@ -10,7 +10,7 @@
 
 // Forward Declarations
 class AtomType;
-class Isotopologues;
+class IsotopologueSet;
 
 // IsotopeMix
 class IsotopeMix
@@ -33,7 +33,7 @@ class IsotopeMix
 
     public:
     // Create mix from Isotopologues
-    void create(const std::vector<Isotopologues> &isotopologues,
+    void create(const std::map<const Species *, double> &speciesPopulations, const IsotopologueSet &isotopologues,
                 const std::vector<std::shared_ptr<AtomType>> &exchangeableTypes);
     // Return types/topes mix
     const KeyedVector<const AtomType *, std::map<Sears91::Isotope, double>> &mix() const;

src/classes/isotopologueSet.cpp

@@ -17,109 +17,127 @@ void IsotopologueSet::clear() { isotopologues_.clear(); }
 // Add Isotopologue with the specified relative weight
 void IsotopologueSet::add(const Isotopologue *iso, double relativeWeight)
 {
-    auto it = std::find_if(isotopologues_.begin(), isotopologues_.end(),
-                           [iso](auto &data) { return data.species() == iso->parent(); });
-    if (it != isotopologues_.end())
-        it->mix().add(iso, relativeWeight);
-    else
-    {
-        isotopologues_.emplace_back(iso->parent(), 0);
-        isotopologues_.back().mix().add(iso, relativeWeight);
-    }
+    if (!isotopologues_.contains(iso->parent()))
+        isotopologues_.set(iso->parent(), ResolvableKeyedVector<const Isotopologue *, double>());
+
+    isotopologues_[iso->parent()].add(iso, relativeWeight);
 }
 
 // Remove specified Species from the list (if it exists)
 void IsotopologueSet::remove(const Species *sp)
 {
-    isotopologues_.erase(
-        std::remove_if(isotopologues_.begin(), isotopologues_.end(), [sp](const auto &data) { return data.species() == sp; }),
-        isotopologues_.end());
+    if (!isotopologues_.contains(sp))
+        isotopologues_.erase(sp);
 }
 
 // Remove any occurrences of the specified Isotopologue
 void IsotopologueSet::remove(const Isotopologue *iso)
 {
-    // Get parent Isotopologues from the contained Species pointer
-    auto it = std::find_if(isotopologues_.begin(), isotopologues_.end(),
-                           [iso](auto &data) { return data.species() == iso->parent(); });
-
-    if (it != isotopologues_.end())
+    if (isotopologues_.contains(iso->parent()))
     {
-        it->mix().erase(iso);
+        isotopologues_[iso->parent()].erase(iso);
 
         // Check for Isotopologues now being empty
-        if (it->mix().size() == 0)
-            isotopologues_.erase(it);
+        if (isotopologues_[iso->parent()].size() == 0)
+            isotopologues_.erase(iso->parent());
     }
 }
 
 // Return whether Isotopologues for the specified Species exists
-bool IsotopologueSet::contains(const Species *sp) const
-{
-    return std::any_of(isotopologues_.cbegin(), isotopologues_.cend(),
-                       [sp](const Isotopologues &mix) { return mix.species() == sp; });
-}
+bool IsotopologueSet::contains(const Species *sp) const { return isotopologues_.contains(sp); }
 
-// Return IsotopologueSet for the specified Species
-OptionalReferenceWrapper<const Isotopologues> IsotopologueSet::getIsotopologues(const Species *sp) const
+// Return Isotopologues with normalised populations for the specified Species
+std::map<const Isotopologue *, double> IsotopologueSet::normalisedIsotopologues(const Species *sp) const
 {
-    auto it =
-        std::find_if(isotopologues_.cbegin(), isotopologues_.cend(), [sp](const auto &data) { return data.species() == sp; });
-    if (it == isotopologues_.end())
-        return {};
+    if (isotopologues_.contains(sp))
+    {
+        auto normalised = isotopologues_.value(sp);
+        auto sum = std::accumulate(normalised.begin(), normalised.end(), 0.0,
+                                   [](const auto acc, const auto &isoWeight) { return acc + isoWeight.second; });
+        for (auto &weight : std::views::values(normalised))
+            weight /= sum;
+        return normalised;
+    }
 
-    return *it;
+    std::map<const Isotopologue *, double> natural;
+    natural[sp->naturalIsotopologue()] = 1.0;
+    return natural;
 }
 
 // Return number of species covered by set
 int IsotopologueSet::nSpecies() const { return isotopologues_.size(); }
 
 // Return vector of all Isotopologues
-std::vector<Isotopologues> &IsotopologueSet::isotopologues() { return isotopologues_; }
+ResolvableKeyedVector<const Species *, ResolvableKeyedVector<const Isotopologue *, double>> &IsotopologueSet::isotopologues()
+{
+    return isotopologues_;
+}
 
-const std::vector<Isotopologues> &IsotopologueSet::isotopologues() const { return isotopologues_; }
+const ResolvableKeyedVector<const Species *, ResolvableKeyedVector<const Isotopologue *, double>> &
+IsotopologueSet::isotopologues() const
+{
+    return isotopologues_;
+}
 
 /*
  * Serialisation
  */
 
-// Read data through specified LineParser
-bool IsotopologueSet::deserialise(LineParser &parser, const CoreData &coreData)
+// Express as a serialisable value
+void IsotopologueSet::serialise(std::string tag, SerialisedValue &target) const
 {
-    clear();
+    if (isotopologues_.size() == 0)
+        return;
 
-    const auto nSpecies = parser.argi(0);
-    for (auto n = 0; n < nSpecies; ++n)
-    {
-        // Add a new isotopologue set and read it
-        isotopologues_.emplace_back();
-        if (!isotopologues_.back().deserialise(parser, coreData))
-            return false;
-    }
-
-    return true;
+    SerialisedValue value;
+    value["set"] = fromVectorToTable(isotopologues_, [](const auto &topes)
+                                     { return fromVectorToTable(topes, [](const auto isoWeight) { return isoWeight; }); });
+    target[tag] = value;
 }
 
-// Write data through specified LineParser
-bool IsotopologueSet::write(LineParser &parser)
+// Read values from a serialisable value
+void IsotopologueSet::deserialise(const SerialisedValue &node, const CoreData &coreData)
 {
-    // Write Configuration name and number of Isotopologues we have defined
-    if (!parser.writeLineF("{}\n", isotopologues_.size()))
-        return false;
-
-    // Write details for each set of Isotopologues
-    for (const auto &topes : isotopologues_)
-        if (!topes.serialise(parser))
-            return false;
+    clear();
 
-    return true;
+    toMap(node, "set",
+          [&](const std::string &speciesName, const SerialisedValue &topes)
+          {
+              auto &set = isotopologues_[speciesName];
+              toMap(topes, [&](const std::string &isoName, const SerialisedValue &population)
+                    { set[isoName] = population.as_floating(); });
+          });
+
+    // Need to resolve species
+    std::map<std::string, const Species *> speciesMap;
+    for (const auto &sp : coreData.species())
+        speciesMap[std::string(sp.get()->name())] = sp.get();
+    resolve(speciesMap);
 }
 
-// Express as a serialisable value
-void IsotopologueSet::serialise(std::string tag, SerialisedValue &target) const { target[tag] = isotopologues_; }
-
-// Read values from a serialisable value
-void IsotopologueSet::deserialise(const SerialisedValue &node, const CoreData &coreData)
+// Resolve internal resolvable name references with supplied data
+void IsotopologueSet::resolve(const std::map<std::string, const Species *> &speciesInScope)
 {
-    toVector(node, [this, &coreData](const auto &item) { isotopologues_.emplace_back().deserialise(item, coreData); });
+    for (auto &[resolvableSpecies, topes] : isotopologues_)
+    {
+        // Find the named species
+        if (speciesInScope.contains(std::string(resolvableSpecies.name())))
+            resolvableSpecies.resolve(speciesInScope.at(std::string(resolvableSpecies.name())));
+        else
+            throw(std::runtime_error(std::format("Species '{}' is used in IsotopologueSet, but no such species is in scope.\n",
+                                                 resolvableSpecies.name())));
+
+        // Resolve Isotopologues
+        auto *sp = resolvableSpecies.raw();
+        for (auto &resolvableIsotopologue : std::views::keys(topes))
+        {
+            auto *iso = sp->findIsotopologue(resolvableIsotopologue.name());
+            if (iso)
+                resolvableIsotopologue.resolve(iso);
+            else
+                throw(std::runtime_error(std::format(
+                    "Isotopologue '{}' from species '{}' is used in IsotopologueSet, but no such isotopologue exists.\n",
+                    sp->name(), resolvableIsotopologue.name())));
+        }
+    }
 }

src/classes/isotopologueSet.h

@@ -5,17 +5,14 @@
 
 #include "base/serialiser.h"
 #include "classes/coreData.h"
-#include "classes/isotopologues.h"
-#include "templates/optionalRef.h"
-#include <vector>
 
 // Forward Declarations
 class Species;
 class Isotopologue;
 class LineParser;
 
 // IsotopologueSet - Isotopologues for one or more Species
-class IsotopologueSet : public Serialisable<const CoreData &>
+class IsotopologueSet : public Serialisable<const CoreData &>, ResolvableContext
 {
     public:
     IsotopologueSet() = default;
@@ -28,7 +25,7 @@ class IsotopologueSet : public Serialisable<const CoreData &>
      */
     private:
     // Isotopologue mixtures for individual Species
-    std::vector<Isotopologues> isotopologues_;
+    ResolvableKeyedVector<const Species *, ResolvableKeyedVector<const Isotopologue *, double>> isotopologues_;
 
     public:
     // Clear all existing data
@@ -41,24 +38,22 @@ class IsotopologueSet : public Serialisable<const CoreData &>
     void remove(const Isotopologue *iso);
     // Return whether Isotopologues for the specified Species exists
     bool contains(const Species *sp) const;
-    // Return Isotopologues for the specified Species
-    OptionalReferenceWrapper<const Isotopologues> getIsotopologues(const Species *sp) const;
+    // Return Isotopologues with normalised populations for the specified Species
+    std::map<const Isotopologue *, double> normalisedIsotopologues(const Species *sp) const;
     // Return number of species covered by set
     int nSpecies() const;
     // Return vector of all Isotopologues
-    std::vector<Isotopologues> &isotopologues();
-    const std::vector<Isotopologues> &isotopologues() const;
+    ResolvableKeyedVector<const Species *, ResolvableKeyedVector<const Isotopologue *, double>> &isotopologues();
+    const ResolvableKeyedVector<const Species *, ResolvableKeyedVector<const Isotopologue *, double>> &isotopologues() const;
 
     /*
      * Serialisation
      */
     public:
-    // Read data through specified LineParser
-    bool deserialise(LineParser &parser, const CoreData &coreData);
-    // Write data through specified LineParser
-    bool write(LineParser &parser);
     // Express as a serialisable value
     void serialise(std::string tag, SerialisedValue &target) const override;
     // Read values from a serialisable value
     void deserialise(const SerialisedValue &node, const CoreData &coreData) override;
+    // Resolve internal resolvable name references with supplied data
+    void resolve(const std::map<std::string, const Species *> &speciesInScope) override;
 };