added a CS function through time to DMR; some code stubs for DMR tests

src/CompartmentalSystems/discrete_model_run.py

@@ -1324,8 +1324,14 @@ def backward_transit_time_quantiles_inputs_only(self, q):
         return res * self.dt
 
     def CS(self, k0, n):
+        """Carbon sequestration from ``k0`` to ``n``."""
         Phi = self._state_transition_operator_matrix
-        return sum([(Phi(n, k) @ self.net_Us[k]).sum() for k in range(k0, n+1, 1)])
+        return sum([(Phi(n, k+1) @ self.net_Us[k]).sum() for k in range(k0, n, 1)])
+
+    def CS_through_time(self):
+        Phi = self._state_transition_operator_matrix
+        return np.array([self.CS(0, n) for n in self.times])
+
 
 #    # return value in unit "time steps x dt[0]"
 #    def backward_transit_time_quantiles_from_masses(self, q, start_age_masses_at_age_bin):

tests/Test_discrete_model_run.py

@@ -662,3 +662,52 @@ def test_age_quantile_at_time(self):
     def test_backward_transit_time_quantile_from_density(self):
         raise
 
+    @unittest.skip("stub")
+    def test_age_densities_vs_mean_age_vector(self):
+        # this is some piece of code that can be reused to test the age densities
+        # versus the mean age system
+
+        def p2_sv(ai, ti):
+            Phi = dmr._state_transition_operator
+
+            if (ai < 0) or (ti <= ai):
+                return np.zeros((dmr.nr_pools,))
+
+            U = dmr.net_Us[ti-ai-1]
+            res = Phi(ti, ti-ai, U)  # age 0 just arrived
+
+            return res # dt=1
+
+        def p1_sv(ai, ti):
+            Phi = dmr._state_transition_operator
+
+            if (ai < 0) or (ai < ti):
+                return np.zeros((dmr.nr_pools,))
+
+            return Phi(ti, 0, p0(ai-ti))
+
+        p_sv = lambda ai, ti: p1_sv(ai, ti) + p2_sv(ai, ti)
+
+        mean_system_age_sv = lambda ti: sum([k * p_sv(k, ti).sum() for k in range(ti+1)]) / dmr.solve()[ti].sum()
+
+        print(mean_system_age)
+        print([mean_system_age_sv(ti) for ti in dmr.times])
+
+        dmr_p1_sv = dmr.age_densities_1_single_value_func(p0)
+        dmr_p2_sv = dmr.age_densities_2_single_value_func()
+        dmr_p_sv = dmr.age_densities_single_value_func(p0)
+
+        print(p_sv(3, 7) - dmr_p_sv(3, 7))
+
+    @unittest.skip("stub")
+    def test_CS(self):
+        self.assertEqual(dmr.Cs(0, 0), 0)
+        self.assertEqual(dmr.Cs(0, 1), dmr.net_Us[0].sum())
+        self.assertEqual(
+            dmr.Cs(0, 2), 
+            (self.Bs[0] @ dmr.net_Us[0] + dmr.net_Us[1]).sum()
+        )
+
+
+
+