JuliaNLSolvers / ManifoldProjections.jl

MIT License
5 stars 3 forks source link

Change symbol of the norm and add a radius to the Sphere #3

Closed benneti closed 5 years ago

benneti commented 5 years ago

A few comments: As the old code would fail for project_tangent for x not on the surface corresponding to the sphere, I changed the code in a way that first projects x onto the surface and then does the projection. If this is the way you think is appropriate, line 24 should be changed to project_tangent(M::Manifold, g, x) = project_tangent!(M, copy(g), copy(x)) else line 79 should be deleted and 78 changed to project_tangent!(S::Sphere,g,x) = (xnorm = normalize(x); g .-= (real(dot(xnorm,g))).*xnorm)

If we want a Contructor that does some more checks we could use the following code

struct Sphere{T} <: Manifold where {T <: Real}
    r::T
    function Sphere(r)
        if r == nothing
            new{Nothing}(nothing)
        elseif typeof(r) <: Real
            r < 0 ? error("radius has to be a positive number!") : new{T}(r)
        else
            error("radius has to be a real number or nothing")
        end
    end
end

If you have any comments let me know
antoine-levitt commented 5 years ago

I don't think project_tangent should enforce x to be on the manifold. It's just assumed it is, otherwise it's hard to avoid doing extra work.

benneti commented 5 years ago

Should there be any error checking or not?

antoine-levitt commented 5 years ago

I would just do it like this:

struct Sphere{T<:Union{Real,Nothing}}
    r::T
    function Sphere(r::T) where {T <: Real}
        new{typeof(r)}(r)
    end
    function Sphere()
        new{Nothing}(nothing)
    end
end

That way you can't call Sphere(nothing) from the outside. I don't usually bother checking things like r < 0 (plus it makes me nervous if the code is type-stable or not), but as you wish.

antoine-levitt commented 5 years ago

Actually now I'm thinking about this let's not bother doing it this complicated way just to avoid dividing by 1. Just have a regular struct Sphere{T<:Real} with Sphere() = Sphere(1) and that's it.

benneti commented 5 years ago

Do you think checking whether r>0 should be done or is not necessary either?

antoine-levitt commented 5 years ago

Cf supra, it's up to you! I just checked, julia is clever or not to find out that error will never return and infers the return type.

benneti commented 5 years ago

if you are happy with the current implementation I'll squash the comits. I think raising an error is better, but if you think it is unnecessary I can remove the check, too. How would you like the dokumentation? I think it would be helul to have documentation for retract and project_tangent, that for example tells the user that x is assumed to lie on the surface of the manifold and not checked in any way.

antoine-levitt commented 5 years ago

That's great! I'll go ahead and merge it (no need to squash). The doc can be done separately (if you're motivated to do a PR it's very welcome! I'd just document the interface in the README.)

@pkofod btw, do you know how to setup CI?

pkofod commented 5 years ago

btw, do you know how to setup CI?

Yes, I'll do that.

benneti commented 5 years ago

Can one of you tell me wether Sphere(r::T) where {T <: Real} = r < 0 ? error("radius has to be a positive number!") : new{T}(r) or Sphere(r::Real) = r < 0 ? error("radius has to be a positive number!") : new{typeof(r)}(r) would have been better? Or whether there is any difference at all?

antoine-levitt commented 5 years ago

There's no difference. You can take a look at the generated code using the @code_ macros.

benneti commented 5 years ago

Cool, thank you