Should now be able to crank up photon counts without running out of m…

…emory (as fast). References weren't being used where they should've been

include/MIDSX/Core/quantity.h

@@ -10,12 +10,12 @@
 class VectorValue {
 public:
     // overloading the + operator
-    VectorValue operator+(const VectorValue& other) const;
+    VectorValue operator+(VectorValue& other) const;
 
     // pretty self-explanatory
     void addValue(double value);
     void addValues(const std::vector<double>& values);
-    std::vector<double> getVector() const;
+    std::vector<double>& getVector();
     double getSum();
     double getSumSTD();
     double getMean();

include/MIDSX/Core/surface_quantity.h

@@ -61,7 +61,7 @@ class VectorSurfaceQuantity {
      * @return The sum of the two VectorSurfaceQuantity objects.
      */
 
-    VectorSurfaceQuantity operator+(const VectorSurfaceQuantity& other) const;
+    VectorSurfaceQuantity operator+(VectorSurfaceQuantity& other) const;
 
     /**
      * @brief Measures the VectorSurfaceQuantity for a TempSurfaceTallyData object.
@@ -78,11 +78,11 @@ class VectorSurfaceQuantity {
     VectorSurfaceQuantityType getType() const;
 
     // Wrapper functions for VectorValue. Self-explanatory.
-    VectorValue getTotalValues();
-    VectorValue getPrimaryValues() const;
-    VectorValue getSingleIncoherentScatterValues() const;
-    VectorValue getSingleCoherentScatterValues() const;
-    VectorValue getMultipleScatterValues() const;
+    VectorValue& getTotalValues();
+    VectorValue& getPrimaryValues();
+    VectorValue& getSingleIncoherentScatterValues();
+    VectorValue& getSingleCoherentScatterValues();
+    VectorValue& getMultipleScatterValues();
 private:
     VectorSurfaceQuantityType type_;
     bool totaled_ = false;

include/MIDSX/Core/surface_quantity_container.h

@@ -40,7 +40,13 @@ class SurfaceQuantityContainer {
      * @param other The SurfaceQuantityContainer to add to this SurfaceQuantityContainer.
      * @return The sum of the two SurfaceQuantityContainer objects.
      */
-    SurfaceQuantityContainer operator+(const SurfaceQuantityContainer& other) const;
+    SurfaceQuantityContainer operator+(SurfaceQuantityContainer& other) const;
+
+    /**
+     * @brief Clears the container of all quantities.
+     *
+     */
+    void clear();
 
     /**
      * @brief Gets the map of vector quantities in the container.

include/MIDSX/Core/volume_quantity.h

@@ -59,7 +59,7 @@ class VectorVolumeQuantity {
      * @param other The VectorVolumeQuantity to add to this VectorVolumeQuantity.
      * @return The sum of the two VectorVolumeQuantity objects.
      */
-    VectorVolumeQuantity operator+(const VectorVolumeQuantity& other) const;
+    VectorVolumeQuantity operator+(VectorVolumeQuantity& other) const;
 
     /**
      * @brief Measures the VectorVolumeQuantity for a TempVolumeTallyData object.
@@ -74,11 +74,11 @@ class VectorVolumeQuantity {
      * @return The type of the VectorVolumeQuantity.
      */
     VectorVolumeQuantityType getType() const;
-    VectorValue getTotalValues();
-    VectorValue getPrimaryValues() const;
-    VectorValue getSingleIncoherentScatterValues() const;
-    VectorValue getSingleCoherentScatterValues() const;
-    VectorValue getMultipleScatterValues() const;
+    VectorValue& getTotalValues();
+    VectorValue& getPrimaryValues();
+    VectorValue& getSingleIncoherentScatterValues();
+    VectorValue& getSingleCoherentScatterValues();
+    VectorValue& getMultipleScatterValues();
 private:
     VectorVolumeQuantityType type_;
     bool totaled_ = false;

include/MIDSX/Core/volume_quantity_container.h

@@ -41,7 +41,13 @@ class VolumeQuantityContainer {
      * @param other The VolumeQuantityContainer to add to this VolumeQuantityContainer.
      * @return The sum of the two VolumeQuantityContainer objects.
      */
-    VolumeQuantityContainer operator+(const VolumeQuantityContainer& other) const;
+    VolumeQuantityContainer operator+(VolumeQuantityContainer& other) const;
+
+    /**
+     * @brief Clears the container of all quantities.
+     *
+     */
+    void clear();
 
     /**
      * @brief Gets the map of vector quantities in the container.

src/MIDSX/Core/derived_quantities.cpp

@@ -3,27 +3,30 @@
 double DerivedQuantity::getPrimaryFluence(SurfaceQuantityContainer &surface_quantity_container, double energy, double energy_width, bool is_cosine_weighted) {
     // fluence (E +- delta E / 2) = dN/dE = total number of photons in energy range / area of surface  [units: 1/cm^2]
     double total = 0;
-    VectorValue incident_energies;
 
     try {
-        incident_energies = surface_quantity_container.getVectorQuantities().at(VectorSurfaceQuantityType::IncidentEnergy).getPrimaryValues();
+        auto& incident_energies = surface_quantity_container.getVectorQuantities().at(VectorSurfaceQuantityType::IncidentEnergy);
     } catch (std::out_of_range& e) {
         throw std::runtime_error("SurfaceQuantityContainer does not contain incident energy vector quantity required for cosine-weighted fluence calculation");
     }
+
+    VectorValue& incident_energies = surface_quantity_container.getVectorQuantities().at(VectorSurfaceQuantityType::IncidentEnergy).getPrimaryValues();
+
     // get vector
-    std::vector<double> incident_energies_vector;
-    incident_energies_vector = incident_energies.getVector();
+    std::vector<double>& incident_energies_vector = incident_energies.getVector();
 
     if (is_cosine_weighted) {
         // sum of entrance cosine
-        VectorValue entrance_cosines;
         try {
-            entrance_cosines = surface_quantity_container.getVectorQuantities().at(VectorSurfaceQuantityType::EntranceCosine).getPrimaryValues();
+            auto& entrance_cosines = surface_quantity_container.getVectorQuantities().at(VectorSurfaceQuantityType::EntranceCosine);
         } catch (std::out_of_range& e) {
             throw std::runtime_error("SurfaceQuantityContainer does not contain entrance cosine vector quantity required for cosine-weighted fluence calculation");
         }
+
+        VectorValue& entrance_cosines = surface_quantity_container.getVectorQuantities().at(VectorSurfaceQuantityType::EntranceCosine).getPrimaryValues();
+
         // get vectors
-        std::vector<double> entrance_cosines_vector = entrance_cosines.getVector();
+        std::vector<double>& entrance_cosines_vector = entrance_cosines.getVector();
 
         for (int i = 0; i < incident_energies_vector.size(); i++) {
             if (incident_energies_vector[i] >= energy - energy_width/2 && incident_energies_vector[i] <= energy + energy_width/2) {
@@ -45,7 +48,7 @@ double DerivedQuantity::getPrimaryAirKerma(SurfaceQuantityContainer &surface_qua
                                            double energy, double energy_width, bool is_cosine_weighted) {
     // AK(E +- delta E / 2) = E * fluence(E +- delta E / 2) * (mu_en(E)/rho)  [units: eV/g]
     int AIR_INDEX = 3;
-    auto air_data = interaction_data.getMaterialFromId(AIR_INDEX).getData();
+    auto& air_data = interaction_data.getMaterialFromId(AIR_INDEX).getData();
     double mu_en = air_data.interpolateMassEnergyAbsorptionCoefficient(energy); // [units: cm^2/g]
     double fluence = getPrimaryFluence(surface_quantity_container, energy, energy_width, is_cosine_weighted); // [units: 1/cm^2]
     return energy * fluence * mu_en; // [units: eV/g]

src/MIDSX/Core/physics_engine.cpp

@@ -129,8 +129,10 @@ std::vector<VolumeQuantityContainer> PhysicsEngine::getVolumeQuantityContainers(
 
     for (int i = 0; i < num_of_threads; ++i) {
         for (int j = 0; j < num_of_volume_tallies; ++j) {
-            auto thread_local_volume_container = thread_local_volume_tallies_[i][j]->getVolumeQuantityContainer();
+            auto& thread_local_volume_container = thread_local_volume_tallies_[i][j]->getVolumeQuantityContainer();
             combined_volume_containers[j] = combined_volume_containers[j] + thread_local_volume_container;
+            // clear the thread local container to save memory
+            thread_local_volume_container.clear();
         }
     }
 
@@ -147,6 +149,8 @@ std::vector<SurfaceQuantityContainer> PhysicsEngine::getSurfaceQuantityContainer
         for (int j = 0; j < num_of_surface_tallies; ++j) {
             auto thread_local_surface_container = thread_local_surface_tallies_[i][j]->getSurfaceQuantityContainer();
             combined_surface_containers[j] = combined_surface_containers[j] + thread_local_surface_container;
+            // clear the thread local container to save memory
+            thread_local_surface_container.clear();
         }
     }
 
@@ -158,6 +162,8 @@ void PhysicsEngine::addVoxelDataToComputationalDomain() {
         for (auto &temp_voxel_data: thread_local_voxel_data) {
             temp_voxel_data.voxel.dose.addValue(temp_voxel_data.energy_deposited);
         }
+        // after adding all the dose values, clear the thread local voxel data to save memory
+        thread_local_voxel_data.clear();
     }
 }
 

src/MIDSX/Core/quantity.cpp

@@ -1,6 +1,6 @@
 #include "Core/quantity.h"
 
-VectorValue VectorValue::operator+(const VectorValue& other) const {
+VectorValue VectorValue::operator+(VectorValue& other) const {
     // combine two vectors
     VectorValue sum = *this;
     sum.addValues(other.getVector());
@@ -20,7 +20,7 @@ void VectorValue::addValues(const std::vector<double>& values) {
     values_.insert(values_.end(), values.begin(), values.end());
 }
 
-std::vector<double> VectorValue::getVector() const {
+std::vector<double>& VectorValue::getVector() {
     return values_;
 }
 

src/MIDSX/Core/surface_quantity.cpp

@@ -34,7 +34,7 @@ VectorSurfaceQuantity::VectorSurfaceQuantity(VectorSurfaceQuantityType type) {
     }
 }
 
-VectorSurfaceQuantity VectorSurfaceQuantity::operator+(const VectorSurfaceQuantity& other) const {
+VectorSurfaceQuantity VectorSurfaceQuantity::operator+(VectorSurfaceQuantity& other) const {
     if (type_ != other.type_) {
         throw std::runtime_error("Cannot add VectorSurfaceQuantity objects of different types");
     }
@@ -66,7 +66,7 @@ VectorSurfaceQuantityType VectorSurfaceQuantity::getType() const {
     return type_;
 }
 
-VectorValue VectorSurfaceQuantity::getTotalValues() {
+VectorValue& VectorSurfaceQuantity::getTotalValues() {
     if (!totaled_) {
         total_values_.addValues(primary_values_.getVector());
         total_values_.addValues(single_incoherent_scatter_values_.getVector());
@@ -77,19 +77,19 @@ VectorValue VectorSurfaceQuantity::getTotalValues() {
     return total_values_;
 }
 
-VectorValue VectorSurfaceQuantity::getPrimaryValues() const {
+VectorValue& VectorSurfaceQuantity::getPrimaryValues() {
     return primary_values_;
 }
 
-VectorValue VectorSurfaceQuantity::getSingleIncoherentScatterValues() const {
+VectorValue& VectorSurfaceQuantity::getSingleIncoherentScatterValues() {
     return single_incoherent_scatter_values_;
 }
 
-VectorValue VectorSurfaceQuantity::getSingleCoherentScatterValues() const {
+VectorValue& VectorSurfaceQuantity::getSingleCoherentScatterValues() {
     return single_coherent_scatter_values_;
 }
 
-VectorValue VectorSurfaceQuantity::getMultipleScatterValues() const {
+VectorValue& VectorSurfaceQuantity::getMultipleScatterValues() {
     return multiple_scatter_values_;
 }
 

src/MIDSX/Core/surface_quantity_container.cpp

@@ -17,7 +17,7 @@ void SurfaceQuantityContainer::measureAll(TempSurfaceTallyData& temp_tally_data)
     }
 }
 
-SurfaceQuantityContainer SurfaceQuantityContainer::operator+(const SurfaceQuantityContainer& other) const {
+SurfaceQuantityContainer SurfaceQuantityContainer::operator+(SurfaceQuantityContainer& other) const {
     // check if either container is empty, if so return the other container
     if (vector_quantities_.empty()) {
         return other;
@@ -53,6 +53,11 @@ SurfaceQuantityContainer SurfaceQuantityContainer::operator+(const SurfaceQuanti
     return new_container;
 }
 
+void SurfaceQuantityContainer::clear() {
+    vector_quantities_.clear();
+    count_quantities_.clear();
+}
+
 std::unordered_map<VectorSurfaceQuantityType, VectorSurfaceQuantity>& SurfaceQuantityContainer::getVectorQuantities() {
     return vector_quantities_;
 }

src/MIDSX/Core/volume_quantity.cpp

@@ -38,7 +38,7 @@ VectorVolumeQuantity::VectorVolumeQuantity(VectorVolumeQuantityType type) {
     }
 }
 
-VectorVolumeQuantity VectorVolumeQuantity::operator+(const VectorVolumeQuantity& other) const {
+VectorVolumeQuantity VectorVolumeQuantity::operator+(VectorVolumeQuantity& other) const {
     if (type_ != other.type_) {
         throw std::runtime_error("Cannot add VectorVolumeQuantity objects of different types");
     }
@@ -70,7 +70,7 @@ VectorVolumeQuantityType VectorVolumeQuantity::getType() const {
     return type_;
 }
 
-VectorValue VectorVolumeQuantity::getTotalValues() {
+VectorValue& VectorVolumeQuantity::getTotalValues() {
     if (!totaled_) {
         total_values_.addValues(primary_values_.getVector());
         total_values_.addValues(single_incoherent_scatter_values_.getVector());
@@ -81,19 +81,19 @@ VectorValue VectorVolumeQuantity::getTotalValues() {
     return total_values_;
 }
 
-VectorValue VectorVolumeQuantity::getPrimaryValues() const {
+VectorValue& VectorVolumeQuantity::getPrimaryValues() {
     return primary_values_;
 }
 
-VectorValue VectorVolumeQuantity::getSingleIncoherentScatterValues() const {
+VectorValue& VectorVolumeQuantity::getSingleIncoherentScatterValues() {
     return single_incoherent_scatter_values_;
 }
 
-VectorValue VectorVolumeQuantity::getSingleCoherentScatterValues() const {
+VectorValue& VectorVolumeQuantity::getSingleCoherentScatterValues() {
     return single_coherent_scatter_values_;
 }
 
-VectorValue VectorVolumeQuantity::getMultipleScatterValues() const {
+VectorValue& VectorVolumeQuantity::getMultipleScatterValues() {
     return multiple_scatter_values_;
 }
 

src/MIDSX/Core/volume_quantity_container.cpp

@@ -17,7 +17,7 @@ void VolumeQuantityContainer::measureAll(TempVolumeTallyData& temp_tally_data) {
     }
 }
 
-VolumeQuantityContainer VolumeQuantityContainer::operator+(const VolumeQuantityContainer& other) const {
+VolumeQuantityContainer VolumeQuantityContainer::operator+(VolumeQuantityContainer& other) const {
     // check if either container is empty, if so return the other container
     if (vector_quantities_.empty()) {
         return other;
@@ -53,6 +53,11 @@ VolumeQuantityContainer VolumeQuantityContainer::operator+(const VolumeQuantityC
     return new_container;
 }
 
+void VolumeQuantityContainer::clear() {
+    vector_quantities_.clear();
+    count_quantities_.clear();
+}
+
 std::unordered_map<VectorVolumeQuantityType, VectorVolumeQuantity>& VolumeQuantityContainer::getVectorQuantities() {
     return vector_quantities_;
 }

src/MIDSX/Core/voxel_grid.cpp

@@ -96,7 +96,7 @@ std::unordered_map<int, VectorValue> VoxelGrid::getEnergyDepositedInMaterials()
 
     for (auto& voxel : voxels_) {
         int materialID = voxel.materialID;
-        std::vector<double> voxel_dose = voxel.dose.getVector();
+        std::vector<double>& voxel_dose = voxel.dose.getVector();
         energyDepositedInMaterials[materialID].addValues(voxel_dose);
     }