[Bug] Wrong behavior for the life extraction

docs/src/index.md

@@ -110,8 +110,9 @@ $$\mathbf{E}(P) = \int_0^\inf S_p (t) dt = \sum_{j=0}^\inf \frac{S_p(t_j)}{\lamb
 
 Two approximations are made when the life gets out of the life table:
 
-- The last line of the ratetable is assumed to last until eternity. Indeed, the last line represents persons that are already 110 years old, and thus assuming that their future death rates are constants is not that much of an issue. 
-- When on the other hand a life exits the ratetable from the right, i.e. into the future but at a young age, we assume the last column of the rate table to define the future for this person.
+- When the person reaches the last line, she continues to move right. 
+- When the person reaches the last column, she continues to moves down.
+- The last cell lasts for eternity. 
 
 All this is implemented as a method for the `Distributions.expectation` function, since Lifes are random variables:
 

src/Life.jl

@@ -3,40 +3,43 @@
 
 This function returns a random variable that correspond to an extracted Life from the `BasicRateTable` at age `a` and date `d`. 
 
-This works by checking if the individual is closer to the oldest age than the last year in the ratetable, calculating at each step the time difference and the hazard values. For the younger individuals, we assume they go through the last column at the end no matter what age they are. 
+This works by checking if the individual is closer to the oldest age than the last year in the ratetable, calculating at each step the time difference and the hazard values. When reaching the borders of the table we follow them until the last cell, which runs till infinity. 
 """
 struct Life<:Distributions.ContinuousUnivariateDistribution
     ∂t::Vector{Float64}
     λ::Vector{Float64}
     function Life(brt::BasicRateTable,a,d)
+        # Current position in the life table: 
         i, j = dty(a, brt.extrema_age...), dty(d, brt.extrema_year...)
-
+
+        # Next position: Happy new year (right, →) or happy birthday (below, ↓) first ?
+        # Cap obtained values to avoid going out of bounds.
         rem_a = RT_DAYS_IN_YEAR - rem(a, RT_DAYS_IN_YEAR)
         rem_d = RT_DAYS_IN_YEAR - rem(d, RT_DAYS_IN_YEAR)
-
-        # Do we go right or below first ? 
-        # Happy birthday (below) or Happy new year (right) first ?
-        k,l = rem_a < rem_d ? (i+1,j) : (i,j+1)
-
-        # Cap obtained values to avoid going out of bounds: 
-        K,L = size(brt.values)
-        k = min(k, K)
-        l = min(l, L)
-
-        # lengths and hazards in the first two cells:  
+        k, l = rem_a < rem_d ? (i+1,j) : (i,j+1)
+        K, L = size(brt.values)
+        k, l = min(k, K), min(l, L)
+
+        # We can extract time spend and daily hazard in the first two cells,
+        # Either going →↓ or ↓→ as needed.
         ∂t = [min(rem_a,rem_d), abs(rem_a - rem_d)]
         λ  = [brt.values[i,j],  brt.values[k,l]]
 
-        while (k < K) && (l < L)
+        # Then go to the next two cells: 
+        while (k < K) || (l < L)
             i,j,k,l = i+1, j+1, k+1, l+1
-            push!(∂t, RT_DAYS_IN_YEAR - ∂t[2], ∂t[2])
-            push!(λ, brt.values[i,j], brt.values[k,l])
-        end
-        if (l >= L) # exit on the right => still young ! 
-            # A good approximation is to go through the last column. 
-            for m in (k+1):K
-                push!(∂t, RT_DAYS_IN_YEAR)
-                push!(λ, brt.values[m,end])
+            if (k < K) && (l < L) # We go →↓ or ↓→
+                push!(∂t, RT_DAYS_IN_YEAR - ∂t[2], ∂t[2])
+                push!(λ, brt.values[i,j], brt.values[k,l])
+            else
+                if (l < L) # We exit →→→→→→→→ along the bottom row
+                    append!(∂t, ones(L-l) .* RT_DAYS_IN_YEAR)
+                    append!(λ, brt.values[end, (l+1):end])
+                else # We exit ↓↓↓↓↓↓↓↓ along the last column
+                    append!(∂t, ones(K-k) .* RT_DAYS_IN_YEAR)
+                    append!(λ, brt.values[(k+1):end, end])
+                end
+                break
             end
         end
         return new(∂t,λ)