odlgroup / odl Goto Github PK

All the cool kids seem to have a logo of some sort, we need one aswell. Any input on what would look cool would be cool.

We need one big logo (to use in presentations etc). And one small logo to use in webpage thumbnails and such.

We need a public documentation page. This would include both the odl programming doc, but also issues #2, #4, #5. First we need to decide on where to put this documentation, some options are

• Use Read the docs
• Create a proper odl webpage
• Use GitHub

Since we have decided to use Sphinx for the documentation, all of these are ofc compatible, and we could for example use all of them. For maintainance, I propose that we try to keep this at a minimum, at least until we get stable releases.

As discussed in #21, we need the following two operators in default_ops

def ConstantOperator(odl.Operator):
    def __init__(self, x):
        self.x = x

    def _call(self, _):
        return self.x

def ResidualOperator(odl.Operator):
    def __init__(self, op, x):
        self.op = op
        self.x = x

    def _call(self, y):
        return self.op(y) - self.x

We need a build server that runs windows.

Let's try to find out if there is a good and predictable way to move documentation from _apply() and _call() of Operator implementations to the respective __call__() method so it appears in the right place in the final doc. In principle, a script would for each subclass of Operator

• take the respective cls._apply.__doc__
• copy it to cls.__call__.__doc__
• modify the line out : something to out : something, optional
• add the note that if out is provided as argument, it is returned
• replace cls._apply.__doc__ and cls._call.__doc__ with a hint not to use the functions directly but rather call operator(x) or operator(x, out), respectively

The question is where to put this functionality. Maybe one could implement it in Operator.__new__()?

When reading the discussion in #49, I realized that we could use a "project on subspace" operator

class Projection(odl.Operator):
    def __init__(self, space, indices):
        assert isinstance(space, ProductSpace)
        self.indices = indices
        super().__init__(space, space[indices], linear=True)

    def _call(self, x):
        return x[self.indices]

I have somewhat lost track about what expression in which combinations of operators, scalars and vectors are currently supported and what they mean mathematically.

I suggest that we write up a documentation chapter about operator expressions, possibly summarized as a table containing (example: operator * scalar)

• The expression: op * scal
• What it means mathematically: x --> op(scal * x)
• Conditions: scal in op.domain.field
• The created expression class: OperatorRightScalarMult

I suggest that we do the following:

We make odlpp a proper sub-package. I.e. we move all odl/space/cu_ntuples and test/space/cu_ntuples out of this repo and into the current odlpp repo.

We then change odlpp to odl-cuda

Finally we make odlpp a Namespace package. What this means is that the users can install it separately, and it will appear just as now under odl.space.cuda.

What do we think?

I think that in the long run, the solvers need more structure than can be handled in a single file. They will probably also add dependencies on third-party software which we don't want in the core ODL package. So unless there are strong reasons against it now, I'll move the code to the odl-solvers repository.

There are a couple of commit messages like "merge" or "bugfix" in the commit history which are actually worse than no message at all. We need to become better at this if we want our Git history to be searchable and useful.

As a very good reference, one can once again at how NumPy does it. I suggest that we follow the same nomenclature, at least in the final versions of pull requests. Intermediate commits on branches may have bad messages, as long as they are cleaned up with git rebase -i before the pull request is sent (or at the very latest before it is merged into master).

Short summary:

• Start each commit message with a capitalized acronym indicating its purpose, e.g. BUG for a bugfix (see the table in the link).
• Write a short summary (max ~50 characters) on the first line.
  • Use imperative style, i.e. it should be possible to append the message (w/o acronym) to "When applied, this commit will ...". In particular, don't write in past tense (add feature xyz instead of added feature xyz). A very simple reason for this style is that you save on average 2 characters for something else.
    • Be as precise as possible (needs thinking).
    • Don't use markdown.
• Optional: write an extended summary if additional information is needed.
  • Separate it by a blank line from the first line.
  • Explain why the commit is necessary, e.g. detail the problem which is fixed.
  • You can use markdown here, e.g. include item lists, tables, links, @mentions, ...
  • Restrict the lines to 72 characters in length.
• Refer to tickets (issues) if existing, e.g. "See #2", "Closes #2" (this will actually close the issue) or similar

Beside the pure implementation, we need to review one design decision, namely where the discrete space functions are specified. Currently, this is completely done in the data space Fn by adding a weight to the inner product, norm and dist.
If we decide to keep it that way, we must add an exponent parameter to the weighting classes, and for p != 2, things become slightly more complicated:

• inner must be disabled
• norm must be re-implemented, e.g. in the matrix-weighted case, the matrix power M^(1/p) has to be (pre-)computed since the weighted norm is |||x|||_p = ||M^(1/p)x||_p. Probably this is not feasible for larger problems, we'll see. For vector- and constant-weighted spaces, it's easy though.

I added intersphinx so that we can link documentation against numpy, it does however seem like Read The Docs has a old version of intersphinx that crashes, we should somehow fix this.

This has started to fail:

================================== FAILURES ===================================
_________ [doctest] odl.discr.discr_mappings.RawGridCollocation._call _________
245
246         >>> coll_op = RawGridCollocation(funcset, grid, rn, order='F')
247         >>> coll_op(func_elem)
248         Rn(6).element([-2.0, -1.0, -3.0, -2.0, -4.0, -3.0])
249
250         Use all free variables in functions you supply, otherwise
251         the automatic broadcasting will yield a wrong shape:
252
253         >>> func_elem = funcset.element(lambda x, y: 2 * x)
254         >>> coll_op(func_elem)
Expected:
    Traceback (most recent call last):
    ...
    ValueError: input shape (2,) not broadcastable to shape (6,).
Got:
    Rn(6).element([2.0, 4.0, 2.0, 4.0, 2.0, 4.0])

E:\Github\ODL\odl\discr\discr_mappings.py:254: DocTestFailure
==================== 1 failed, 295 passed in 4.28 seconds =====================

My suggestion is that we should not have exception tests in a doctest.

See odlgroup/odl-solvers#7 and odlgroup/odl-solvers#9.

This file is autogenerated and is way to long, we should be able to clear up a lot of these things.

From #36 it seems like the cg method (found here) can get a division-by-zero error.

I think (but I'm not sure) that the problem is that if to many iterations are used, the exact solution is found. In this case, in the next iteration the search direction p found is the zero-direction, and thus p.inner(Ap) is zero, which is used in the code as alpha = sqnorm_r_old / p.inner(Ap).

Keeping make html warning free enables us to notice when we get new warnings. We should thus try to make it warning free. I am currently on my way to do this.

I tired to push to the repo and got the error message:
remote: Permission to odlgroup/odl.git denied to aringh.
fatal: unable to access 'https://github.com/odlgroup/odl.git/': The requested URL returned error: 403

Do I miss some rights, or is there somethings else wrong?

It's become relatively common that we check in things that break existing tests and examples to the main branch. This is OK right now, but as we proceed this will not be feasible.

I therefore propose that we (in some way) create a build server. This could preferably be run on some decently good computer at KTH, and should preferably run the tests with different configurations (Windows, Linux, python 2, python 3).

Solvers are currently absent from docs.

We should gather common questions and put together a FAQ. I started one in 5fe6d3a.

With :automodule: it seems class magic methods, most notably __init__ does not get documented. This causes some obvious issues.

The docstrings in the __init__.py files do not end up in the Sphinx documentation. I think they should, since they are a good place to put an introduction to the subpackage contents.

In contrast, I'd argue that one only writes a minimal docstring at the top of a module so that one does not have to scroll all the way down to see the code.

The off-the-shelf gitwash documentation is a good start, but some aspects the work flow are a bit unspecific or missing. We need to adapt those sections, e.g. coding style, commit messages.

After the discussions we had about functionals and the properties they have in the light of (convex) optimization, the question is how to solve this in code. I'll make a quick list of what is needed, feel free to edit:

Suppose phi:X --> R is a functional, x in X

• grad phi(x) = element in X. This is special to Hilbert spaces - in Banach spaces it can be identified with, but is not equal to, a dual space element. Implementation-wise, gradient will probably be an Operator from the functional domain to itself.
• Hessian H_phi(x) = linear functional X --> R. Here it becomes more interesting: should it be something bi-linear (x, y) --> H_phi(x, y) in R or rather an Operator mapping a domain element to a linear functional?
  Update: probably better to regard it as the derivative of the derivative.
• Fenchel conjugate phi* = functional X* --> R. We need to know that phi maps to R. This can also be defined on general Banach spaces, probably harder to find closed forms though. Unless we introduce the notion of a dual space, we can only support Hilbert spaces (or maybe a few other cases where the dual can be identified with a simpler space, e.g. (L^p)* ~ L^q).
• Lagrange (Fenchel-Lagrange?) conjugate phi^ = functional in X?
• Bregman distance?

The fact that we need the real numbers as range would cause slight headaches in the current Operator framework, and most definitions only seem to be relevant for functionals. So I don't really think it makes sense to add this stuff to the Operator class. Better to make a subclass which gets these additional attributes.

Copied from odl-solvers issue 5

I suggest that we write a common class for step length calculation which tells the algorithm what it needs (which part of the history), such that an algorithm can determine which iterates and gradients (what else?) to keep.

A simple approach would be to give the StepLength class attributes

• iter_idcs (iterates)
• func_idcs (function evaluations)
• grad_idcs (gradient evaluations / approximations)
• hess_idcs (Hessian evaluations / approximations)
• step_idcs (step lengths)

which all return lists of indices of items to keep. An algorithm could manage lists of all these objects and pass its lists (possibly kept in a dictionary) to the step length call, which in turn does its thing to calculate the step length.

Just a thought, currently many numpy users do things like:

a = x.dot(A.dot(x))

In odl this would be

a = x.inner(A(x))

or using the Vector.T syntax

a = x.T(A(x))

do we want to somehow support the numpy syntax? This could also allow odl to be used with duck-typing in cases where a code expects a numpy array.

Currently the following gives you 50 figures:

for i in range(50):
    #do stuff
    x.show()

a syntax like

fig=None
for i in range(50):
    #do stuff
    fig=x.show(figure=fig)

or possibly something better, would be very nice.

Since operators like the ResidualOp in the solver tests probably occur rather often, I suggest that we support adding an operator and an element in its range to create a new AffineOperator (may suggest linearity) or ShiftedOperator whose derivative is the one of the original operator, of course.

As noted in implementing #53 ProductSpace.__getitem__ currently returns a list if indices is a slice:

return self.spaces[indices]

In the case where indices is a slice, ProductSpace.__getitem__ should do:

return ProductSpace(self.spaces[indices])

to be sure, this should fail if the space has a custom norm (since this does not transfer trivially).

if indices is a int, we should simply

return self.spaces[indices]

We need some advice for people who want to contribute but do not exactly know how. Intended audience are users and developers. We should describe

• How to best get acquainted with the code
• How to work with Git[Hub] - fork, branch, pull request
• How we would like contributed code to look like

As seen in this CI log the optimizer test occasionally fails.

In the setup.py there is a version number given. I think it would be useful, especially later on, to have this version number also explicit in the README file (or somewhere else, easily accessible).

I know that for large problems this is a BAD idea. However, for small test problems this might be a good things, especially to check things like conditioning of the problem (in my case it would be a 510x510 matrix, so it should be possible to handle).

What lead to the thought was a discussion with @ozanoktem. I have a nonlinear operator and want to solve f(x) = 0. I'm trying to do this with Newton's method, however the CG used there does not seem to find a solution to gradf(x)p = -f(x), although gradf(x) is a positive definite operator for all x. It would therefore be useful to be able to investigate the matrix, to see if it is ill-conditioned.

We need a way for users to install odl without needing to point them towards git.

We need something that users can look at to see where to get started, similar to numpys quickstart

Lets make a short list of what this should contain (obviously inspired by numpy)

• A short note on pre-reqs, mainly some programming knowledge and a note that we try to follow numpys conventions where possible, also some words on the math needed and some links.
• Introduce a LinearSpace, quite obviously Rn, show that you can do what you'd expect with it.
  • Space creation
  • Vector creation, element, zero.
  • Workhorse methods (+, -, /, *), indexing
  • Core methods lincomb, multiply.
  • Show how to interact with numpy. asarray, __array_wrap__ to use ufunc's
  • printing arrays
• Introduce a Set, I'd prefer Interval(0, 1), do some simple things with it.
  • __contains__
• Define L2(Interval(0, 1))
  • create continuous function
  • show function arithmetic
• Discretize l2_uniform_discretization(L2(Interval(0, 1)))
  • Show that norms are mapped as expected etc
  • Show Rn/CudaRn possibility
  • Show complex
  • Show the odl.diagnostics possibilities.
• Define a Operator on the space
  • Introduce Operator
  • Define a very simple operator, perhaps the convolution
  • Show some operator magic (perhaps addition, or scaling)
• Solve an inverse problem
  • I realy like the deconvolution problem here.

The merging of solvers has caused several warnings in the doc, these should be fixed.

Moved from https://github.com/odlgroup/odl-solvers/issues/3

I initialized a structure, but it is obviously far from optimal. We need to find sensible categories for solvers, e.g.

• methods for finding a zero of a goal function(al), e.g. quasi-Newton schemes
• simple methods for minimizing goal functions
• "more complex" methods involving e.g. Fenchel/Lagrange duals and methods combining primal and dual problems. These probably need a more structured approach.

We need to reflect these categories in classes in a sensible manner.

For some reason, we needed the requirements in text files instead of a string in setup.py. What about extra dependencies? For example, the Fourier transform operator will depend on pyfftw, and the wavelet transform will use PyWavelets. I'd want to be able to pull the dependencies by writing

pip install --user -e .[fft,wavelets]

This actually works and asks for the extras_require dictionary for the dependencies, so the dependency info goes there.

Question is: also into a file? Or in a string? Or is it perhaps possible to create extra sections in the requirements file? I could not find anything on the web treating the last topic.

Since we're building a new package, we have the luxury of choosing whether to have unicode literals enabled by default or not. Since they are standard in Python 3, I strongly suggest that we make an attempt at it.

A known issue is the ABCMeta class complaining when unicode_literals are imported. Back then, I tried to solve it by explicitly passing a byte string, but that didn't help. There is, however, good documentation available which may point to an elegant solution of this problem.

If some external module really screws up and cannot be convinced at all to work, it is still possible to remove unicode_literals from imports selectively.

we should add a way for users to determine the version of odl

>>> import odl
>>> odl.ODL_MAJOR_VERSION
0
>>> odl.ODL_MINOR_VERSION
9

There are a few Operator's on product spaces that we could standardize and that would be usefull, see #239 for example.

The ones i can think of are the following (naming obviously TBD).

Distributive

The distributive operator, which calls a series of "sub" operators given the same argument. Think of the matrix:

[A,
 B] 

DistributiveOperator(A,B)(x) = [A(x), B(x)]

Diagonal

The "diagonal" operator, which calls a series of operators "per space", think of the matrix

[A, 0,
 0, B]

DiagonalOperator(A,B)((x,y)) = [A(x), B(y)]

Sum

Finally we can think of calling a sequence of operators with given (different) arguments and somehow reducing the result. Reduction by summation is the obvious choice, but other options could work, think of the matrix:

[A, B]

SumOperator(A,B)((x, y)) = A(x) + B(y)

All of these operators could obviously be generalized to tuples of n operators, for example:

DistributiveOperator(A, B, C, D)(x) = [A(x), B(x), C(x), D(x)]

I think we should add these to odl.operator.default_ops, but first I need some input on the names.

References

The first one ("Distributive") is called vstack in numpy

The second one seems to be scipy.linalg.block_diag

and the last one is hstack

The more I think of it, it would maybe be a good idea to not implement three separate classes for this, but instead to create one base class and then (maybe) have these as special cases. Something like a sparse matrix class (but for operators). That way the user could do complicated shit like

op = GeneralProductOperator([[A   , None, None, B   ],
                             [A   , None, C   , D   ],
                             [A   , C   , None, None]])

which would have the obvious meaning:

op([x, y, z, w]) == [A(x) + B(w), A(x) + C(z) + D(w), A(x) + C(y)]

This is similar to the scip.sparse.bmat method.

We need a buildserver that is configured to run cu_ntuples

The build server should do pep8 checking and report fail if it does not work.

Copied from odl-solvers issue 4

Likewise step length rules, it will be equally useful to define a common interface for termination rules of optimization schemes, provided that there is some common structure.

Somehow 98 % of our issues end up having the proposal label. Slightly pointless. My suggestion is that we follow the commit message label distinction made in NumPy. I'll post it here so we don't have to follow a link:

API: an (incompatible) API change
BLD: change related to building numpy
BUG: bug fix
DEP: deprecate something, or remove a deprecated object
DEV: development tool or utility
DOC: documentation
ENH: enhancement
MAINT: maintenance commit (refactoring, typos, etc.)
REV: revert an earlier commit
STY: style fix (whitespace, PEP8)
TST: addition or modification of tests
REL: related to releasing numpy

Out of these, I'd say we take

• api change
• bug (already there)
• deprecation
• documentation
• maintenance
• style
• testing
• release
• instead of enhancement, new feature or feature request, (more specific)

And I would also vote for labels for different parts of the library,

• core
• solvers
• trafos
• cuda
• ...

I'll label this as a proposal despite the irony.

I'll do a simple test with Travis CI.

Something screwed up the build process of ReadTheDocs. The API part is empty.

According to this travis, the test test_matrix_inner occationally fails.

This is because the result in almost_equal is a float, and our automated precision selector will select a higher precision (since float is 64 bit), but the calculation was actually in 32 bit!

It seems having such a file is not pythonic, instead users should simply do

$ py.test

which should run the "usual tests". We need to update the readme aswell if we do this.