pmelchior / scarlet2

Scarlet, all new and shiny
MIT License
15 stars 4 forks source link

Deimos #65

Closed pmelchior closed 4 months ago

pmelchior commented 5 months ago

This PR implements the moment deconvolution method (closes #48), so that an initial Gaussian from adaptive_gaussian_morph is computed for the deconvolved model frame rather then the convolved observed frame.

It also updates the spectrum using the zeroth moment, which is the integrated flux. This is somewhat dangerous in the presence of nearby objects (in that case pixel_spectrum is the safest choice), but it has lower variance (more pixel contribute).

Because of the deconvolution, the moments might be ill-defined, which will throw a ValueError, so the recommended pattern is:

try:
    spectrum, morph = init.adaptive_gaussian_morph(obs, center, min_corr=0.99)
except ValueError:
    spectrum = init.pixel_spectrum(obs, center)
    morph = init.compact_morphology()
SampsonML commented 5 months ago

When running on this branch using the above supplied initialization routine I run into this error

File [~/Research/Melchior/2024_scarlet_version/scarlet2/scarlet2/spectrum.py:18](https://file+.vscode-resource.vscode-cdn.net/Users/mattsampson/Research/Melchior/2024_scarlet_version/scarlet2/~/Research/Melchior/2024_scarlet_version/scarlet2/scarlet2/spectrum.py:18), in ArraySpectrum.__init__(self, data)
     16 def __init__(self, data):
     17     self.data = data
---> 18     self.bbox = Box(self.data.shape)

AttributeError: 'Parameter' object has no attribute 'shape'

The code I am calling is this

from numpyro.distributions import constraints
from scarlet2 import relative_step
from functools import partial 

spec_step = partial(relative_step, factor=5e-3)
morph_step = partial(relative_step, factor=5e-3)

with Scene(model_frame) as scene:
    for center in centers:

        try:
            spectrum, morph = init.adaptive_gaussian_morph(obs, center, min_corr=0.99)
        except ValueError:
            spectrum = init.pixel_spectrum(obs, center)
            morph = init.compact_morphology()

        print(f"morph shape is {np.shape(morph)}")
        print(f"spectrum shape is {np.shape(spectrum)}")

        Source(
            center,
            ArraySpectrum(
                Parameter(spectrum, constraint=constraints.positive, stepsize=spec_step)
            ),
            ArrayMorphology(Parameter(morph, prior=prior, stepsize=morph_step)),
        )
b-remy commented 4 months ago

Yes since #56, Parameters are to be set after the scene initialization. This should work instead:

from numpyro.distributions import constraints
from scarlet2 import relative_step
from functools import partial 

spec_step = partial(relative_step, factor=5e-3)
morph_step = partial(relative_step, factor=5e-3)

with Scene(model_frame) as scene:
    for center in centers:

        try:
            spectrum, morph = init.adaptive_gaussian_morph(obs, center, min_corr=0.99)
        except ValueError:
            spectrum = init.pixel_spectrum(obs, center)
            morph = init.compact_morphology()

        print(f"morph shape is {np.shape(morph)}")
        print(f"spectrum shape is {np.shape(spectrum)}")

        Source(
            center,
            ArraySpectrum(spectrum),
            ArrayMorphology(morph),
        )

parameters = scene.make_parameters()
for i in range(len(scene.sources)):
    parameters += Parameter(scene.sources[i].spectrum.data, name=f"spectrum.{i}", constraint=constraints.positive, stepsize=spec_step)
    parameters += Parameter(scene.sources[i].morphology.data, name=f"morph.{i}", prior=prior, stepsize=morph_step)
SampsonML commented 4 months ago

This works for me now.