# positivemax(m, sq; tol) drops sq and passes expand_connections' tol (1e-10)
# as si. Outer-connector stream mixing is then an unweighted average.
using InteractiveUtils
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using SciMLBase

@connector function SPort(; name)
    @variables p(t) m(t) [connect = Flow] h(t) [connect = Stream]
    System(Equation[], t, [p, m, h], []; name)
end

function Src(; name, rate, h0)
    @named port = SPort()
    @variables x(t) = h0
    System([port.m ~ -rate, port.h ~ x, D(x) ~ 0], t; name, systems = [port])
end

function Recv(; name)
    @named port = SPort()
    @variables z(t) = 0.0
    System([port.p ~ 1.0, port.h ~ z, D(z) ~ instream(port.h) - z], t;
           name, systems = [port])
end

# Two inner sources with unequal flows; one outer port. The generated
# outer-stream equation is gp.h ~ mix(s1.h, s2.h).
function Group(; name)
    @named gp = SPort()
    @named s1 = Src(rate = 1.0, h0 = 2.0)
    @named s2 = Src(rate = 2.0, h0 = 8.0)
    System([connect(gp, s1.port, s2.port)], t; name, systems = [gp, s1, s2])
end

@named g = Group()
@named c = Recv()
sys = mtkcompile(System([connect(g.gp, c.port)], t; name = :hier, systems = [g, c]))
prob = ODEProblem(sys, [], (0.0, 1.0))

got = SciMLBase.getsym(prob, sys.g.gp.h)(prob)
# Modelica §15.2 outer connection equation (si = 1+2 = 3 ≫ 1e-4):
# gp.h = (1*2 + 2*8) / (1+2) = 18/3 = 6
expected = 6.0
# Bug: _positivemax(m, 1e-10) ≈ 1e-4 for every positive m, so
# gp.h ≈ (2+8)/2 = 5 (unweighted).

println("bug2: positivemax drops sq; outer stream mix is unweighted")
println("ModelingToolkit ", pkgversion(ModelingToolkit))
println("observed gp.h equation:")
for eq in observed(sys)
    s = string(eq)
    if occursin("gp", s) && occursin("h", s) && occursin("positivemax", s)
        println("  ", s)
    end
end
println("got      g.gp.h = ", got)
println("expected g.gp.h = ", expected, "  (flow-weighted  (1*2 + 2*8)/3 )")
println("unweighted average would be ", (2.0 + 8.0) / 2)
versioninfo()