Herculens / herculens_workspace

Example notebooks and other resources for Herculens
MIT License
5 stars 1 forks source link

Fit using 2 analytic sources or lenses possible? #5

Closed PrayaagKatta closed 6 months ago

PrayaagKatta commented 7 months ago

I wanted to know if adding 2 analytic sources or analytic profiles of light from 2 lenses is possible, and if so how to do it. I dont want to switch over to a pixellated model as of now, but do let me know if that's the only option! Thank you.

aymgal commented 7 months ago

Hi @PrayaagKatta , yes you can simply add a second profile in the list, e.g., LightModel('SERSIC_ELLIPSE', 'SERSIC_ELLIPSE'), and make sure to have the associated kwargs_light to be consistent. Let me know if you still have issues!

PrayaagKatta commented 7 months ago

How do I change the ProbModel class to account for this? I tried doing this, will it work? TIA!

def model(self):
        # Parameters of the source
        prior_sources = []
        for i in range(len(previous_kwargs_result['kwargs_source'])):
            prior_sources.append({
                'amp': numpyro.sample(f'source_amp_{i}', dist.LogNormal(np.log(previous_kwargs_result['kwargs_source'][i]['amp']), 1)),
                'R_sersic': numpyro.sample(f'source_R_sersic_{i}', dist.TruncatedNormal(previous_kwargs_result['kwargs_source'][i]['R_sersic'], 0.5, low=0.05)),
                'n_sersic': numpyro.sample(f'source_n_{i}', dist.Uniform(1., 4.)),
                'e1': numpyro.sample(f'source_e1_{i}', dist.TruncatedNormal(previous_kwargs_result['kwargs_source'][i]['e1'], 0.5, low=-0.5, high=0.5)),
                'e2': numpyro.sample(f'source_e2_{i}', dist.TruncatedNormal(previous_kwargs_result['kwargs_source'][i]['e2'], 0.5, low=-0.5, high=0.5)),
                'center_x': numpyro.sample(f'source_center_x_{i}', dist.Normal(previous_kwargs_result['kwargs_source'][i]['center_x'], 0.1)),
                'center_y': numpyro.sample(f'source_center_y_{i}', dist.Normal(previous_kwargs_result['kwargs_source'][i]['center_y'], 0.1))
            })
...
def params2kwargs(self, params):
    for i in range(len(previous_kwargs_result['kwargs_source'])):
            source_params_dict = {
                'amp': params[f'source_amp_{i}'],
                'R_sersic': params[f'source_R_sersic_{i}'],
                'n_sersic': params[f'source_n_{i}'],
                'e1': params[f'source_e1_{i}'],
                'e2': params[f'source_e2_{i}'],
                'center_x': params[f'source_center_x_{i}'],
                'center_y': params[f'source_center_y_{i}']
            }

            kw['kwargs_source'].append(source_params_dict)
aymgal commented 6 months ago

Your update to the ProbModel looks good to me! The important point is that the model() and params2kwargs() methods are consistent with each other, and consistent with the lists of profiles you choose for your different model components.

If you source model is LightModel('SERSIC_ELLIPSE', 'SERSIC_ELLIPSE'), then it looks like your ProbModel is correct. Did that work?

The other important thing when using JAX (and thus numpyro and Herculens) is that any for loop should not depend on the parameters values, so your code should be fine.

aymgal commented 6 months ago

@PrayaagKatta Do you still have issues with the model? If not, you can close the issue. Otherwise, please let me know!

PrayaagKatta commented 6 months ago

Thanks for your help, I closed the issue! It is running, but the final model is not close to the mock data (Pic attached below). I'm assuming that I just have to play around with the width of the distributions it is sampling from? Or do I use another solver method for the larger number of parameters? (Which would that be?) Sorry for the additional questions, let me know if you need any info from my side. Thanks again for your help! Model_PJ_mass-9_1

aymgal commented 6 months ago

Thanks @PrayaagKatta . Ah I see. It a rather difficult task to fit for two Sérsic blobs in source plane, as there are probably lots of degeneracies between the two, and with the mass model. Do you attempt to fit the whole parameter at once, or do you first obtain a initial mass model using a single Sérsic in the source, and only then increasing the source complexity?

For continuing the discussion if needed, what about opening a discussion on the Herculens repo? This is better suited for this kind of matter, and might be useful to other people. Thanks!