rusandris / StateTransitionNetworks.jl

Toolkit for dynamics on state transition networks
1 stars 0 forks source link

Incorrect results for `q=0` with `renyi_entropy` #23

Closed rusandris closed 3 months ago

rusandris commented 3 months ago

Incorrect result for the q = 0 case. MWE:

using DynamicalSystems
using StateTransitionNetworks

ds = Systems.henon()

orbit,t = trajectory(ds,10^6;Ttr=10^3)
grid_size = 5
sts = timeseries_to_grid(orbit, grid_size)
P,Q,x = calculate_transition_matrix(sts)

renyi_entropy(P,0.0)
renyi_entropy(P,0.001)
rusandris commented 3 months ago

Broadcasting also does the requested operation (P_q = P .^ q) on the zero elements of the sparse matrix. A solution is to make a copy of the matrix, and do the operation only on the nonzeros(P) (zeroth power makes the matrix dense with ones):

P_q = deepcopy(P)
nonzeros(P_q) .= nonzeros(P_q) .^ q

Now the results are also correct for q=0:

julia> renyi_entropy(P,0.0)
0.7987545412073304

julia> renyi_entropy(P,0.001)
0.7983849420037713

See the fix here 27989a1

rusandris commented 3 months ago

Also, eigsolve can return/ target certain eigenvalues. With setting the which argument to :LM, eigenvalues of largest magnitude are returned first (which is the default btw). This way, the line in renyi_entropy

λ_max, = findmax(abs.(λs))

can be omitted.

rusandris commented 3 months ago

In

function renyi_entropy(P::SparseMatrixCSC{Float64, Int64}, q::Float64; x::Vector{Float64}=stationary_distribution(P)...

x is not needed at all => should be deleted

rusandris commented 3 months ago

Closing issue after d1c7021 .