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

WIP real-valued functional #1333

Open
wants to merge 17 commits into
base: master
Choose a base branch
from
Open

Conversation

miallo
Copy link
Contributor

@miallo miallo commented Apr 27, 2018

I started working on the implementation of the range of Functional as mentioned in #1328.

If I am not mistaken the old implementation of FunctionalRightVectorMult leads to the functional “forgetting” its grad_lipschitz and linear properties.

I tagged some places in functional.py with “# HELP” where I am especially unsure about what I am doing.

I did NOT implement any kind of typecasting anywhere yet but just tried to make sure that the range property gets passed along.

@pep8speaks
Copy link

pep8speaks commented Apr 27, 2018

Checking updated PR...

No PEP8 issues.

Comment last updated on June 14, 2018 at 22:12 Hours UTC

@@ -686,7 +693,12 @@ def __init__(self, func, vector):
''.format(func))

OperatorRightVectorMult.__init__(self, operator=func, vector=vector)
Functional.__init__(self, space=func.domain)
# HELP: before it was Functional.__init__(self, domain=func.domain).
# But doesn't this loose the information on grad_lipschitz and linear?
Copy link
Member

Choose a reason for hiding this comment

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

Yes, it seems like the previous implementation forgets the linearity and lipschitz constant for the gradient. However, wouldn't the latter change for this functional? Anyway, we are thinking of removing grad_lipschitz since it is not currently used anywhere. See #1300

Copy link
Contributor Author

Choose a reason for hiding this comment

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

You are right and if I am not mistaken, it would change in a nontrivial way. My guess for the upper bound would be something like func.grad_lipschitz * ||constant||_\infty? Should I implement that or will it be scrapped anyways?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Wait... This only works for linear operators, correct? So just ignore it and make it NaN?

@@ -1020,11 +1033,11 @@ def __init__(self, func, quadratic_coeff=0, linear_term=None,
raise ValueError(
"Complex-valued `constant` coefficient is not supported.")
self.__constant = constant.real

#HELP: Need to check for range? The + should take care of it, right?
Copy link
Member

Choose a reason for hiding this comment

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

Since FunctionalQuadraticPerturb inherits from Functional the super call is going to execute the code in Functional.__init__. So if the test there is sufficient (which I think it is), then this should be fine.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It has been since forever that I worked with subclasses... 😅 (Which is probably a good indicator for my terrible style of programming)

@@ -1151,10 +1164,10 @@ def __init__(self, left, right):
if not isinstance(right, Functional):
raise TypeError('`right` {} is not a `Functional` instance'
''.format(right))

#HELP: Check ranges?
Copy link
Member

Choose a reason for hiding this comment

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

Should not be needed; it is done in OperatorPointwiseProduct.__init__.

@@ -1216,9 +1229,11 @@ def __init__(self, dividend, divisor):

self.__dividend = dividend
self.__divisor = divisor


#HELP: Check ranges?
Copy link
Member

Choose a reason for hiding this comment

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

Yes, here you would have to check the range.

Copy link
Member

Choose a reason for hiding this comment

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

Actually, let me be a bit more specific. What needs to be checked is that the range of the two functionals are the same. Since this only inherits from Functional, there is no check elsewhere that will make sure that this is the case. In the same way, there is already a check that the domains are the same.

@@ -1303,9 +1318,9 @@ def __init__(self, func):
if not isinstance(func, Functional):
raise TypeError('`func` {} is not a `Functional` instance'
''.format(func))

#HELP: The range should be equal?!
Copy link
Member

Choose a reason for hiding this comment

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

Actually, I think that here is where it becomes tricky.... if the range is complex then how should one define the convex conjugate? This is related to #590. So my current best answer is: I am not sure of what the range of this functional should be to make most sense.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Now we could just check for functional.range == RealNumbers and otherwise throw an error

Copy link
Member

@kohr-h kohr-h Apr 27, 2018

Choose a reason for hiding this comment

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

Yes we have had this issue from the beginning and just decided not to do anything about it. So IMO let's not raise any exceptions yet, instead do the right thing for RealNumbers and otherwise do whatever we were doing all along. At some point this whole complex functionals thing needs a thorough analysis and solution (or maybe not?), but this is not the time.

if range is None:
range = domain.field
elif range != domain.field:
linear = False
Copy link
Member

Choose a reason for hiding this comment

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

After some very recent discussions (see #1331), we decided for practical reasons to go with linearity in the sense "C = R^2". That would mean here that this elif case should go away.

@miallo miallo changed the title WIP realvalued functional WIP real-valued functional Apr 27, 2018
super(IndicatorLpUnitBall, self).__init__(domain=domain, linear=False,
range=range)
# HELP: The LpNorm should always have a range in the real numbers so
# that it is comparable to 1, right?
Copy link
Member

Choose a reason for hiding this comment

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

Yeah, so you can override the Functional default before the super call:

if range is None:
    range = RealNumbers()

And then see what happens 😆

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Isn't this more like what we want to have: self.__norm = LpNorm(domain, exponent, range=RealNumbers())
Then the range of the class itself still can be complex if needed

Copy link
Member

Choose a reason for hiding this comment

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

Oh I thought you were in LpNorm, didn't realize it was the indicator. There, yes, I agree with your suggestion.

@kohr-h
Copy link
Member

kohr-h commented Apr 27, 2018

Regarding checks, general remark:

  • If __init__ inputs are simply passed on to other __init__ methods, the checks done there are sufficient, unless a subclass needs its inputs to satisfy additional conditions.
    def __init__(self, domain, range):
        super(...).__init__(domain, range)   # OK, parent class does the checking
  • If before passing them to another __init__, inputs are compared to other things with is or == it's also okay not to check types and let the __eq__ implementation handle it.
    def __init__(self, domain, range=None):
        if domain == RealNumbers():  # OK
            # do stuff
        if range is None:  # OK
            # do other stuff
        super(...).__init__(domain, range)
  • If before passing them to another __init__, attributes of the inputs are used, we usually check beforehand to avoid errors like 'Foo' has no attribute 'bar':
    def __init__(self, domain, range=None):
        if range is None:
            # not great, if domain has no `field` attribute, error message is bad
            range = domain.field
        if not isinstance(domain, LinearSpace):
            # also not great, too restrictive (excludes `RealNumbers` for instance)
            raise ...
        # Better: duck-typing
        default_range = getattr(domain, 'field', RealNumbers())  # version 1, default RealNumbers
        default_range = getattr(domain, 'field', None)  # version 2, no default
        if range is None:
            if default_range is None:  # for variant 2 where there's no default
                raise ...
            range = default_range
    
        super(...).__init__(domain, range)

Domain of the functional.
exponent : float
Exponent for the norm (``p``).
"""
super(LpNorm, self).__init__(
space=space, linear=False, grad_lipschitz=np.nan)
domain=domain, linear=False, grad_lipschitz=np.nan, range=range)
Copy link
Member

Choose a reason for hiding this comment

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

So this is step 1: allow to manually set the range. I'm wondering about step 2: the sensible default for norms and such. What about

if range is None and getattr(domain, 'is_complex', False):
    range = RealNumbers()
super(...)

@@ -500,6 +503,7 @@ def convex_conj(self):
[BC2011] Bauschke, H H, and Combettes, P L. *Convex analysis and
monotone operator theory in Hilbert spaces*. Springer, 2011.
"""
# HELP: I guess it needs `range`, too
Copy link
Member

Choose a reason for hiding this comment

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

Yes, you'll want to use self.range (it makes sense IMO, and since it's a @property you can't pass arguments).


if (operator is not None and vector is not None and
vector not in operator.domain):
raise ValueError('domain of `operator` and space of `vector` need '
'to match')

# HELP: operator.range == operator.domain to be able to evaluate .inner
# => functional.range = operator.domain.field ?
Copy link
Member

Choose a reason for hiding this comment

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

True, this check is missing (or delayed until evaluation)

@@ -2154,9 +2168,10 @@ def __init__(self, space, outer_exp=1, singular_vector_exp=2):
>>> norm(space.one())
inf
"""
# HELP: range should be RealNumbers since comparable?
Copy link
Member

Choose a reason for hiding this comment

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

It's a norm, so yes

@@ -2262,8 +2277,10 @@ def __init__(self, functional, sigma=1.0):
>>> l1_norm = odl.solvers.L1Norm(space)
>>> smoothed_l1 = MoreauEnvelope(l1_norm)
"""
# HELP: Does this only makes sense, if
# functional.range == RealNumbers()? Otherwise raise Error?
Copy link
Member

Choose a reason for hiding this comment

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

No, there's a definition adapted for complex numbers, where the minimum is taken over the real part of the inner product only. Leave this as it is for now.

Copy link
Member

Choose a reason for hiding this comment

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

Leave this as it is for now.

That is, do add the range parameter.

@@ -326,7 +326,7 @@ def __init__(self, field):
"".format(name))

linear = name in LINEAR_UFUNCS
Functional.__init__(self, space=field, linear=linear)
Functional.__init__(self, domain=field, linear=linear)
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Why is there an explicit keyword passed in here for space/domain?
I have absolutely no experience with unit testing, so someone else needs to write the tests...

Copy link
Member

Choose a reason for hiding this comment

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

It's just for clarity.

If you write Functional.__init__(domain, range, linear) it's obvious what each of the arguments does, and keywords just add noise: Functional.__init__(domain=domain, range=range, linear=linear).
On the other hand, something like Functional.__init__(field, True) doesn't convey any information about the purpose of the arguments, so here Functional.__init__(domain=field, linear=True) is much clearer.

But why is this related to tests?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I stumbled across it in the test results and wanted to make sure that someone else takes over when I am done with implementing most of the stuff

Copy link
Member

@adler-j adler-j left a comment

Choose a reason for hiding this comment

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

I did a quick review. Overall nice changes, but some stylistic stuff etc needs to be adressed.

The main issue with this is that it's a sweeping change but without any tests added. I'd prefer proper testing to be performed.

@@ -223,6 +230,9 @@ def __init__(self, vfspace, exponent=None):
0 and 1 are currently not supported due to numerical
instability. Infinity gives the supremum norm.
Default: ``vfspace.exponent``, usually 2.
range : `Field`, optional
Range of the functional. The default `None` will be treated as
`vfspace.field`.
Copy link
Member

Choose a reason for hiding this comment

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

We need to decide on a default value here. Either we go for RealNumbers for all functionals that can only return real values anyway, or we go for domain.field. I'm personally not quite sure which convention is the best, but I guess the msot restrictive would be best, which indicates that we should go for RealNumbers whenever the functional is purely real valued.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I would agree, that probably the default for something that is real-valued should be real from a users perspective. But you should take into account that this leads to changes in the behaviour of existing code. (But as in stated in #1055 it might be reasonable to do it for the future)

Oh, and by the way: I don't want to write any tests because 1) I have never done it before and 2) Someone else might find flaws in my reasoning that I wouldn't


@property
def exponent(self):
"""Exponent corresponding to the norm."""
return self.__exponent

@property
def range_(self):
Copy link
Member

Choose a reason for hiding this comment

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

This is not needed. IndicatorLpUnitBall inherits from Operator and hence already has a range property, which is the one that should be used.

Domain of the functional.
range : `Field`, optional
Range of the functional. The default `None` will be treated as
Copy link
Member

Choose a reason for hiding this comment

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

The standard notation is

Range of the functional. Default: `RealNumbers`

@@ -1547,7 +1599,7 @@ def proximal(self):
See Also
--------
odl.solvers.nonsmooth.proximal_operators.\
proximal_convex_conj_kl_cross_entropy :
proximal_convex_conj_kl_cross_entropy :
Copy link
Member

Choose a reason for hiding this comment

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

This will break the link. It has to stay un-indented

@@ -1655,7 +1707,8 @@ def __init__(self, *functionals):
domain = ProductSpace(*domains)
linear = all(func.is_linear for func in functionals)

super(SeparableSum, self).__init__(space=domain, linear=linear)
super(SeparableSum, self).__init__(domain=domain, linear=linear,
range=functionals[0].range)
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 about this choice. One would perhaps expect that real+complex should be complex, but with this implementation it depends on the order.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I understand your reasoning and personally I would agree with you. I just thought I went with the convention of operators, where the domain (including dtype) has to be the same. But you are certainly right, that this should not depend on the order of the ops. I will think about it

"""

def __init__(self, functional, sigma=1.0):
def __init__(self, functional, sigma=1.0, range=RealNumbers()):
Copy link
Member

Choose a reason for hiding this comment

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

Should be None

Copy link
Member

Choose a reason for hiding this comment

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

When did we make that decision? I vaguely remember a discussion.

@@ -2364,6 +2441,7 @@ def __init__(self, space, gamma):
>>> abs(huber_norm(x) - l1_norm(x)) < tol
True
"""

Copy link
Member

Choose a reason for hiding this comment

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

We typically do not add a newline after the docstring like this

odl/set/sets.py Outdated
return float(inp)
else:
return 0.0
return float(getattr(inp, 'real', 0.0))
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Is this the way to go?

Copy link
Member

Choose a reason for hiding this comment

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

I wouldn't do it this way. The issue is that this maps anything which hasn't got a .real attribute to 0.0. Core types like int or bool will play nicely with this, but anything else that could be mapped to a float wouldn't. For instance strings or custom classes that implement __float__, but not real.

So I'd prefer this:

if inp is None:
    return 0.0
else:
    return float(getattr(inp, 'real', inp))

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I hadn't thought about that... That is a nice solution!

for functional in functionals):
range = RealNumbers()
else:
range = Integers()
Copy link
Contributor Author

@miallo miallo Apr 28, 2018

Choose a reason for hiding this comment

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

Should we even consider Integer functionals? Is there any use for them?

@@ -223,6 +229,8 @@ def __init__(self, vfspace, exponent=None):
0 and 1 are currently not supported due to numerical
instability. Infinity gives the supremum norm.
Default: ``vfspace.exponent``, usually 2.
range : `Field`, optional
Range of the functional. Default: `vfspace.field`.
Copy link
Member

Choose a reason for hiding this comment

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

For this kind of thing you need to use double backticks.
With single backticks, Sphinx will try to generate a hyperlink, which fails because Sphinx doesn't know the type of vfspace. Single backticks are usable for class names and glossary terms mostly, other things that are not referable need to go in double backticks.

return '{}({!r}, exponent={})'.format(self.__class__.__name__,
self.domain,
return '{}({!r} -> {!r}, exponent={})'.format(self.__class__.__name__,
self.domain, self.range,
Copy link
Member

@kohr-h kohr-h Apr 29, 2018

Choose a reason for hiding this comment

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

We try, whenever possible, to make repr return a string that is a valid Python expression. So please add range as a second (comma-separated) argument.
Note that there's a quite huge overhaul of all this repr stuff going on, see #1276 for reference. You don't need to spend too much effort on this part.

Copy link
Member

Choose a reason for hiding this comment

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

Actually, after looking at the failing doctest (output L2Norm(rn(3) -> RealNumbers()) doesn't match expected output L2Norm(rn(3)), something like this), I'd suggest not changing anything in __repr__. That's suboptimal for the non-standard case but at least doesn't break the printout for the default case.

@@ -1106,8 +1105,14 @@ def __init__(self, left, right, tmp_ran=None, tmp_dom=None):
rn(3).element([ 2., 4., 6.])
"""
if left.range != right.range:
raise OpTypeError('operator ranges {!r} and {!r} do not match'
''.format(left.range, right.range))
if isinstance(left.range, Field) and \
Copy link
Member

Choose a reason for hiding this comment

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

Please no backslash line continuation. To split across lines, use parentheses:

if (isinstance(left.range, Field) and
        isinstance(right.range, Field)):
    # do stuff

odl/set/sets.py Outdated
@@ -310,6 +310,25 @@ def field(self):
"""
return self

def contains_set(self, other):
raise NotImplementedError('field {!r} does not have method '
'`contains_set`'.format(self))
Copy link
Member

Choose a reason for hiding this comment

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

This error would already be raised by the parent, wouldn't it?

odl/set/sets.py Outdated
return other
else:
raise ValueError('Fields {!r} and {!r} do not include '
'each other'.format(self, other))
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 skeptical towards this. It doesn't fulfill the rules for set addition, and just as a shortcut for something used in OperatorSum, it shouldn't introduce an API change for Field.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I am adding fields, not sets and so I thought I went with the Minkowski addition A + B = {a + b | a in A, b in B} (which at least for the implemented numbers is trivial). I actually haven't thought about other sets.

What would be your suggestion to approach this? At first I did something like if isinstance(a.range, ComplexNumbers) or isinstance(b.range, ComplexNumbers): return ComplexNumbers(); elif ... but that was rather unpleasant to look at...

Copy link
Member

Choose a reason for hiding this comment

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

I am adding fields, not sets and so I thought I went with the Minkowski addition A + B = {a + b | a in A, b in B} (which at least for the implemented numbers is trivial). I actually haven't thought about other sets.

I was thinking about that definition as well. My point is that IMO, it makes little sense to implement an operation that would be valid (but hard to implement in general) for any Set only for a subclass. So anything short of extending Set, or at least Field, in a generic way is not really desirable. Sure, if all we are dealing with are either real or complex numbers, adding them is trivial. But what if someone decides to add a finite field class -- adding them might not even produce a new field, depending on the definition of addition.

The main issue is that the means you chose are way too generic and sweeping for the problem you're trying to solve. Assuming that this problem actually needs to be solved (which it might not, see my other comments), the least invasive solution would be to introduce an ordering of the sets you're looking at and then just use the max function:

# At module level
_SET_ORDER = {Integers(): 1, RealNumbers(): 2, ComplexNumbers(): 3}

# Somewhere else
field = max(field1, field2, key=lambda f: _SET_ORDER[f])

# Or making a function (only for internal use)
def _largest_field(*fields):
    _SET_ORDER = {Integers(): 1, RealNumbers(): 2, ComplexNumbers(): 3}
    return max(fields, key=lambda f: _SET_ORDER[f])

But this would be just a workaround, really.

Copy link
Contributor Author

@miallo miallo Apr 30, 2018

Choose a reason for hiding this comment

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

[look at the comment below first]
Can you think about this for a moment:
Let us define addition as (a_1⊕b_1)+(a_2⊕b_2)=(a_1+a_2)⊕(b_1+b_2) and (a_1⊕b_1)*(a_2⊕b_2)=(a_1*a_1)⊕(b_1*b_2) (a_i in A, b_i in B). Let the order of A be larger than the order of B. Therefore the mapping (a,b)->a should (according to my calculations for GF(2)⊕GF(3) (yes, I wrote down the tables 😂) and some thoughts) lead to a valid field again. And since all finite fields are isomorphic to Galois Fields then the implementation of contains_set(self, other) for them should just return self.order >= other.order and the addition works perfectly.

From my point of view this would be the "obvious" addition for fields.

Copy link
Contributor Author

@miallo miallo Apr 30, 2018

Choose a reason for hiding this comment

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

The point I was trying to make is that I think, that this diagram commutes. And F is also a homomorphism if the order of the isomorphic Galois field of A is larger, than of B:
diagram
(But the days of my maths classes are all long gone)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I see... The spaces are homomorph but not isomorph. You were right all along...

Copy link
Member

Choose a reason for hiding this comment

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

🤣 You went on quite a journey there, my mention of finite fields was more like an example.. In some sense you proved my minimal point that it's not trivial to implement set addition and keep it in line with mathematics for arbitrary fields.

But it doesn't seem to be necessary here anyway.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I know. But it was fun to try to remember all the maths stuff I learned 5 years ago and never used again. And despite (or maybe because) studying physics I sometimes enjoy to do "recreational maths" 😂

raise OpTypeError('operator ranges {!r} and {!r} do not match'
''.format(left.range, right.range))
else:
range = left.range
Copy link
Member

Choose a reason for hiding this comment

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

Actually, I'd prefer not to have any clever mechanism for "fixing" cases that "obviously" should work, like

L1Norm(cn(2), range=ComplexNumbers()) + L1Norm(cn(2), range=RealNumbers())

IMO this should fail, and users should explicitly have to use RealPart or ComplexEmbedding to state their intentions (of course, for that, those operators need to work with RealNumbers and ComplexNumbers, resp.).

''.format(left.range, right.range))
if isinstance(left.range, Field) and \
isinstance(right.range, Field):
range = left.range + right.range
Copy link
Member

Choose a reason for hiding this comment

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

Just as an aside, you can make this a one-liner without implementing __add__:

range = max(left.range, right.range, key=lambda s: 0 if s == RealNumbers() else 1)

Slightly hackish since it assumes that else means ComplexNumbers, but it works 😉 .

"""
if range is None:
range = RealNumbers()
Copy link
Member

Choose a reason for hiding this comment

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

How about

def __init__(self, domain, exponent, range=RealNumbers()):
    ...

? I feel that the current variant is unnecessarily indirect. None is usually a proxy for some_value_that_depends_on_other_stuff_or_must_be_computed

Copy link
Member

Choose a reason for hiding this comment

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

This only works if we want RealNumbers to be forever immutable. Are we sure about that? (Id be fine with it)

Copy link
Member

Choose a reason for hiding this comment

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

Yeah, we can make that clearer by adding __slots__ to RealNumbers and ComplexNumbers.

Domain of the functional.
outer_exp : {1, 2, inf}, optional
Exponent for the outer norm.
singular_vector_exp : {1, 2, inf}, optional
Exponent for the norm for the singular vectors.

range : `Field`, optional
Range of the functional. Default: `domain.field`.
Copy link
Member

Choose a reason for hiding this comment

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

Blank line needed here

Functional.__init__(self, space=func.domain)
# TODO: Implement grad_lipschitz? If func.is_linear:
# grad_lipschitz = func.grad_lipschitz * ||constant||_\infty
# else: grad_lipschitz = np.nan
Copy link
Member

Choose a reason for hiding this comment

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

Nah, I wouldn't bother. The odds are that grad_lipschitz will be removed.

@@ -929,8 +948,10 @@ def __init__(self, left, right):
raise TypeError('`func` {} is not a `Functional` instance'
''.format(right))

# HELP: What should happen if left.range != right.range?
Copy link
Member

Choose a reason for hiding this comment

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

IMO either nothing, or raise early whatever error adding the convex conjugates would raise. I tend towards not doing anything and letting convex_conj fail, at least if the error message is good enough to tell users what exactly went wrong.

Another way to play this would be to strictly require RealNumbers as range for both, since from the mathematical definition it wouldn't make sense otherwise. However, as in the convex conjugate there's probably an adapted definition for complex-valued functionals, so we probably shouldn't close the door here.

@@ -1151,10 +1171,11 @@ def __init__(self, left, right):
if not isinstance(right, Functional):
raise TypeError('`right` {} is not a `Functional` instance'
''.format(right))
range = left.range + right.range
Copy link
Member

Choose a reason for hiding this comment

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

See earlier comment, let's not add such cleverness.

@@ -1211,14 +1232,14 @@ def __init__(self, dividend, divisor):
raise TypeError('`divisor` {} is not a `Functional` instance'
''.format(divisor))

if dividend.domain != divisor.domain:
raise ValueError('domains of the operators do not match')
range = dividend.range + divisor.range
Copy link
Member

Choose a reason for hiding this comment

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

Same here

return 0.0
else:
return float(getattr(inp, 'real', inp))
Copy link
Member

Choose a reason for hiding this comment

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

Besides the implementation, I'm wondering if we really want this implicit conversion. I'd actually prefer to let this fail for objects that don't implement __float__ but do implement real.
There's also a weird side effect for objects that implement both __float__ and real, namely that real is preferred over __float__. That doesn't seem quite right to me.
Other opinions, @aringh, @adler-j?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I just needed something to convert the complex numbers to real ones. You could instead have elif dtype(inp)==complex: return inp.real; else: return float(inp). Or do you have a different idea?

Copy link
Member

Choose a reason for hiding this comment

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

Yes, my point was that IMO this shouldn't happen in the first place. If you call RealNumbers().element(1 + 1j) it should fail rather than taking the real part.

To me the better approach would be to make the RealPart and ImagPart operators (and their friends) work for RealNumbers and ComplexNumbers, respectively. Likely this requires little less than implementing .real and .imag on those classes, but that's not a guarantee 😉

Copy link
Member

Choose a reason for hiding this comment

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

In this last case, should ImagPart(RealNumbers().element(1)) give 0 or should it fail? The latter represents seeing RealNumbers and ComplexNumbers as two different entities, while the former somehow implicitly assumes that the real numbers are seen as embedded in the complex numbers, just with imaginary part 0. Which one is preferred?

Copy link
Member

Choose a reason for hiding this comment

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

I'd take Python (and NumPy) as a guide here:

>>> x = 1.0
>>> x.imag
0.0
>>> x = np.array(1.0)
>>> x.imag
array(0.)
>>> float(1 + 1j)
...
TypeError: can't convert complex to float
>>> np.array(1 + 1j, dtype=float)
...
TypeError: can't convert complex to float

Copy link
Member

Choose a reason for hiding this comment

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

Some quick trials gives:

>>> (1).imag
0

Hence, ImagPart(RealNumbers().element(1)) should return 0.

functionals = [functionals[0]] * functionals[1]
for func1 in functionals:
Copy link
Contributor Author

Choose a reason for hiding this comment

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

I guess

from functools import reduce
...
range = reduce((lambda x, func: x if func.range.contains_set(x) else func.range), functionals)

would be a more elegant way, but needs the import. Any comments?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Now that I think about it, I guess this is the way to do it (because it doesn't need the import):

range = RealNumbers()
for func in functionals:
    range = func.range if func.range.contains(range) else range

(Or maybe set range = None in the beginning?)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Or maybe it is reasonable to implement a method for the unification of sets? Something like

class Set {...
    def unified_with(self, other):
        """This is the default implementation for the unification of two sets
        It only supports equal sets
        """
        if type(self) == type(other):
            return other
        else:
            raise NotImplementedError('Cannot get the unified field for {} and'
                                      '{}'.format(self, other))

And the Number implementation

    def unified_with(self, other):
        """Returns the unified Field of `self` and `other`"""
        return other if other.contains_set(self) else self

That way new fields could easily implement their own method.
If we would end up doing so you might want to discuss the name...

Copy link
Member

Choose a reason for hiding this comment

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

I'd suggest the simplest possible solution:

range = RealNumbers()
for func in functionals:
    if func.range.contains(range):
        range = func.range

@@ -182,7 +182,7 @@ class LinCombOperator(Operator):
LinCombOperator(a, b)([x, y]) == a * x + b * y
"""

def __init__(self, space, a, b):
def __init__(self, space, a, b, rangeType=None):
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Oh, that wasn't meant to be pushed... I used it for a mapping-oparator from R2 to C...

@kohr-h
Copy link
Member

kohr-h commented Jun 15, 2018

Merged from odl main branch

Two things:

  1. We don't do merge from master, please rebase instead. That keeps a nice and clean history. (If you need to know how it works, just ask.)
  2. We have some guidelines regarding commit messages that we try to follow ourselves (though sometimes we get sloppy, too 😉). It makes the list of commits much more digestible if you see at a glance which commit does what.
    Short version:
    • Prefix with BUG:, MAINT:, STY:, TST: etc.
    • Use imperative style ("Fix bug", not "Fixed bug")
    • Fit first line in ~50 characters, then go wild after a blank line

I guess this PR is gonna end up as one big squashed commit since it's been going a bit back and forth (no issue with that), so we'll likely fix it then.

@adler-j
Copy link
Member

adler-j commented Jun 28, 2018

So what's the status of this? I'd love for this to get finalized!

@miallo
Copy link
Contributor Author

miallo commented Jul 12, 2018

So what's the status of this?

As I mentioned somewhere before: I have never written any tests in python and I would rather have someone else do them since I am not 100% confident with what is to be expected for the default values (expecially for things like IndicatorLpUnitBall.convex_conj). Otherwise I hope I found most of the 'obvious' bugs for now...

@adler-j
Copy link
Member

adler-j commented Jul 12, 2018

I'm currently very busy leading up to the vacation, but I can do that after summer. Thanks again for the contribution!

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.

5 participants