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

Type stability #21

Closed
wants to merge 5 commits into from
Closed

Type stability #21

wants to merge 5 commits into from

Conversation

matbesancon
Copy link

A first step towards fixing #15 #19 #20

Still type unstable:

@code_warntype QuadGK.preprocess_real
@code_warntype QuadGK.do_quadgk

seemed ok.

To fix it, the API should be changed to take into account the different types

E
struct Segment{TB<:Number,TI,TE}
a::TB
b::TB
Copy link
Member

Choose a reason for hiding this comment

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

Note that this is a different constraint: currently we're saying that the types of a and b are both subtypes of Number, but here you're saying they have to be the same subtype of Number. Is this intentional?

Copy link
Author

Choose a reason for hiding this comment

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

yes indeed, I considered the boundaries of the domain would be of the same type, so this was intentional. Would there be some use cases where this is not the case?

Copy link
Member

Choose a reason for hiding this comment

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

No, not necessarily, I just wanted to make sure the change was intentional.

end
(s0,si) = inf1 ? (s2,s1) : (s1,s2)
if si < 0 # x = s0 - t/(1-t)
return do_quadgk_common(t -> begin den = 1 / (1 - t);
Copy link
Member

Choose a reason for hiding this comment

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

You can use a do block rather than t -> begin.

Copy link
Author

Choose a reason for hiding this comment

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

it was written like this, I found it a bit dense, I'll rewrite this part

@ararslan ararslan requested a review from stevengj June 3, 2018 17:14
src/QuadGK.jl Outdated
for i in 1:length(s) - 1
heappush!(segs, evalrule(f, s[i],s[i+1], x,w,gw, nrm), Reverse)
end
segs = sort([evalrule(f, s[i], s[i+1], x, w, gw, nrm) for i in 1:length(s)-1], order = Reverse)
Copy link
Member

Choose a reason for hiding this comment

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

shouldn't this be heapify! and not sort?

Copy link
Author

Choose a reason for hiding this comment

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

fixed

@stevengj
Copy link
Member

stevengj commented Jun 3, 2018

If the integrand is not type-stable, at some point you may need to reallocate the segs array (if an integrand evaluation yields a type that doesn't fit into the existing array).

@matbesancon
Copy link
Author

@stevengj a solution would be to make the integrand type stable then? Otherwise the instability will keep propagating through the call stack

src/QuadGK.jl Outdated
for i in 1:length(s) - 1
heappush!(segs, evalrule(f, s[i],s[i+1], x,w,gw, nrm), Reverse)
end
segs = heapify([evalrule(f, s[i], s[i+1], x, w, gw, nrm) for i in 1:length(s)-1], Reverse)
Copy link
Member

Choose a reason for hiding this comment

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

heapify!

mid = (seg.a + seg.b) * 0.5
seg1 = evalrule(f, seg.a, mid, x,w,gw, nrm)
seg2 = evalrule(f, mid, seg.b, x,w,gw, nrm)
heappush!(segs, seg1, Reverse)
Copy link
Member

Choose a reason for hiding this comment

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

Seems to need some kind of check here for seg1 isa eltype(segs) (similarly for seg2), and reallocate segs if not (will also require a refactored do_quadgk_common routine that can take an existing segs heap as a starting point).

Copy link
Author

Choose a reason for hiding this comment

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

this would mean all segments don't have the same type? Wouldn't that be strange, boundaries all have the same type?

@stevengj
Copy link
Member

stevengj commented Jun 6, 2018

a solution would be to make the integrand type stable then?

No way to do enforce this other than to throw an error, and it would be better to keep the routine generic to an arbitrary user-defined integrand. Type-stability should be an optimization, not a requirement.

Otherwise the instability will keep propagating through the call stack

No, basically the solution is similar to how map(f, x) is implemented for a non type-stable f. Each time f returns a new type, you widen the type of the segs array (reallocating) and re-call the inner integration routine recursively. In practice, non type-stable functions return only a few types (usually only 2), so the reallocation is not triggered often. For the common case of type-stable integrands it would still be efficient.

@matbesancon matbesancon mentioned this pull request Aug 21, 2018
@stevengj stevengj mentioned this pull request Jun 7, 2019
@stevengj stevengj closed this in #36 Jun 7, 2019
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