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

svector points #112

Open
wants to merge 20 commits into
base: main
Choose a base branch
from
Open

Conversation

skygering
Copy link
Collaborator

Switch functions that return geometries to have points of the form GI.Point{SVector} rather than a Tuple.

This addresses #106

I did not update barycentric.jl. I don't really understand what is going on there and it seems to depend on GeometryBasics functions. Do we also want that swapped over? I think it is fine because it returns the lambda values, not an actual object, correct?

@asinghvi17
Copy link
Member

Hmm, good point on barycentric - that code is a little messy but was basically trying to roundtrip the z-values. In this case I think tuples are probably fine (or at least no better or worse than static vectors).

const Edge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
const TupleEdge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T

const SVPoint{T} = GI.Point{false, false, SA.SVector{2, T}, Nothing} where T <: AbstractFloat
Copy link
Member

Choose a reason for hiding this comment

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

I'm not sure how to enable SVPoint{Z, M, T} - tried to do it the obvious way with:

const SVPoint{Z, M, T, CRS} = GI.Point{Z, M, SA.SVector{2+Z+M}, T, CRS}

but that complains about Z and M being type variables. Any thoughts on how we could accomplish this would be great :D

Copy link
Collaborator Author

@skygering skygering Apr 16, 2024

Choose a reason for hiding this comment

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

Yeah... this is also going to be needed for #97 so that we can actually store the extents... If we can't figure it out we could put that into the next PR.

I also don't know how to do it 😆 which is why I hard-coded for now

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Also, the functions that use it only currently work for 2D so it isn't the end of the world right now, but definitely a needed improvement.

src/methods/polygonize.jl Outdated Show resolved Hide resolved
function (::ClosedRing)(::GI.PolygonTrait, polygon)
exterior = _close_linear_ring(GI.getexterior(polygon))

function (::ClosedRing)(::Type{T}, ::GI.PolygonTrait, polygon) where T
Copy link
Member

Choose a reason for hiding this comment

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

especially here I'm not sure if this is necessary, since we're directly copying the point. Maybe enforcing standard point datatypes could be another fix, though!

Copy link
Collaborator Author

@skygering skygering Apr 16, 2024

Choose a reason for hiding this comment

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

The function was type unstable before since one branch of the LienarRing dispatch used tuples and the other didn't. So I just changed it so all branches use SVPoints.

Copy link
Member

Choose a reason for hiding this comment

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

Ah I meant the addition of T here. Changing to svpoints does make sense.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Then I would need to do eltype or get rid of the Float64 assumption or else it will just default anyways. Do we want to do that?

end
else
return apply(PointTrait(), geom; kw...) do p
(GI.y(p), GI.x(p))
GI.Point(SA.SVector{2}(GI.y(p), GI.x(p)))
Copy link
Member

Choose a reason for hiding this comment

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

should this be SVPoint?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

We could probably make something like that work with a basic constructor. I think I didn't do it since we didn't have a generalizable way so I wanted all of the branches to stay the same until we got that figured out...

Copy link
Member

Choose a reason for hiding this comment

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

Ah yeah a function svpoint{Z, M, [T]}(x, y, [z, [m]]) should work

if _is3d(geom)
if _ismeasured(geom)
return apply(PointTrait(), geom; kw...) do p
(T(GI.x(p)), T(GI.y(p)), T(GI.z(p)), T(GI.m(p)))
Copy link
Member

Choose a reason for hiding this comment

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

It isn't necessary that a geometry has both Z and M. This would have to be a little more complex (though I feel like we should define some form of framework for this, even if it's just conceptual and in the docs)

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Hmmmm.... It does say in the GI docs that for a point to be measured, it MUST be a 4-element tuple.

Copy link
Member

Choose a reason for hiding this comment

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

That's only for tuple and vector points:

Length 3 Tuple and Vector points

@@ -104,7 +104,7 @@ end
])
c1, area1 = GO.centroid_and_area(m1)
c1_from_LG = LG.centroid(m1)
@test c1[1] ≈ GI.x(c1_from_LG)
Copy link
Member

Choose a reason for hiding this comment

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

Ah is this because the points are now GI.Point? We should still allow getindex behaviour on those IMO @rafaqz, just for convenience's sake.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Yep. Would that be implemented in GI?

@rafaqz
Copy link
Member

rafaqz commented Apr 17, 2024

I dont 100% get the reason to put SVector in points... its only useful for dispatch and wrapping it loses that.

We could make points <: StaticVector ...

@asinghvi17
Copy link
Member

asinghvi17 commented Apr 17, 2024

I guess at this point we would have to make our own point type? In that case it may as well be backed by a static vector. (Since GI points can wrap any collection I'm not sure how fast we can get basic matrix-vector multiplication etc., given that we'd have to use GI.x and so on on them).

@rafaqz
Copy link
Member

rafaqz commented Apr 17, 2024

Well GI.Point is my own point type ;)

But the SVector dep might be too heavy for GI. We can dup it here as <: SVector if we really need that.

@skygering
Copy link
Collaborator Author

I dont 100% get the reason to put SVector in points... its only useful for dispatch and wrapping it loses that.

We could make points <: StaticVector ...

I thought that you wanted this @rafaqz ? If not, I honestly also don't like it. I think it is annoying to generate the polygons this way, since GI doesn't automatically do the wrapping... I would have rather just done SVectors but I thought this is what we talked about in the slack and #106.

Do we just want SVectors? Does having the point wrapper do anything that we want it to do?

@skygering
Copy link
Collaborator Author

From our conversation on Slack today, it seems like the best course of action is to create our own point type:

struct Point{N,T,Z,M} <: StaticVector{N,T}
    vals::NTuple{N,T}
end

This will allow dispatch off of Z and M, but still give us the benefits of StaticVectors, like indexing.

@skygering
Copy link
Collaborator Author

Okay! It is now working! Within types.jl I have added the new SVPoint type and have 2D, 3D, and 4D constructors to make points of the correct type. Should Z and M be booleans or BoolsAsTypes?

Should there be a more general constructor? We could dispatch off of the type of NTuple, but if anything but an NTuple is passed in, we will need to do a _ismeasured and a _is3d check.

This then leads into my next question. Should we try to generalize the _to_points and _to_edges code? Should we have different functions for TuplePoints and SVPoints? Because right now they just return 2D tuples. However, that really isn't apparent from the name and requires having the _sv_points function. I also tried to collect and consolidate the various functions that turn points and geometries into different point types within the utils.jl file.

Finally, I have it so that we are passing T around a lot so that the points can be a user-specified type. I think this is nice for type stability... but it doesn't work for the GEOS functions, for example, as LibGEOS only works with Float64...

Also, I know I need to add comments on the new types, but I wanted to make sure everything is correct first.

@asinghvi17
Copy link
Member

asinghvi17 commented Jun 13, 2024

Yes, I think a general constructor would be nice. Generally, it probably should be parametrized to avoid a bunch of if statements.

Z and M can be booleans as in the GI structs - the values are encoded in the type domain so it's not a problem for dispatch and specialization. We only really need the BoolsAsTypes when passing around values.

Passing T around should be fine since we do it a bunch...they are probably not applicable to GEOS but when we add GEOSConvert then they could be useful.

I have to read the PR more thoroughly (and test it out) but this is my first impression.

@skygering
Copy link
Collaborator Author

I will say - this doesn't even fix my problem with using transform. Here, the point type is SVPoint but after going through a transform the geometry will have point that are of type StaticVector. My big annoyance is that type change, which this doesn't even fix...

Thoughts on that?

@rafaqz
Copy link
Member

rafaqz commented Jun 13, 2024

I think I mentioned in the past that there is not so much use to having a wrapped SVector, as it hides all the StaticArrays.jl dispatch anyway. I think the options are:

  1. Tuple - basic
  2. Point - its a geometry at least
  3. SVector - it has nice math defined on it, and its AbstractArray

But not a Point{Svector}

@skygering
Copy link
Collaborator Author

It isn't a Point{SVector}. It is a subtype of SVector.

However, I am very, very happy to abandon this effort, as I like the tuples and only work in 2D. I just want transforms to also return tuples, which is an easy fix at the point level. So if that is an alternative PR for now, I would be glad to do that so at least transforms isn't changing the point type of geometries.

@asinghvi17
Copy link
Member

The fix for transform would just be changing StaticArray to SVPoint there, IMO.

return apply(PointTrait(), geom; kw...) do p
f(StaticArrays.SVector{3}((GI.x(p), GI.y(p), GI.z(p))))
SVPoint_4D(f(StaticArrays.SVector{4}((GI.x(p), GI.y(p), GI.z(p), GI.m(p)))))
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
SVPoint_4D(f(StaticArrays.SVector{4}((GI.x(p), GI.y(p), GI.z(p), GI.m(p)))))
SVPoint_4D(f(SVPoint_4D((GI.x(p), GI.y(p), GI.z(p), GI.m(p)))))

@rafaqz
Copy link
Member

rafaqz commented Jun 13, 2024

Ok lol I should read it first I misunderstood what SVPoint is.

Do we just need to add some similar methods? I've extended static vectors before and it is possible. (maybe not similar... but I remember some wrangling to get the same types back from base methods...)

Probably we want that anyway for other use cases.

@skygering
Copy link
Collaborator Author

Do you have any example code, @rafaqz ? I have never extended an SVector so I am a little lost.

I am also not sure how to make to_points and to_edges general without doing a conditional check for is_measured and is_3d. Thoughts on that? Do we just do the check?

@asinghvi17
Copy link
Member

May as well do the test - it should ideally compile away if that is known at compile time.

@rafaqz
Copy link
Member

rafaqz commented Jun 13, 2024

Is3d/ismeasured still won't compile away for some backends where it's a runtime check, so we want the checks as far out as possible... Like nested if blocks and four separate calls to apply for each combination will be the most stable

@rafaqz
Copy link
Member

rafaqz commented Jun 14, 2024

I'm not sure for SArray in this case, it could be missing constructors as well. Static arrays has so many fallbacks.

You could use Cthulhu and see where the type changes, and add that method.

Could be similar_type ?

https://github.com/rafaqz/LandscapeChange.jl/blob/main/src%2Fnamed_vector.jl#L81

This static NamedVector works really well, but a bunch of those methods are for NamedTuple behaviour so a bit confusing for you.



=#
struct SVPoint{N, T, Z, M} <: GeometryBasics.StaticArraysCore.StaticVector{N,T}
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Thoughts on this structure update @asinghvi17 and @rafaqz ?

Copy link
Member

Choose a reason for hiding this comment

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

The struct looks good to me!

Copy link
Member

@rafaqz rafaqz Jun 21, 2024

Choose a reason for hiding this comment

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

Yep looks good. N is kinda redundant but we need it for StaticArrays

@skygering
Copy link
Collaborator Author

Also all of my tests are passing locally so I am confused as to why they are failing on GitHub.

src/utils.jl Outdated
Comment on lines 75 to 84
function find_point_constructor(x, ::Type{T})
if _ismeasured(x)
SVPoint_4D
elseif _is3d
SVPoint_3D
else
SVPoint_2D
end

end
Copy link
Member

@asinghvi17 asinghvi17 Jun 20, 2024

Choose a reason for hiding this comment

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

Suggested change
function find_point_constructor(x, ::Type{T})
if _ismeasured(x)
SVPoint_4D
elseif _is3d
SVPoint_3D
else
SVPoint_2D
end
end
function find_point_constructor(x, ::Type{T})
M = _ismeasured(x)
Z = _is3d(x)
N = 2 + M + Z
return SVPoint{N, T, Z, M}
end

This is what we should have, IMO, but it does overlap in purpose with the SVPoint_4D etc constructors defined earlier. I wonder whether we should simply define a default and a fallback for three-dimensional points, since that's where all the uncertainty is...

Here, the constructor for SVPoint{N, T, Z, M} would be overloaded to provide fallbacks (zero(T)) in case the input is not measured or 3d. Usually, that does compile out in most cases where it can.

Copy link
Collaborator Author

@skygering skygering Jun 21, 2024

Choose a reason for hiding this comment

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

What do you mean? I don't understand what you mean. I like your constructor - it is very elegant, but seems type unstable. I don't really understand your fallback concept.

These (M = _ismeasured(x), Z = _is3d(x)) return false which is equivalent to zero(T) when it isn't measured or 3d. I don't understand how we would have a fallback without running those lines. And then if we are running those lines, we are basically within find_point_constructor that you defined.

The reason I defined the SVPoint_4D etc. constructors was to eliminate type instability within functions where we know what types we are working with (i.e. clipping).

@asinghvi17
Copy link
Member

All the errors are in the table integration tests. I'll take a look locally - it's possible that there's a missing sneaking in somewhere.

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.

None yet

3 participants