Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Detect cases where there's fewer brownians than equations, but noise still 'diagonal' #2914

Merged
merged 9 commits into from
Aug 2, 2024

Conversation

MasonProtter
Copy link
Contributor

@MasonProtter MasonProtter commented Jul 31, 2024

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

Additional context

Followup to #2882 and #2886

Comment on lines 754 to 777
@testset "Diagonal noise, less brownians than equations" begin
@parameters σ ρ β
@variables x(tt) y(tt) z(tt)
@brownian a b
eqs = [D(x) ~ σ * (y - x) + 0.1a * x, # One brownian
D(y) ~ x * (ρ - z) - y + 0.1b * y, # Another brownian
D(z) ~ x * y - β * z] # no brownians -- still diagonal
@mtkbuild de = System(eqs, tt)

u0map = [
x => 1.0,
y => 0.0,
z => 0.0
]

parammap = [
σ => 10.0,
β => 26.0,
ρ => 2.33
]

prob = SDEProblem(de, u0map, (0.0, 100.0), parammap)
@test solve(prob, SOSRI()).retcode == ReturnCode.Success
end
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Tested this example and it failed before this PR.

@MasonProtter
Copy link
Contributor Author

The downstream failure in Catalyst is caused by this, not entirely sure why it failed though it 's a bit weird. I'll try a more explicit version of the check for whether the noise is diagonal.

@MasonProtter
Copy link
Contributor Author

MasonProtter commented Aug 1, 2024

Okay, I think I figured out the problem, I was iterating over the axes of the noise equations in the wrong order, effectively looking at the transpose of the noise equations.

This bug was masked in the test suite, but it reared its head in the Catalyst tests because it creates some reaction systems where the noise matrix is of the form

[a 0 0 0 e 0 0 0
 0 b 0 0 0 f 0 0
 0 0 c 0 0 0 g 0
 0 0 0 d 0 0 0 h]

and since I was iterating over the transpose of that, I was mistakenly detecting it as diagonal noise.

Should be fixed now, and I've updated the non-diagonal noise test to watch for this case and ones like it

@MasonProtter
Copy link
Contributor Author

MasonProtter commented Aug 1, 2024

Hm, so now the failing Catalyst test is happening because Catalyst does this:

    p_matrix = zeros(length(get_unknowns(sde_sys)), numreactions(rs))
    return SDEProblem(sde_sys, u0map, tspan, pmap, args...; check_length,
        noise_rate_prototype = p_matrix, kwargs...)

in https://github.com/SciML/Catalyst.jl/blob/v14.1.0/src/reactionsystem_conversions.jl#L767. What this ends up doing is hard-coding into the SDEProblem that the noise is going to be a matrix, rather than diagonal, but that interacts badly with the change added in #2886, which will unconditionally make g assume the noise matrix is a vector, regardless of what the user is doing with the noise equations.

I've removed that change from #2886 and things seems to work now, we'll see how this test-run goes. I don't think that change was necessary if I understand correctly.

@isaacsas
Copy link
Member

isaacsas commented Aug 1, 2024

If calling SDEProblem on an SDESystem now infers the prototype optimally we could drop that in Catalyst.

@isaacsas
Copy link
Member

isaacsas commented Aug 1, 2024

By drop I just mean we won’t pass a prototype.

@MasonProtter
Copy link
Contributor Author

MasonProtter commented Aug 1, 2024

Yeah, dropping that would be be good because some cases in Catalyst would then get a performance boost, but it's also important that we not break code that does assume a specific form of the prototype like that, so I think the change I made in fe1a3ee is still necessary (though perhaps I misunderstood something, and @AayushSabharwal can chime in as he was the one who added it).

@isaacsas
Copy link
Member

isaacsas commented Aug 1, 2024

I’ll try to update that next week when I’m back from vacation.

@MasonProtter
Copy link
Contributor Author

@ChrisRackauckas can you merge this?

@ChrisRackauckas ChrisRackauckas merged commit e64c479 into SciML:master Aug 2, 2024
23 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants