Potential issue in extending `tipping_probabilities` to attraction basins computed on non-uniform grid

[mdepitta]

@Datseris @awage , there is a potential issue with how tipping_probabilities are computed at the moment. The result is correct if the grid steps are unitary in all directions. If it is not, then the computation of the basin measure must consider the true grid size. Here is my temporary workaround:

function my_tipping_probabilities(basins_before, basins_after, grid)
    @assert size(basins_before) == size(basins_after)

    bid, aid = unique.((basins_before, basins_after))
    P = zeros(length(bid), length(aid))
    # Make -1 last entry in bid, aid, if it exists
    put_minus_1_at_end!(bid); put_minus_1_at_end!(aid)

    # Notice: the following loops could be optimized with smarter boolean operations,
    # however they are so fast that everything should be done within milliseconds even
    # on a potato
    for (i, ι) in enumerate(bid)
        B_i = findall(isequal(ι), basins_before)
        μ_B_i = 0
        for ind in B_i
            # μ = measure on non-unitary grid (in my specific case grid is from ranges/vectors of a state space of 4 variables)
            ## Maybe there is a way to directly sum from the vectors extracted from CartesianIndexes -- in python this would be immediate: but I do not know how to do it in Julia for now...
            μ_B_i += grid[1][ind[1]]*grid[2][ind[2]]*grid[2][ind[3]]*grid[4][ind[4]]
        end
#         μ_B_i = length(B_i) # μ = measure ## This was the original code instead
        for (j, ξ) in enumerate(aid)
           ## Directly compute the overlap
            overlap = findall(x->x.>0,(basins_pre.==ξ).&(basins_pst.==ξ))
            μ_overlap = 0
            for ind in overlap
                μ_overlap += grid[1][ind[1]]*grid[2][ind[2]]*grid[2][ind[3]]*grid[4][ind[4]]
            end           
            P[i, j] = μ_overlap/μ_B_i
        end
    end
    return P
end

function put_minus_1_at_end!(bid)
    if -1 ∈ bid
        sort!(bid)
        popfirst!(bid)
        push!(bid, -1)
    end
end

[awage]

Hi! Do you have a reference of what you are trying to achieve? If I understand correctly you want to extend the computation to non-uniform grid steps.

What does this multiplication of grid points implies: grid[1][ind[1]]*grid[2][ind[2]]*grid[2][ind[3]]*grid[4][ind[4]] ? I don't see the relation with the grid step.

[mdepitta]

@awage You are right. Indeed the grid should have been replaced by the equivalent step. Here is a temporary edit. I need to test it yet. Unfortunately, I am running out of memory as mentioned before, and I won't have access to HPC resources before the next three weeks. In my case, I have 4D system, and I am calling the grid steps by d<var_name>. However, if correct, the code below should be easy to generalize to ND systems. Unfortunately, not by my knowledge of the language at the moment.

function my_tipping_probabilities(basins_before, basins_after, grid)
    @assert size(basins_before) == size(basins_after)

    bid, aid = unique.((basins_before, basins_after))
    P = zeros(length(bid), length(aid))
    # Make -1 last entry in bid, aid, if it exists
    put_minus_1_at_end!(bid); put_minus_1_at_end!(aid)
     
    ## We now compute all the differences across all elements of the grid, so that we can compute all voxels accordingly
    ## Compute incremental steps of grid
    dg = diff(grid[1])
    df = diff(grid[2])
    da = diff(grid[3])
    dc = diff(grid[4])
    
    # Notice: the following loops could be optimized with smarter boolean operations,
    # however they are so fast that everything should be done within milliseconds even
    # on a potato
    for (i, ι) in enumerate(bid)
        B_i = findall(isequal(ι), basins_before)
        μ_B_i = 0
        for ind in B_i
            # μ = measure on non-unitary grid
            ## The condition is required because diff is length-1 the original vectors for the state variables upon which grid is defined 
            if (ind[1]<length(dg))&&(ind[2]<length(df))&&(ind[3]<length(da))&&(ind[4]<length(dc))
                μ_B_i += dg[ind[1]]*df[ind[2]]*da[ind[3]]*dc[ind[4]]
            end
        end
#         μ_B_i = length(B_i) # μ = measure
        for (j, ξ) in enumerate(aid)
            overlap = findall(x->x.>0,(basins_pre.==ξ).&(basins_pst.==ξ))
            μ_overlap = 0
            for ind in overlap
                if (ind[1]<length(dg))&&(ind[2]<length(df))&&(ind[3]<length(da))&&(ind[4]<length(dc))
                    μ_overlap += dg[ind[1]]*df[ind[2]]*da[ind[3]]*dc[ind[4]]
                end
            end           
            P[i, j] = μ_overlap/μ_B_i
        end
    end
    return P
end

[awage]

I see. I will think of a way to extend the computation to Ndim non uniform grids.

[awage]

I have a solution but it is untested and will not work at first. Can you please think of a test unit? I leave the code.

The idea is to restrict the size of the basins (with view) instead of cheking the bounds. Also I have written the computation of the volume using broadcast on the array. Maybe we can discuss about having this feature in the codebase with @Datseris.

Another solution is to have a convention for the element of the grid step missing. In some application the last step is repeated (in filtering for example).

function my_tipping_probabilities(basins_before, basins_after, grid)
    @assert size(basins_before) == size(basins_after)

    bid, aid = unique.((basins_before, basins_after))
    P = zeros(length(bid), length(aid))
    # Make -1 last entry in bid, aid, if it exists
    put_minus_1_at_end!(bid); put_minus_1_at_end!(aid)
     
    ## We now compute all the differences across all elements of the grid, so that we can compute all voxels accordingly
    ## Compute incremental steps of grid
    dg = diff.(grid)
    ng = ntuple(i -> 1:length(grid[i])-1, length(grid))
    bb = view(basins_before, ng...)
    ba = biew(basins_after, ng...)
    # Notice: the following loops could be optimized with smarter boolean operations,
    # however they are so fast that everything should be done within milliseconds even
    # on a potato
    for (i, ι) in enumerate(bid)
        B_i = findall(isequal(ι), bb)
        μ_B_i = 0
        for ind in B_i
            μ_B_i += prod(getindex.(dg, Tuple(ind))) 
        end
        for (j, ξ) in enumerate(aid)
            overlap = findall(x->x.>0,(bb.==ξ).&(ba.==ξ))
            μ_overlap = 0
            for ind in overlap
                    μ_overlap += prod(getindex.(dg, Tuple(ind))) 
            end           
            P[i, j] = μ_overlap/μ_B_i
        end
    end
    return P
end

[awage]

Corrected some typos and put a dumb example:

function my_tipping_probabilities(basins_before, basins_after, grid)
    @assert size(basins_before) == size(basins_after)

    bid, aid = unique.((basins_before, basins_after))
    P = zeros(length(bid), length(aid))
    # Make -1 last entry in bid, aid, if it exists
    put_minus_1_at_end!(bid); put_minus_1_at_end!(aid)
     
    ## We now compute all the differences across all elements of the grid, so that we can compute all voxels accordingly
    ## Compute incremental steps of grid
    dg = diff.(grid)
    ng = ntuple(i -> 1:length(grid[i])-1, length(grid))
    bb = view(basins_before, ng...)
    ba = view(basins_after, ng...)
    # Notice: the following loops could be optimized with smarter boolean operations,
    # however they are so fast that everything should be done within milliseconds even
    # on a potato
    for (i, ι) in enumerate(bid)
        B_i = findall(isequal(ι), bb)
        μ_B_i = 0
        for ind in B_i
            μ_B_i += prod(getindex.(dg, Tuple(ind))) 
        end
        for (j, ξ) in enumerate(aid)
            overlap = findall(x->x.>0,(bb.==ξ).&(ba.==ξ))
            μ_overlap = 0
            for ind in overlap
                    μ_overlap += prod(getindex.(dg, Tuple(ind))) 
            end           
            P[i, j] = μ_overlap/μ_B_i
        end
    end
    return P
end


function put_minus_1_at_end!(bid)
    if -1 ∈ bid
        sort!(bid)
        popfirst!(bid)
        push!(bid, -1)
    end
end

# Random basins 3 dim grid 
bas1 = rand(1:4,(10,10,10))
bas2 = rand(1:4,(10,10,10))
grid = (1:10, 1:10, 1:10) 

P = my_tipping_probabilities(bas1, bas2, grid)

[mdepitta]

@awage thanks. I need to test it on my data sets to be able to give you consistent feedback. I am currently building the data, but, as I mentioned earlier, I will have access to reasonable computing resources only in a couple of weeks. Thanks for your patience.

[mdepitta]

@awage Thank you for your patience. I got my laptop stolen and had to rewrite the code from scratch. I am testing your code on dummy data, and it works so far. Unfortunately, I am waiting to build the real attractors for my problem. But I am experiencing some issue that I am reporting in a separate thread.

[mdepitta]

@awage There is a typo in the computation of the overlapping basin. The overlapping indexes were originally computed by:

overlap = findall(x->x.>0,(bb.==ξ).&(ba.==ξ))

should instead be estimated by:

overlap = findall(x->x.>0,(bb.==ι).&(ba.==ξ))

based on the definition of tipping probabilities by Kaszás, Feudel & Tél. Tipping phenomena in typical dynamical systems subjected to parameter drift. Scientific Reports, 9(1).

[Datseris]

Yes, you are right, please do a PR that fixes it?

[mdepitta]

@Datseris I am not sure how to do a PR nor what a PR is. This code is not in the current Attractors.jl since it implements the tipping_probabilities that allows for basins defined on irregularly-spaced grids.

[Datseris]

okay no worries, I now understand you were referring to the code Alex pasted here, not to the code in the package.

[awage]

Thanks! Do you have an example that can be tested?

Also to submit a Pull Request (PR) you need to:

• fork the repository in your github environment.
• clone the repository locally on your computer.
• create a new branch in this repository.
• push your changes.
• Submit the PR in the github environment.

There are many youtube videos that will guide you through the process.