Rewritten rescale in lyapunovs.jl #214
Conversation
There was a problem hiding this comment.
That's a good start. However, now each time rescale! is called it allocates a new vector r. The best is to pre-allocate r = copy(u). Then, you can use broadcasting r .= u .+ (...)... so that you overwrite an existing r in-place instead of creating a new one each time.
src/chaosdetection/lyapunovs.jl
Outdated
| @@ -269,16 +269,18 @@ end | |||
| function lyapunov(pinteg, T, Ttr, Δt, d0, ut, lt) | |||
| # transient | |||
| t0 = pinteg.t | |||
| # array to preallocate state for rescale | |||
| r = [get_state(pinteg,1)] | |||
There was a problem hiding this comment.
I don't understand. Why is this a vector of vectors, instead of just copy(get_state(pinteg, 2))...?
There was a problem hiding this comment.
I was trying this copy(get_state(pinteg,2)) but then I get an r as an SVector. With that when I try to do
@. r=get_state(integ,1) + (get_state(integ,2)-get_state(integ,1))/a I get a set_index error as you cannot change the contents of an SVector.
src/chaosdetection/lyapunovs.jl
Outdated
| set_state!(integ,r,2) | ||
| function rescale!(integ,a,r) | ||
| # r =[get_state(pinteg,1)] for preallocated computation | ||
| r[1] = get_state(integ,1) + (get_state(integ,2)-get_state(integ,1))./a |
There was a problem hiding this comment.
This is still allocating. A lot. You need to broadcast the operations by adding . after the opereators, r .= u .+ v....
(and you will replace r[1] with r, as there is no reason to define r as a vector of vectors.)
|
Huh, why was this PR closed? It was going good with only a small amount of tweaking needed. |
|
Oh sorry, I was trying to make a PR on my forked repo and I think I messed up and closed this PR. I'll reopen it. |
SudoMishra
force-pushed
the
edit_rescale
branch
from
76a19d5
to
5e7da91
Compare
|
@SudoMishra test fail so something is going wrong with the way you rescale, can you have a closer look? |
src/chaosdetection/lyapunovs.jl
Outdated
| function rescale!(integ,a,r) | ||
| r .= get_state(integ,1) .+ (get_state(integ,2).-get_state(integ,1))./a | ||
| set_state!(integ,r,2) |
There was a problem hiding this comment.
You need two versions of this function that dispatch via multiple dispatch on r. If r::SVector then you cannot use .= and you just type r = (no dots in any formula). If instead r::AbstractArray then you use the above version.
There was a problem hiding this comment.
I've added the multiple dispatched rescale function. So test 5 is failing in lyapunov_exponents.jl and its a bit strange why that is happening. It is the case where r::Abstract_Array
When I do rescale with no preallocation i.e.
r = get_state(integ,1) .+ (get_state(integ,2).-get_state(integ,1))./a
the test is passing and I'm getting the correct lyapunov exponent.
However, when I simply change the above line to use the preallocated r, i.e.
r .= get_state(integ,1) .+ (get_state(integ,2).-get_state(integ,1))./a
The test fails as the lyapunov value that I get is wrong (around 18.3..).
I'm trying to figure out why this is happening.
There was a problem hiding this comment.
Hi @Datseris , I was trying to resolve this issue but its getting a bit strange. So, in the lyapunov function I've created a new variable r=copy(get_state(pinteg,2)). As I've used copy() for creating r, it should now be independent of pinteg. So it should only change when the rescale! function is called, where it used as a placeholder for preallocation. However, that is not the case, and it even changes outside the rescale!. Due to this, during the set_state! in rescale! a wrong value is being passed on which results in the test failing. I'm not sure why this is happening, I've also used deepcopy instead of copy but the problem still persists. Am I missing some basic thing here? I would be grateful if you could give some advice on how to resolve this.
There was a problem hiding this comment.
Can you repeat the same with deepcopy on r?
There was a problem hiding this comment.
Yeah, it is happening even with deepcopy. I was visualising the variables using the debugger and even with deepcopy the r variable gets updated during step! . And apparently this is only affecting one of the test cases. It should only get updated during rescale because that is the only place where I use it.
There was a problem hiding this comment.
Actually, I"m still struggling to understand what exactly is your problem with respect to the algorithm. Why does it matter if r changes? r is a dummy variable. Its correct value is set exactly before giving this to set_state. Even if r and the 2nd state are referencing each other, we still don't care. ...?
There was a problem hiding this comment.
I'm actually a bit confused myself. Even I think that since the correct value of r is getting changed just before being sent to set_state it shouldn't matter and the algorithm should work. But for some reason, test 5 in lyapunov_exponents was failing. I tried to figure out why this was happening and the only wierd thing that I saw happening was that r was changing outside rescale. I still think my code should be working just by set_state!(integ,r,2) but its failing for that specific case. I'll try to add the dispatch method to DynamicalSystemsBase.jl .
There was a problem hiding this comment.
Okay. I'll look over this locally over the weekend and will identify the error. Thanks for investigating as well.
There was a problem hiding this comment.
I have fixed the set_state! and opened a PR in DynamicalSystemsBase.jl. Using the corrected set_state, set_state!(integ,r,2) is working properly and its passing all tests. I can commit the proper rescale function once the DynamicalSystemsBase.jl PR is merged and ChaosTools uses that.
There was a problem hiding this comment.
I have added the updated rescale function without copy(r) in set_state. One test case will not pass currently as ChaosTools.jl currently uses the wrong set_state function from DynamicalSystemsBase.jl . I have added a PR to fix that in DynamicalSystemsBase.jl. For u::AbstractVector, set_state should use .=.
SudoMishra
force-pushed
the
edit_rescale
branch
from
bd2822f
to
afa85ba
Compare
I rewrote the rescale function in lyapunov.jl with get_state and set_state functions as you had mentioned.
In the current form of the rescale function, we no longer need multiple rescale functions as that can be handled by the multiple dispatches of get_state and set_state for different scenarios.
I'm unsure about the type of declarations that would go with the new rescale function but I think as now it handles both vector and matrix scenarios, the type declaration might not be needed.
I'm new to Julia and git. This is my first PR. I'm quite interested in contributing more to this package.
Thank you.