JuliaApproximation / DomainSets.jl

A Julia package for describing domains as continuous sets of elements
MIT License
72 stars 12 forks source link

WIP Add Segment #68

Closed dlfivefifty closed 6 months ago

daanhb commented 3 years ago

Shall we put this in 0.5.1?

daanhb commented 3 years ago

Looking at this, I was just reminded that ApproxFun defines a canonicaldomain function, as well as tocanonical and fromcanonical. These functions are also defined in DomainSets 0.5 (and exported). I'm pretty sure I've started using them because I saw the name in ApproxFun first a while ago, and it is a really good idea which enables lots of flexibility and automated procedures. But I forgot about where it came from later and did not check the expected interface, so there might be a clash now. Hopefully we can easily sort that out here, it doesn't look too bad actually.

I'll describe what the functions do in DomainSets:

In addition DomainSets defines a two-argument version of canonicaldomain where the second type indicates the kind of canonical domain we're after (there might be more than one). This allows a more generic way of checking whether two domains are equal, a way to find simpler domains and a generic description of parameterizations (typically these are maps from lower-dimensional domains).

Currently, there is:

There are corresponding functions fromcanonical(d, Parameterization()) etcetera. These would clash with fromcanonical(d, x), so I'll rethink them. The two-argument version fromcanonical(d, x) makes sense to avoid having to construct a map and I'd actually like to use it. It would probably be sufficient to put the canonical type first and then we can disambiguate, i.e., canonicaldomain(Parameterization(), d). This would lead to fromcanonical(d::Domain, x) and fromcanonical(::Parameterization, d::Domain, x). No clash. Or, the more general variant could have a different name. It is not really user-friendly to have to type canonicaldomain(d, DomainSets.Parameterization()) anyway, I defined parameterdomain(d) for that purpose.

In DomainSets itself, the generic implementation allows a very flexible mapto(d1,d2) function, which tries to find a map from d1 to d2:

julia> using DomainSets, StaticArrays

julia> m = mapto(2..3, Sphere(4.0, SA[0.5; 0.6]))
(x -> 4.0 * x + [0.5, 0.6]) ∘ F ∘ (x -> 1.0 * x + -2.0)

F = DomainSets.UnitCircleMap{Float64}()

The interval is first mapped to [0,1], which maps to the unit circle, which maps to the circle with different center and radius. This is tricky to do without ambiguities, but I think I resolved it. Composite maps define jacobians and all the components do as well, so we can do:

julia> jacobian(m, 2.7)
2-element SVector{2, Float64} with indices SOneTo(2):
 23.90265731793245
 -7.766444154901849

julia> (m(2.7+1e-6)- m(2.7))/1e-6
2-element SVector{2, Float64} with indices SOneTo(2):
 23.902681719922825
 -7.766369063855194

julia> jacobian(m)
x -> F₂(x) * F₄(x) * F₃(x)

F₁ = DomainSets.UnitCircleMap{Float64}()
F₂ = (x -> [4.0 0.0; 0.0 4.0]) ∘ F₁ ∘ (x -> 1.0 * x + -2.0)
F₃ = x -> 1.0
F₄ = LazyJacobian(F₁) ∘ (x -> 1.0 * x + -2.0)

julia> jacobian(m)(2.7)
2-element SVector{2, Float64} with indices SOneTo(2):
 23.90265731793245
 -7.766444154901849

Anyway, that's just to say why I like the canonical domain functionality.

Packages can define their own canonical type to make different choices from the default. I've done that in DomainIntegrals, to make intervals map to [-1,1] in fact :-) I'm also using parameterizations there to automate the evaluation of some integrals (hence why I needed jacobians), but I know of no other uses so far. I think it may still be fine to just change the interface at this stage.

A short term solution would be to stop exporting the canonical methods in DomainSets until the interface has settled down and is useful for ApproxFun?

dlfivefifty commented 3 years ago

To be honest I never really liked tocanonical and fromcanonical...

Perhaps canonicalmap? I think a Map can be thought of as an invertible function, whose inverse is specified by invmap (or whatever its called). Then tocanonical(d,x) becomes canonicalmap(d)(x) and fromcanonical(d,x) becomes invmap(canonicalmap(d))(x)... though perhaps this is even uglier...

daanhb commented 3 years ago

It could be mapto_canonical and mapfrom_canonical?

I found it works well if a single argument function returns a map and the two-argument version evaluates the map at a point. In some cases I've made the single-argument function return a lazy map by default (like LazyJacobian, LazyInverse), such that one only has to implement inverse(m,x) or jacobian(m,x). That's often cheaper to compute and it doesn't require having to define some new map type.

The inverse of a map right now is inverse(m), but a map might not be invertible. It could be invertible on its range (like a parameterization) or for a subset of its domain (like a projection), so I've also defined leftinverse and rightinverse. The map from a canonical domain to d is typically a parameterization, so there is the default: tocanonical(d) = leftinverse(fromcanonical(d)) and tocanonical(d, x) = leftinverse(fromcanonical(d), x).

A segment in 2D or 3D could be the map of a 1D interval, which would be a non-invertible rectangular map. That map does have a left inverse, which maps a point on the segment back to the interval. That would be my plan here.

I've added a lot of logic for vector-valued functions which probably ultimately belongs elsewhere, but most packages I've checked work with square maps only.

daanhb commented 3 years ago

I've changed the names to mapfrom_canonical and mapto_canonical in v0.5.2. It seems the remaining issues for getting ApproxFun to work with v0.5 are minor. The biggest one is the canonical interval, which currently is [0,1] in DomainSets and [-1,1] in ApproxFun. That gave me errors in the test suite of ApproxFunOrthogonalPolynomials. I've also removed the line that said ()(f::Function, g::Function) = x->f(x)g(x), because it strikes me as a dubious thing to do that may affect other packages. But apparently ApproxFun had a test for that (is it relevant?) and that test now fails too.

daanhb commented 6 months ago

Closing because this is outdated