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

Issue 59 coerce #95

Closed
wants to merge 5 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 9 additions & 12 deletions odl/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -33,26 +33,23 @@
# Propagate names defined in __all__ of all submodules into the top-level
# module

from . import diagnostics

from . import discr
from .discr import *
__all__ += discr.__all__
# Import modules in order of dependency
from . import set
from .set import *
__all__ += set.__all__

from . import operator
from .operator import *
__all__ += operator.__all__

from . import set
from .set import *
__all__ += set.__all__

from . import space
from .space import *
__all__ += space.__all__

from . import trafos
from .trafos import *
__all__ += trafos.__all__
from . import discr
from .discr import *
__all__ += discr.__all__

from . import trafos
from . import diagnostics
from . import solvers
4 changes: 2 additions & 2 deletions odl/discr/discretization.py
Original file line number Diff line number Diff line change
Expand Up @@ -270,7 +270,7 @@ def copy(self):
"""Create an identical (deep) copy of this vector."""
return self.space.element(self.ntuple.copy())

def asarray(self, out=None):
def asarray(self, *args, **kwargs):
"""Extract the data of this array as a numpy array.

Parameters
Expand All @@ -279,7 +279,7 @@ def asarray(self, out=None):
Array in which the result should be written in-place.
Has to be contiguous and of the correct dtype.
"""
return self.ntuple.asarray(out=out)
return self.ntuple.asarray(*args, **kwargs)

def __eq__(self, other):
"""``vec.__eq__(other) <==> vec == other``.
Expand Down
54 changes: 50 additions & 4 deletions odl/discr/lp_discr.py
Original file line number Diff line number Diff line change
Expand Up @@ -149,7 +149,7 @@ def element(self, inp=None):
if arr.ndim > 1 and arr.shape != self.grid.shape:
raise ValueError('input shape {} does not match grid shape {}'
''.format(arr.shape, self.grid.shape))
arr = arr.flatten(order=self.order)
arr = arr.ravel(order=self.order)
return self.element_type(self, self.dspace.element(arr))

@property
Expand Down Expand Up @@ -229,14 +229,25 @@ class DiscreteLpVector(DiscretizationVector):
"""Representation of a `DiscreteLp` element."""

def asarray(self, out=None):
"""Extract the data of this array as a numpy array.
"""Extract the data of this vector as a numpy array.

Parameters
----------
out : `numpy.ndarray`, optional
Array in which the result should be written in-place.
Has to be contiguous and of the correct dtype and
shape.

Returns
-------
out : `numpy.ndarray`
This vector as an array, if ``out`` parameter was given, they are
equal.

See also
--------
asflatarray
FnBaseVector.asarray
"""
if out is None:
return super().asarray().reshape(self.space.grid.shape,
Expand Down Expand Up @@ -266,14 +277,49 @@ def asarray(self, out=None):
super().asarray(out=out.ravel(order=self.space.order))
return out

def asflatarray(self, *args, **kwargs):
"""Extract the data of this vector as a 1d numpy array.

Parameters
----------
out : `numpy.ndarray`, optional
Array in which the result should be written in-place.
Has to be contiguous and of the correct dtype.

Returns
-------
out : `numpy.ndarray`
This vector as an array, if ``out`` parameter was given, they are
equal.

See also
--------
asarray
FnBaseVector.asflatarray
"""
return super().asarray(*args, **kwargs)

@property
def ndim(self):
"""Number of dimensions."""
"""Number of dimensions.

See also
--------
TensorGrid.ndim
"""
return self.space.grid.ndim

@property
def shape(self):
# override shape
"""Size in each dimension

For example, a 3d space can have size (3, 5, 10), indicating that it
is 3 elements big in the "x" direction, 5 in the "y" and "10" in the z.

See also
--------
TensorGrid.shape
"""
return self.space.grid.shape

def __setitem__(self, indices, values):
Expand Down
4 changes: 4 additions & 0 deletions odl/operator/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -33,3 +33,7 @@
from . import pspace_ops
from .pspace_ops import *
__all__ += pspace_ops.__all__

from . import embedding
from .embedding import *
__all__ += embedding.__all__
133 changes: 133 additions & 0 deletions odl/operator/embedding.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,133 @@
# Copyright 2014, 2015 Jonas Adler
#
# This file is part of ODL.
#
# ODL is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# ODL is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with ODL. If not, see <http://www.gnu.org/licenses/>.

"""Operators to cast between spaces."""

# Imports for common Python 2/3 codebase
from __future__ import print_function, division, absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import super

from odl.space.base_ntuples import FnBase
from odl.discr.lp_discr import DiscreteLp
from odl.set.pspace import ProductSpace
from odl.operator.operator import Operator

__all__ = ('Embedding', 'EmbeddingFnInFn', 'EmbeddingPowerSpaceInFn',
'EmbeddingFnInPowerSpace')


class Embedding(Operator):
"""An operators to cast between spaces."""

# TODO: is this needed?

def __init__(self, origin, target):
super().__init__(origin, target)


class EmbeddingFnInFn(Embedding):
Copy link
Member

Choose a reason for hiding this comment

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

Small thing: Into instead of In

Copy link
Member Author

Choose a reason for hiding this comment

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

Sure!

""" Embed a `FnBase` space into another """

def __init__(self, origin, target):
assert isinstance(origin, FnBase)
assert isinstance(target, FnBase)
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 wondering if we could replace all our obvious type checks with asserts? At least when the docstring says what you can provide, the error message is easy to relate to the error source.

Copy link
Member Author

Choose a reason for hiding this comment

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

I thought about this before, but I'm not quite sure. Anytime a user could give the wrong input we should give an error imo. Asserts are to make sure we havnt made any algorithmic error.

Copy link
Member

Choose a reason for hiding this comment

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

I just tried it in the console. assert is really very short on information, you almost never get away without looking into the code. Better proper messages.


super().__init__(origin, target)

def _apply(self, x, out):
out[:] = x.asflatarray()

def _call(self, x):
return self.range.element(x.asflatarray())

@property
def adjoint(self):
return EmbeddingFnInFn(self.range, self.domain)


class EmbeddingDiscreteLpInDiscreteLp(Embedding):
""" Embed a `DiscreteLp` space into another """

Copy link
Member

Choose a reason for hiding this comment

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

This operator could also handle extension by zero (constant? nearest value?) to a larger domain.

Copy link
Member Author

Choose a reason for hiding this comment

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

Only with some strict constraints, for example the values need to coicide in some points at least approximately if we want to do this on a dspace level. Another option would be to walk the long walk with restriction and extensions, but I fear for the performance.

Copy link
Member

Choose a reason for hiding this comment

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

On a dspace level, only some restricted operations are possible, that's for sure. The best you can hope for is that your old data stays together as a chunk in the new flat array.

But in general, it's hard to get around using restriction and extension, and then one can implement a much wider range of operations going beyond pure embedding, e.g. subsample by interpolation. With vectorization, this should not be too inefficient. The good part here is also that the interpolation comes along with some options for how to handle values outside.

Conclusion: alright, keep it simple and implement the fancy stuff in a different operator.

Copy link
Member Author

Choose a reason for hiding this comment

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

Agree.

def __init__(self, origin, target):
assert isinstance(origin, DiscreteLp)
assert isinstance(target, DiscreteLp)
assert origin.grid.shape == target.grid.shape
Copy link
Member

Choose a reason for hiding this comment

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

This should only be possible if target.exponent <= origin.exponent for mathematical reasons.

Copy link
Member Author

Choose a reason for hiding this comment

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

Not sure about that, embedding (to me) means, "Do as your told and ignore the consequences". Kindoff like how mathematicians often "identify" things with others. Perhaps the vectors given to this will only be in a subspace where it still works? I think we need a way for users to "force" this through.

Copy link
Member

Choose a reason for hiding this comment

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

This guy could take a parameter strict=True which would check the mathematical condition. To me, an embedding is the identity mapping from a space to a superspace, not simply closing your eyes. If users want to go for it, they can set strict=False. We can also call it force=False.

Copy link
Member Author

Choose a reason for hiding this comment

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

This guy could take a parameter strict=True which would check the mathematical condition.

Sure, but I'd default it to False.

To me, an embedding is the identity mapping from a space to a superspace, not simply closing your eyes.

Then this needs to be renamed. Very few of these are actually to superspaces, and the adjoint would then never be defined (which obviously isn't that useful).

Copy link
Member

Choose a reason for hiding this comment

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

Then this needs to be renamed. Very few of these are actually to superspaces, and the adjoint would then never be defined (which obviously isn't that useful).

It's a projection then. In that case I'm good with not checking the exponent.


super().__init__(origin, target)

def _apply(self, x, out):
out[:] = x.asarray()

def _call(self, x):
return self.range.element(x.asarray())

@property
def adjoint(self):
return EmbeddingDiscreteLpInDiscreteLp(self.range, self.domain)


class EmbeddingPowerSpaceInFn(Embedding):
""" Embed a `PowerSpace` of `FnBase` space into a `FnBase` """

def __init__(self, origin, target):
# TODO: tests goes here

super().__init__(origin, target)

def _apply(self, x, out):
out[:] = x.asflatarray()

def _call(self, x):
return self.range.element(x.asflatarray())

@property
def adjoint(self):
return EmbeddingFnInPowerSpace(self.range, self.domain)


class EmbeddingFnInPowerSpace(Embedding):
""" Embed a `FnBase` into `PowerSpace` of `FnBase`'s"""

def __init__(self, origin, target):
# TODO: tests goes here

super().__init__(origin, target)

def _apply(self, x, out):
index = 0
for sub_vec in out:
sub_vec[:] = x.asflatarray(index, index + sub_vec.size)
index += sub_vec.size

def _call(self, x):
out = self.range.element()
index = 0
for sub_vec in out:
sub_vec[:] = x.asflatarray(index, index + sub_vec.size)
index += sub_vec.size
return out

@property
def adjoint(self):
return EmbeddingPowerSpaceInFn(self.range, self.domain)


if __name__ == '__main__':
from doctest import testmod, NORMALIZE_WHITESPACE
testmod(optionflags=NORMALIZE_WHITESPACE)
Loading