Since the measurement operator and its adjoints should be functions, we can't always use Eigen's CG.
A first implementation might just wrap it though.

• implement L1 proximal algorithm #63
• implement algorithm in pdf below

Carrillo et al (2015) - Why CLEAN when you can PURIFY? [BASP2015 presentation].pdf

aka sopt_l1_solver_padmm in sopt_l1.c

• that fft one
• with sara rather than a single wavelet family

• Framework to access previous version of sopt
• sdmm in-painting test
• that other test with an FFT

f(x, γ, Ψ) = min_{z} 0.5||x - z||_2^2 + λ||Ψ z||_1

aka sopt_prox_l1 in sopt_prox.c

This is just to make sure we can instanciate the same algorithms as in the original sopt.
Needed to finish off #30.

Some sort of framework in cmake so that we can easily test the new code vs the old code. It might
mean copying some of the C files to the test directory and keeping duplicates for now. Or adding
some way to download an old commit of the code, compile it and run against that.

At present the l1 norm is recomputed outside of the l1 prox call. This could be avoided by having the l1 prox function return the last l1 norm, therefore avoiding the need to reapply Psit to the current image.

• sopt_utility_projposr: fancy speak for A = A.real; [A.real <= 0] = 0e0
• sopt_utility_l1norm: sum w_i ||x_i||
• sopt_utility_softth: sgn(x) * max(0, ||x|| - threshhold)
• optionally, that sampling measurement thingie: sopt_meas_urandsamp
• CG: see #26

• Google
• Hayai
• Others...

spdlog doesn't look too bad.
Should be optional, though.

None of the files have copyright statements for now.
We could add one, although I'm not sure who holds the copyright (basp-group?)
The copyright could be updated to GPL3
An AUTHORS file could be added

Copy-pasta from purify job which use JUnits XML format spat by pytest. Catch has something as well. Need to check.

Catch has a fairly rapid development cycle. So checking MD5
when downloading the file fails unconveniently often. Since this is for tests only, I'll try to use
some other method to verify download.

Implementation does not care what the $g_i$ objective functions are. L1 will be one of many possible
instanciation when creating the SDMM operator. So not a good name.

Right now, I'm just using iter < itermax and a user specified functions.
But maybe we could use the $z ~ \sum_i x_i - prox(x_i)$ + user function.

It should mean less copies of inputs and outputs

in the file sopt_l1.c, in padmm, we have:

// Copy r in such a way so that as it can be read as real data
// in the real-complex case (xsol used for temporary storage).
if (param.real_out == 1 && param.real_meas == 0) {
for (i = 0; i < nx; i++){
*((double*)xsol + i) = creal(*((complex double*)r + i));

*((complex double*)r + i) = 0.0 +0.0*I;
}
for (i = 0; i < nx; i++)
*((double*)r + i) = *((double*)xsol + i);
*((double*)xsol + i) = 0.0;
}

The second to last line should likely be in the last for-loop.

Since we've added and are adding a number of new solvers, we probably need to develop a consistent naming approach. For example, we currently have sopt_l1_solver (Douglas-Rachford), sopt_l1_sdmm (SDMM) and sopt_l1_solver_aadmm (approximate ADMM), and also reweighted versions. In future, we might also like to support the unconstrained, as well as the constrained problem (i.e. + lambda * data fidelity term).

Renaming will be a bit of a pain, since it would be a breaking change, but it's probably best to sort out now rather than later (we can also add in old interfaces that could be depreciated in time).

How about sopt_l1_solver_? For example, sopt_l1_solver_constrained_aadmm_rw or sopt_l1_solver_unconstrained_sdmm. Or is that convention too heavy?

@rafael-carrillo @mdavezac @vijaykartik Thoughts? Suggestions?

For instance, SARA returns a larger vector

General framework seeks to minimize $\sum_i g_i(L_i*x)$. It takes as input a list of operators
consisting of:

• L_i
• L_i^\dagger
• the proximal of g_i

It might/can also take as input a conjugate gradient method, a convergence criteria, itermax, gamma,
... something else?

Small bug on google/gbenchmark prevents compilation.
Once google/gbenchmark#152 is fixed, we should move back to tracking it rather than my fork.

Required for #38

FISTA only helps with many sub-iterations are required. For the inexact ADMM we look for inexact solutions so we don't need to run so many sub-iterations. We could get a performance saving by turning FISTA off since it requires a lot of additional memory.

Some general way to handle random numbers.
Right-now, only used in sampling operator.
It might be the sampling operator should take an argument for the random number generator.

In [2]: import wavelets.wavelets as wv
In [3]: import numpy as np
In [4]: input = np.random.random(128)
In [5]: wv.dwt(input,"xyz",1)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-5-71b4d0f139bc> in <module>()
----> 1 wv.dwt(input,"xyz",1)
wavelets.pyx in wavelets.dwt (python/wavelets/cython_wavelets.cc:2989)()
wavelets.pyx in wavelets._rdwt (python/wavelets/cython_wavelets.cc:2173)()
RuntimeError: std::exception

That really does not explain what went wrong, neither in python nor in C++.
For #38 or later.

Probably some in-painting example with random images of different sizes, and different wavelets
basis.

Currently WSClean exposes these functions:

void wsclean_initialize(void** userData, const struct purify_domain_info* domain_info, struct purify_domain_data_format* data_info);
void wsclean_read(void* userData, DCOMPLEX* data, double* weights);
void wsclean_write(void* userData, const double* image);
void wsclean_deinitialize(void* userData);
void wsclean_operator_A(void* dataIn, void* dataOut, void* userData);
void wsclean_operator_At(void* dataIn, void* dataOut, void* userData);

For neatness, I suggest two more changes. Maybe not really necessary for now during the hack week, but nice to think about in the future:

• The initialize method can return function pointers to the other functions. That way, the only function call you need to know about is the initialization function. These functions pointers can be stored in what is currently the 'purify_domain_data_format' structure, but maybe rename it. So, Purify would do something like this:
  wsclean_initialize(&userdata, &dinfo, &formatinfo);
  read_function=formatinfo.read_function
  deinitialize_function=formatinfo.deinitializate_function
  ...
  (_read_function) (userdata, parameters)
  ...
  (_deinitialize_function) (userdata)
• The measurement operator functions (for wsclean; wsclean_operator_A and wsclean_operator_At) should be fully typed, be const correct and for consistency start with the userdata parameter:
  void wsclean_operator_A(void* userData, const double* dataIn, double complex* dataOut);
  void wsclean_operator_At(void* userData, const double complex* dataIn, double* dataOut);
  Purify should cast the function pointers when necessary to e.g. support double->double measurement operators, instead of casting the parameters everytime in both the callee and caller.

Means less maintenance for sopt proper.

• Catch
• gbenchmark
• spdlog

$ make
make: *** No rule to make target `../sopt/src/c/sopt_error.o', needed by `../sopt/lib/libsopt.a'.  Stop.

probably it should be removed.

In purify example, changing the logging level in initialized does not seem to work...

At this juncture there are still further enhancements that might be beneficial to the (refactored)
wavelet code:

• remove const_cast and use rvalue reference overloads instead. const_cast is the canonical
  Eigen way of doing things, but it is quite ugly.
• 2d and 1d code could be consolidated if the 2d operations are always colwise and rowwise (which I
  think they are)
• periodic_scalar_product could be made to use eigen functions only, without an explicit loop.
• nominally, we only need as much memory as N by rows, where N is the number of wavelet
  coefficients.

Tested with L_i = Id for simplicity. Now need to test general case.
Blocks #30.

It should be possible to install sopt to disk somewhere. And it should register itself so other
CMake projects can find it.
This is done in the "old" version. Needs to be updated for refactoring.