currently we use proxmins gradient implementation instead of the faster LinearFilters. simply using gradients should speed these up

AsinhNorm uses the algorithm from the Lupton et. al 2004 paper but does not include modifications to the algorithm included in the LSST stack and astropy that allow more control over the stretch. So the output using AsinhNorm is

while make_lupton_rgb can give a much more reasonable

with the same beta by altering the stretch of the image.

We should either fix AsinhNorm or slightly adjust our classes so that we can use make_lupton_rgb from astropy in channel_to_rgb. Basically Norm.call would have to accept an argument list *img instead of an array of images.

My preference is to use make_lupton_rgb but if someone else wants to fix AsinhNorm I'm fine with that too.

The code currently introduces a bias when at least one of the PSF image dimensions is even. SCARLET should produce a warning when a PSF with an even shape is used, and potentially add an extra row/column of zero values to make the shape odd.

Currently there is only a warning printed by proxmin when the fit hasn't converged. It would be useful to set an internal flag in such cases.

I think there is something wrong with the codes in 'blend.py' in latest version. That is on the 137th row. There is:

c._gamma(pos, self._img.shape, offset_int=source.center_int) for c in self.components 

While in the gamma definition, there is:

def __init__(self, psfs=None, center=None, dy=0, dx=0): 

The "pos" is correspongding to the "center", but the input shape of dy and dx does not match

There is an error when I am using blend_psf.fit funtion. Concisely, the error is like this:

/export/home/lsst_stack/python/miniconda3-4.2.12/lib/python3.5/site-packages/scarlet-0.2.fd9e460-py3.5-linux-x86_64.egg/scarlet/blend.py in get_model(self, m, combine, combine_source_components, use_sed, flat)
240 model_img = np.zeros((source.K,) + (self._img.shape))
241 for k in range(source.K):
--> 242 model_img[k][source.bb] = model[k][model_slice]
243 return model_img
244

ValueError: could not broadcast input array from shape (3,15,14) into shape (3,15,0)

The input is 400 by 400 image, but it seems that if I use 200 by 200 image, I can successfully run the program. The whole error report is in the attachment.

Error.txt

Scene is an ambiguous word, so it's better to clarify.

In addition, I propose two changes

• Blend and Observation have Frame members instead of being derived from Scene
• Blend.match() happens during initialization to avoid ambiguity of who's responsible for calling it.

When running scarlet on a typical example such as the one provided in the documentation, the script crashes with the following error:

I fixed it to make it execute in my version by replacing line 193 of blend.py by:

converged = res

Traceback (most recent call last):
  File "/Users/remy/Desktop/LSST_Deblending/Test_scarlet.py", line 35, in <module>
    blend.fit(200)
  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scarlet-0.3.b34b32b-py2.7-macosx-10.6-intel.egg/scarlet/blend.py", line 193, in fit
    converged, error = res
ValueError: too many values to unpack

About 1% of my images report a singular matrix error. However, I cannot notice any difference between those images and images that deblend successfully. Detailed bug report and the case are attached.

err.zip

We think that deblender may be overfitting images by including some of the noise belonging to the original image, in the model. Here are two images where we see this effect (with background rms = 10). The residual from after the model has been subtracted from the original image looks a bit smoothed out.

We use metacalibration on the 'Original - Model of Neighbor' image to measure shear and we also calculate any bias in the process. We measured a bias of -0.0215581 +/- 0.00180595 from a set of 1e6 images with background rms = 10. In a control test with 1e6 images and a much lower background rms of 0.001, we measured a smaller bias of -0.00534828 +/- 9.46265e-05. Finally, we also ran another control in which we didn't use deblender and ran metacalibration on images of only a single object in the center. We used 1e6 images and a background rms = 10 for this test and measured a smaller bias of 0.000828241 +/- 0.00177186.

Do you have any suggestions as to how we can prevent this 'overfitting'?

It would be great to provide constraints per objects, e.g. ["MS", "S", None, ...] instead of the same for each objects. This only requires better parsing, no changes to the algorithms.

As part of #78, test if the de/re-convolution can be expressed with simple FFTs instead of a scarlet run to match the PSF kernels.

PR #64 introduced normalization options. Technically these are just another type of constraint, but they are currently handled differently. This issue is meant to correct that in the minimally invasive way.

The normalization is currently only used in the Source initialization, and it only changes the behavior of SimpleConstraint, the only thing that needs to be done is to remove the normalization class member of Source.

This is an extension of #18.

With #29 merged, the size updates are taken from a fixed set of box sizes, which can be configured as argument of Blend. However, Source is still initialized in some way, having some box size. That defies the purpose of having a limited set that can efficiently be cached.

Proposal: Define global variable scarlet.source_sizes with a small helper function to decide which is the best available size. Both is currently part of Blend, but it needs to move up to package level.

Bonus points: Blend carries several other config parameters that change the behavior of the fitting process, but are unlikely to be altered between consecutive instantiations of Blend. I would therefore rather define all of them at the package level. This then cleans up the Blend class.

Thus, add the following to __init__.py:

config = {
  refine_skip=10, 
  source_sizes = [15,25,45,75,115,165],
  edge_flux_thresh=1.,
  exact_lipschitz=False
}

and have Source and Blend load parameters from there.

There a re a lot of things happening when one source is resized. It is clear that the performance suffers from a lot of resizing, but not so clear why.

• resizing all linear operators (gradients, PSFs, etc) for the source
• catching the exception
• reinitializing bSDMM

My suspicion is it's the first item. If the third item is the bottleneck, it's possible doable to rewrite bSDMM to allow for a transparent resizing on the fly. This would also render the exception obsolete.

The documentation should include an example that uses different box sizes to illustrate passing a Config instance to both Source and Blend classes.

This came up during a PR.

In essence, observations should only have psfs=None if there's actually no PSF. We've used this setting as default in cases when no PSF matching is used, but we may want to tweak the PSF of our internal model, e.g. by oversampling Blend with a narrower PSF, which would require the user to remember that they now have to specify the observation PSF.

My proposal is minimal. Instead of https://github.com/fred3m/scarlet/blob/master/scarlet/observation.py#L129 (or the same location in LowResObservation), use if self.psfs is not None or self.psfs != scene.psfs

This will work as intended in all cases apart from this one: observation has a PSF, but scene doesn't, which should never occur.

Multi-component sources are simply components that are co-centered. We could reduce the complexity in Blend and Source considerably if instead of the multi-component implementation, we had a mechanism to tie the centroid updates for multiple sources/components together during the fitting procedure.

Side remark: we want to ability to put constraints on the Blend-level by relating sources to each other, e.g. by requiring that sed or morphology are the same across multiple exposures, or that the center moves on a particular track. We don't have that capability yet, but it would further extend the control over the Source behavior.

gcc (Ubuntu 6.2.0-5ubuntu12) 6.2.0 20161005
Python 3.5.2 |Continuum Analytics, Inc.| (default, Jul 2 2016, 17:53:06)

[esheldon@forest deblender] python setup.py install
Packages: ['deblender']
Packages: ['deblender']
/home/esheldon/miniconda3/lib/python3.5/site-packages/setuptools-27.2.0-py3.5.egg/setuptools/dist.py:340: UserWarning: The version specified ('0.0.8b960d2') is an invalid version, this may not work as expected with newer versions of setuptools, pip, and PyPI. Please see PEP 440 for more details.
running install
running bdist_egg
running egg_info
writing deblender.egg-info/PKG-INFO
writing requirements to deblender.egg-info/requires.txt
writing top-level names to deblender.egg-info/top_level.txt
writing dependency_links to deblender.egg-info/dependency_links.txt
reading manifest file 'deblender.egg-info/SOURCES.txt'
writing manifest file 'deblender.egg-info/SOURCES.txt'
installing library code to build/bdist.linux-x86_64/egg
running install_lib
running build_py
running build_ext
gcc -pthread -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esheldon/miniconda3/include/python3.5m -c /home/esheldon/data/tmp/tmp9putaogg.cpp -o home/esheldon/data/tmp/tmp9putaogg.o -std=c++14
cc1plus: warning: command line option ‘-Wstrict-prototypes’ is valid for C/ObjC but not for C++
gcc -pthread -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esheldon/miniconda3/include/python3.5m -c /home/esheldon/data/tmp/tmpdrzmkkoz.cpp -o home/esheldon/data/tmp/tmpdrzmkkoz.o -fvisibility=hidden
cc1plus: warning: command line option ‘-Wstrict-prototypes’ is valid for C/ObjC but not for C++
building 'deblender.proximal_utils' extension
gcc -pthread -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/esheldon/miniconda3/include/python3.5m -I/home/esheldon/.local/include/python3.5m -I/home/esheldon/miniconda3/include/python3.5m -c deblender/proximal_utils.cc -o build/temp.linux-x86_64-3.5/deblender/proximal_utils.o -DVERSION_INFO="0.0.8b960d2" -std=c++14 -fvisibility=hidden
cc1plus: warning: command line option ‘-Wstrict-prototypes’ is valid for C/ObjC but not for C++
deblender/proximal_utils.cc: In function ‘void prox_monotonic(pybind11::array_t<double>&, const double&, const std::vector<int>&, const std::vector<int>&, const double&)’:
deblender/proximal_utils.cc:16:14: error: ‘class pybind11::array_t<double>’ has no member named ‘mutable_unchecked’
   auto x = X.mutable_unchecked<2>();
              ^~~~~~~~~~~~~~~~~
deblender/proximal_utils.cc:16:35: error: expected primary-expression before ‘)’ token
   auto x = X.mutable_unchecked<2>();
                                   ^
error: command 'gcc' failed with exit status 1

The following code fails when a blend is created and fit:

constraints = (
    sc.SimpleConstraint() 
    & sc.MonotonicityConstraint(use_nearest=True)
    & sc.SymmetryConstraint()
)
sources = [scarlet.source.PointSource(
    (src['y'],src['x']),
    images,
    (15,15),
    constraints=constraints,
) for src in catalog]

because all of the sources are initialized the same ConstraintList object,
but not all of the objects will be the same size.

If we add a copy method to Constraint and ConstraintList, each source can have
a unique ConstraintList object that will be the correct size.

The are currently two conflicting ways to init a source: one that just makes a stamp of a given shape, and one that uses the data to define the source. The latter one is more powerful, and we know have a call_back mechanism to allows users to specify the init style.

This issue should get rid of the first/default source init, and remove the shape keyword as it'll be obsolete then.

The proposal is to call init_func before doing anything else, that is before the sizes are set and containers and constraint matrices are created. The callback can know the image data or whatever else might be needed, so it breaks the dependency between Source and Blend that gave rise to the delayed initialization of Source through Blend. It thus makes the entire thing both flexible and transparent to the user.

The thresholding of the PSF pixels to get a sparser representation is a bad idea. Instead sparisifying should be done for instance by a radial apodization function. But this needs to be left to the user to decide. I suggest that we drop this argument.

The inversion is unstable when a PSF is set. This can very likely be prevented by a ridge regression penalty, set at several times 1/cond(P).

with #78 a lot of functionality has been removed. scipy is only still used in two methods:

• operator.prox_strict_monotonic
• psf_match.fit_target_psf
  If possible, it would be nice to remove the dependency

This appears for me when I run the deblender notebook with psfs=psfs instead of psfs=None:

/Users/pmelchior/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/deblender-0.0.821b962-py2.7-macosx-10.6-x86_64.egg/deblender/nmf.pyc in deblend(img, peaks, constraints, weights, psf, max_iter, sky, l0_thresh, l1_thresh, e_rel, e_abs, monotonicUseNearest, traceback, translation_thresh, prox_A, prox_S, slack, update_order, steps_g, steps_g_update, truncate, fit_positions, txy_diff, max_shift, txy_thresh, txy_wait, txy_skip, Translation)
    352     S = init_S(N, M, K, img=_img, peaks=peaks)
    353     T = Translation(peaks, (N,M), B, P_, txy_diff, max_shift,
--> 354                     txy_thresh, fit_positions, txy_wait, txy_skip)
    355 
    356     # constraints on S: non-negativity or L0/L1 sparsity plus ...

/Users/pmelchior/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/deblender-0.0.821b962-py2.7-macosx-10.6-x86_64.egg/deblender/operators.pyc in __init__(self, *args, **kwargs)
    129         # Create the initial translations
    130         for k, (px, py) in enumerate(self.peaks):
--> 131             self.translate_psfs(k, update=True)
    132 
    133     def build_Tx(self, int_dx):

/Users/pmelchior/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/deblender-0.0.821b962-py2.7-macosx-10.6-x86_64.egg/deblender/operators.pyc in translate_psfs(self, k, ddx, ddy, update)
    212         """
    213         self.get_translation_ops(k, ddx, ddy, update)
--> 214         self.build_Gamma(k, update=update)
    215         return self.Gamma[k]
    216 

/Users/pmelchior/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/deblender-0.0.821b962-py2.7-macosx-10.6-x86_64.egg/deblender/operators.pyc in build_Gamma(self, k, Tx, Ty, update)
    201         else:
    202             Gamma_k = []
--> 203             for b in range(B):
    204                 g = Ty.dot(self.P[b].dot(Tx))
    205                 Gamma_k.append(g)

NameError: global name 'B' is not defined

So far, we've always assumed the A is normalized and the we get the flux of a component by summing of S. That's not necessarily true anymore. Problems arise e.g. here:
https://github.com/fred3m/scarlet/blob/master/scarlet/source.py#L62

We need to provide a per-component method, e.g. get_flux() that is aware of the normalization. This is a problem with the current scheme because the constraints are held by the components, but the normalization is held by the sources. Operationally both are constraints and should be treated equally, in which case both are held by the components. It's true that it's unlikely that someone would ever mix two components with different normalizations, but we've decided to retain that option.

So, this issue is for

• creating Component.get_flux()
• checking that we consistently call this method whenever the component flux is required.

This leaves one problem: what do we mean by flux? I argue it needs to be bolometric, i.e. flux = morph.sum() * sed.sum(). This form is -- yeah -- independent of normalization, which simplifies the problem considerably: components don't need to carry the normalization around at all.

Another issue will address the "normalization is a constraint problem".

Acceleration creates more worries than it helps. Plus the new S-type normalizations speed up convergence already. They will break backward compatibility somewhat in that the normalization isn't where users have come to expect it, but it is the right move.

Hi Fred and Peter,

I am trying to rerun the tests that I previously ran to see how the plots now look with the updates to proxmin. After I pulled, I tried executing the code in the following picture and got these errors:

Because the creation of the transformation operators (i.e. L matrices) is fairly costly, and we want to store them per component (i.e. for each small cutout, #14), we should create them once for a grid of image shapes. I suggest the following ranges:

• 19x19
• 29x29
• 39x39
• 59x59
• 79x79
• 99x99
• 149x149
• 199x199

Those matrices would be stored in a new directory scarlet/cache such that they can easily be read back: with scipy.sparse.save_npz and scipy.sparse.load_npz.

Hi,
I am currently using Galsim to create images of two galaxies and then using deblender to create a model of the two. I enter the true coordinates of the two objects when I create the sources, and although a model is created for both objects, sometimes the second one is empty (as shown in image below).

This is how I run the deblender. Note that img is the original image created by Galsim and coords are the coordinates of the two objects in the original. I am using bg_rms = 100.

def make_model(img,bg_rms,B,coords):
    #constraints on morphology:                                                   
    # "S": symmetry                                                               
    # "m": monotonicity (with neighbor pixel weighting)                           
    # "+": non-negativity                                                         
    constraints = {"S": None, "m": {'use_nearest': False}, "+": None}

    # initial size of the box around each object                                  
    # will be adjusted as needed                                                  
    shape = (B, 15, 15)
    sources = [scarlet.Source(coord, shape, constraints=constraints) for coord in\
 coords]
    blend = scarlet.Blend(sources, img, bg_rms=[bg_rms])
    # if you have per-pixel weights:                                              
    #weights = np.empty((B,Ny,Nx))                                                
    #blend = scarlet.Blend(sources, img, bg_rms=bg_rms, weights=weights)          

    # if you have per-band PSF kernel images:                                     
    # Note: These need to be difference kernels to a common minimum               
    #pdiff = [PSF[b] for b in range(B)]                                           
    #psf = scarlet.transformations.GammaOp(shape, pdiff)                          
    #blend = scarlet.Blend(sources, img, bg_rms=bg_rms, psf=psf)                  
    # run the fitter for 200 steps (or until convergence)                         

    blend.fit(200)#, e_rel=1e-1)                                                  
    # render the multi-band model: has same shape as img                          
    model = blend.get_model()
    mod2 = blend.get_model(m=1)
    return model,mod2

The second issue I am having is that sometimes the code crashes with a ValueError or numpy.linalg.linalg.LinAlgError (shown below). Thanks in advance for your help!

ValueError

/Users/lorena/anaconda3/lib/python3.6/site-packages/proxmin-0.4.3-py3.6.egg/proxmin/operators.py:37: RuntimeWarning: invalid value encountered in true_divide
/Users/lorena/anaconda3/lib/python3.6/site-packages/scarlet-0.1.d2e9b44-py3.6-macosx-10.7-x86_64.egg/scarlet/source.py:149: RuntimeWarning: invalid value encountered in less
  morph[morph<0] = 0
/Users/lorena/anaconda3/lib/python3.6/site-packages/scarlet-0.1.d2e9b44-py3.6-macosx-10.7-x86_64.egg/scarlet/source.py:161: RuntimeWarning: invalid value encountered in greater
  ypix, xpix = np.where(morph>blend._bg_rms[band]/2)
/Users/lorena/anaconda3/lib/python3.6/site-packages/scarlet-0.1.d2e9b44-py3.6-macosx-10.7-x86_64.egg/scarlet/source.py:163: RuntimeWarning: invalid value encountered in greater
  ypix, xpix = np.where(morph>0)
Traceback (most recent call last):
  File "galsim_deblend.py", line 194, in <module>
    model,mod2 = make_model(img,bg_rms,B,coords)
  File "galsim_deblend.py", line 80, in make_model
    blend = scarlet.Blend(sources, img, bg_rms=bg_rms)
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/scarlet-0.1.d2e9b44-py3.6-macosx-10.7-x86_64.egg/scarlet/blend.py", line 76, in __init__
    self.init_sources()
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/scarlet-0.1.d2e9b44-py3.6-macosx-10.7-x86_64.egg/scarlet/blend.py", line 255, in init_sources
    self.sources[m].init_source(self, self._img)
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/scarlet-0.1.d2e9b44-py3.6-macosx-10.7-x86_64.egg/scarlet/source.py", line 579, in init_source
    self.init_func(self, blend, img)
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/scarlet-0.1.d2e9b44-py3.6-macosx-10.7-x86_64.egg/scarlet/source.py", line 164, in init_templates
    Ny = np.max(ypix)-np.min(ypix)
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/numpy/core/fromnumeric.py", line 2272, in amax
    out=out, **kwargs)
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/numpy/core/_methods.py", line 26, in _amax
    return umr_maximum(a, axis, None, out, keepdims)
ValueError: zero-size array to reduction operation maximum which has no identity

LinAlgError

/Users/lorena/anaconda3/lib/python3.6/site-packages/proxmin-0.4.3-py3.6.egg/proxmin/utils.py:340: RuntimeWarning: divide by zero encountered in double_scalars
/Users/lorena/anaconda3/lib/python3.6/site-packages/proxmin-0.4.3-py3.6.egg/proxmin/utils.py:340: RuntimeWarning: invalid value encountered in multiply
/Users/lorena/anaconda3/lib/python3.6/site-packages/proxmin-0.4.3-py3.6.egg/proxmin/utils.py:310: RuntimeWarning: divide by zero encountered in double_scalars
/Users/lorena/anaconda3/lib/python3.6/site-packages/proxmin-0.4.3-py3.6.egg/proxmin/utils.py:310: RuntimeWarning: invalid value encountered in multiply
/Users/lorena/anaconda3/lib/python3.6/site-packages/proxmin-0.4.3-py3.6.egg/proxmin/utils.py:358: RuntimeWarning: divide by zero encountered in double_scalars
/Users/lorena/anaconda3/lib/python3.6/site-packages/proxmin-0.4.3-py3.6.egg/proxmin/utils.py:360: RuntimeWarning: divide by zero encountered in double_scalars
/Users/lorena/anaconda3/lib/python3.6/site-packages/proxmin-0.4.3-py3.6.egg/proxmin/utils.py:360: RuntimeWarning: invalid value encountered in true_divide
Traceback (most recent call last):
  File "galsim_deblend.py", line 194, in <module>
    model,mod2 = make_model(img,bg_rms,B,coords)
  File "galsim_deblend.py", line 92, in make_model
    blend.fit(200)#, e_rel=1e-1)
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/scarlet-0.1.d2e9b44-py3.6-macosx-10.7-x86_64.egg/scarlet/blend.py", line 229, in fit
    e_rel=self._e_rel, e_abs=self._e_abs, accelerated=accelerated, traceback=traceback)
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/proxmin-0.4.3-py3.6.egg/proxmin/algorithms.py", line 487, in bsdmm
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/scarlet-0.1.d2e9b44-py3.6-macosx-10.7-x86_64.egg/scarlet/blend.py", line 541, in _steps_f
    self._stepAS = [self._cbAS[block](block) for block in [0,1]]
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/scarlet-0.1.d2e9b44-py3.6-macosx-10.7-x86_64.egg/scarlet/blend.py", line 541, in <listcomp>
    self._stepAS = [self._cbAS[block](block) for block in [0,1]]
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/proxmin-0.4.3-py3.6.egg/proxmin/utils.py", line 222, in __call__
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/scarlet-0.1.d2e9b44-py3.6-macosx-10.7-x86_64.egg/scarlet/blend.py", line 486, in _one_over_lipschitz
    LA = np.real(np.linalg.eigvals(SSigma_1S).max())
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/numpy/linalg/linalg.py", line 889, in eigvals
    _assertFinite(a)
  File "/Users/lorena/anaconda3/lib/python3.6/site-packages/numpy/linalg/linalg.py", line 217, in _assertFinite
    raise LinAlgError("Array must not contain infs or NaNs")
numpy.linalg.linalg.LinAlgError: Array must not contain infs or NaNs

Currently the model of each component covers all pixels in the image. That's overkill for the evaluation of the likelihood gradients. Instead one could work with a subset of pixels, loosely based on previously detected footprints.

One can think of a few ways of doing that:

1. Each component exists only in its box, so that constraints only have to check for that box (we'd have to insist that the models are zero outside of that bounding box otherwise each component can do whatever it pleases outside of its own box). For the all procedures that require the actual model, e.g. likelihood gradients, the components need to be translated accordingly. That means that S is internally always represented by a vector (not a matrix), which we already partially implement in the dot products of block-diagonal constraint matrices, and the generation of the model is more involved. The downside is that we need to build the constraint matrices separately for each component (or work with fixed box size increments to reduce that number). It seems to me that this will increase initialization time but reduce runtime.
2. Use a sparse matrix form for S that allows us to quickly update a set of predefined elements (for every component: those inside of the box, the rest is zero) as long as we don't change the structure of the matrix. The models would automatically have the right format, and the changes to the code would be minimal. There is probably a slightly larger runtime cost. This approach seems fairly straightforward to me, but I can't seem to find a way to do this within scipy.sparse.

In both cases, we'd automatically prevent the degeneracy between two remote objects having the same color.

Here is the output from the logger showing that one component is shifted by a very large amount:

DEBUG:scarlet.blend:shifting component (0, 0) by (-0.048/0.218) to (32.887/32.236) in it 190
DEBUG:scarlet.blend:shifting component (1, 0) by (-89806566564.740/-75030309405.693) to (-89806566533.021/-75030309379.929) in it 190
DEBUG:scarlet.blend:no flux in component (1, 0), skipping recentering in it 200
DEBUG:scarlet.blend:shifting component (0, 0) by (-0.024/0.018) to (32.863/32.255) in it 200```

The fitter occasionally stops after 10 steps. This is not a convergence issue, it's very likely related to the refine_skip parameter and the RestartException being thrown but not properly treated.

The resizing operation is expensive, so we should attempt to minimize them. One way is to be clever about predicting the final size of the source boxes instead of expanding the boxes several times.

@fred3m suggested to extrapolate from the peak over the edge until one reaches the bg_rms. That's a step forward for large objects. Problems might be that the slope in the tails of galaxies is shallower than linear, and that one has to incorporate variable noise levels across the bands. Nonetheless, this should help.

We need to allow users to construct their own PSF difference kernels to match PSF variations between bands. Galsim has all we need to doing the heavy lifting.

The function deblender.proximal.strict_monotonicity takes peaks and components as separate arguments. This should be modified to accept a deblender.nmf.PeakCatalog instead to ensure that all of the components that require monotonicity have the proximal operator applied.

Going from gamma = None to gamma = GammaOp(psf) within a session means that the cache of the old gamma is used for the new. Since the cache is global, we need to

1. destroy the cache when GammaOp goes out of scope
2. make cache property of GammaOp.

The latter option is much safer and more transparent. We don't have to do this for all transformation operators because only GammaOp can be created by the user in a regular fashion.

The explicit estimation of the Lipschitz constants is expensive and doesn't fully consider all elements of the observations (convolutions, spectral resolution etc). Instead, try to determine them from the change in gradients.

This is for the new fullscene code:

The convergence test we perform is purely relative, i.e. the L2 norm of changes morphology/SED vectors need to be less than e.g. 1e-2 of the L2 norm of the vector itself, but that's hard to do if the vector is small to begin with.

The complication is that the normalization option tells us which vector is normalized to unity, making the other one small. But normalization is under user control, so we don't know that within Blend.

The second complication is that an absolute tolerance would have to be set in relation to the noise levels, but this is very hard to propagate from the observation to the model.

It would nonetheless be useful to check if there are meaningful changes from further iterations, so that we don't waste iterations in incremental improvements.

One option would be to check Blend.mse for changes and terminate if it doesn't improve by some amount, to be set by the user e.g. Blend.fit(100, e_rel=1e-2, e_mse=1e-3).

When installing scalet with python2.7, the installation crashes as shown below.
This is fixed easily by either replacing the <> by "" in the include lines, or by installing with python3.

joseph:scarlet remy$ sudo python2.7 setup.py install
/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/setuptools/dist.py:406: UserWarning: The version specified (u'0.3.b34b32b') is an invalid version, this may not work as expected with newer versions of setuptools, pip, and PyPI. Please see PEP 440 for more details.
  "details." % self.metadata.version
running install
running bdist_egg
running egg_info
writing requirements to scarlet.egg-info/requires.txt
writing scarlet.egg-info/PKG-INFO
writing top-level names to scarlet.egg-info/top_level.txt
writing dependency_links to scarlet.egg-info/dependency_links.txt
reading manifest file 'scarlet.egg-info/SOURCES.txt'
writing manifest file 'scarlet.egg-info/SOURCES.txt'
installing library code to build/bdist.macosx-10.13-x86_64/egg
running install_lib
running build_py
running build_ext
/usr/bin/clang -fno-strict-aliasing -fno-common -dynamic -pipe -Os -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -I/opt/local/Library/Frameworks/Python.framework/Versions/2.7/include/python2.7 -c /tmp/tmpTW9gRr.cpp -o tmp/tmpTW9gRr.o -std=c++14
/usr/bin/clang -fno-strict-aliasing -fno-common -dynamic -pipe -Os -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -I/opt/local/Library/Frameworks/Python.framework/Versions/2.7/include/python2.7 -c /tmp/tmpkRvneG.cpp -o tmp/tmpkRvneG.o -fvisibility=hidden
building 'scarlet.operators_pybind11' extension
/usr/bin/clang -fno-strict-aliasing -fno-common -dynamic -pipe -Os -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -Iinclude -Iinclude -I/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/peigen-0.0.9-py2.7.egg/peigen/include -I/opt/local/Library/Frameworks/Python.framework/Versions/2.7/include/python2.7 -c scarlet/operators_pybind11.cc -o build/temp.macosx-10.13-x86_64-2.7/scarlet/operators_pybind11.o -stdlib=libc++ -mmacosx-version-min=10.7 -DVERSION_INFO="0.3.b34b32b" -std=c++14 -fvisibility=hidden
scarlet/operators_pybind11.cc:1:10: error: 'pybind11/pybind11.h' file not found with <angled> include; use "quotes" instead
#include <pybind11/pybind11.h>
         ^~~~~~~~~~~~~~~~~~~~~
         "pybind11/pybind11.h"
scarlet/operators_pybind11.cc:2:10: error: 'pybind11/numpy.h' file not found with <angled> include; use "quotes" instead
#include <pybind11/numpy.h>
         ^~~~~~~~~~~~~~~~~~
         "pybind11/numpy.h"
scarlet/operators_pybind11.cc:3:10: error: 'pybind11/stl.h' file not found with <angled> include; use "quotes" instead
#include <pybind11/stl.h>
         ^~~~~~~~~~~~~~~~
         "pybind11/stl.h"
scarlet/operators_pybind11.cc:4:10: error: 'pybind11/eigen.h' file not found with <angled> include; use "quotes" instead
#include <pybind11/eigen.h>
         ^~~~~~~~~~~~~~~~~~
         "pybind11/eigen.h"
scarlet/operators_pybind11.cc:72:1: warning: 'pybind11_init' is deprecated: PYBIND11_PLUGIN is deprecated, use PYBIND11_MODULE
      [-Wdeprecated-declarations]
PYBIND11_PLUGIN(operators_pybind11)
^
scarlet/pybind11/detail/common.h:261:20: note: expanded from macro 'PYBIND11_PLUGIN'
            return pybind11_init();                                            \
                   ^
scarlet/operators_pybind11.cc:72:1: note: 'pybind11_init' has been explicitly marked deprecated here
scarlet/pybind11/detail/common.h:256:5: note: expanded from macro 'PYBIND11_PLUGIN'
    PYBIND11_DEPRECATED("PYBIND11_PLUGIN is deprecated, use PYBIND11_MODULE")  \
    ^
scarlet/pybind11/detail/common.h:90:41: note: expanded from macro 'PYBIND11_DEPRECATED'
#  define PYBIND11_DEPRECATED(reason) [[deprecated(reason)]]
                                        ^
1 warning and 4 errors generated.
error: command '/usr/bin/clang' failed with exit status 1

Currently, the constraints on S are evaluated in the "centered" frame, and then the components are shifted to the peak location. This means the model can assume arbitrary values in the regions that get shifted out of the observed postage stamp. Nominally, the monotonicity constraint should prevent greater damage, but it doesn't seem to do it.
Alternatively, I suggest to add another strict penalty to constrain all pixels that are shifted out of the observed image to be zero: prox(x) = x if in frame; 0 if not.

Continuation of this issue pmelchior/proxmin#14

After a successful compilation with gcc 4.8.22 with std=c++11, there's an import error:

>>> import deblender
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/u/mjerdee/.local/lib/python2.7/site-packages/deblender-0.0.1f00041-py2.7-linux-x86_64.egg/deblender/__init__.py", line 2, in <module>
    from . import nmf
  File "/u/mjerdee/.local/lib/python2.7/site-packages/deblender-0.0.1f00041-py2.7-linux-x86_64.egg/deblender/nmf.py", line 13, in <module>
    from .proximal import build_prox_monotonic
  File "/u/mjerdee/.local/lib/python2.7/site-packages/deblender-0.0.1f00041-py2.7-linux-x86_64.egg/deblender/proximal.py", line 8, in <module>
    from . import proximal_utils
ImportError: /u/mjerdee/.local/lib/python2.7/site-packages/deblender-0.0.1f00041-py2.7-linux-x86_64.egg/deblender/proximal_utils.so: undefined symbol: _ZTVN10__cxxabiv120__function_type_infoE

please advise.

this can be handled by a clever treatment of & in Constraint itself

a new segment in my continuing mission to align code structures with natural language

Commit 1d87a84 updated the code to use weights for A and S that are normalized differently. This change needs to be ported to the scarlet branch.