Slam script incorrect ResultInterferometer and ResultImaging adapt_image calls about pyautolens HOT 58 CLOSED

I will raise an issue on PyAutoFit to have more detailed logging of what the parallelization is doing.

Thanks!

In the sh file for starting the run I write:

#!/bin/bash -l
#SBATCH --nodes=1
#SBATCH --partition=riechers-exclusive
#SBATCH --mail-type=ALL
#SBATCH --no-requeue
#SBATCH --ntasks-per-node=16
#SBATCH --cpus-per-task=1
#SBATCH --time=24-12:00:00
#SBATCH --mem=128GB

source /projects/ag-riechers/users/aaron/ALMACE/PyAutoLensNN/activate.sh

export CPUS_PER_TASK=1
export OPENBLAS_NUM_THREADS=$CPUS_PER_TASK
export MKL_NUM_THREADS=$CPUS_PER_TASK
export OMP_NUM_THREADS=$CPUS_PER_TASK
export VECLIB_MAXIMUM_THREADS=$CPUS_PER_TASK
export NUMEXPR_NUM_THREADS=$CPUS_PER_TASK

python3 /projects/ag-riechers/users/aaron/ALMACE/PyAutoLensNN/source_pixelized_2LensesNew.py

So that should be fine.

from pyautolens.

I dont think we ever implemented parallel initializer, got an issue up now:

rhayes777/PyAutoFit#955

from pyautolens.

Thanks again for the quick reply! It's so helpful!

adapt_image certainly makes more sense to me now, thanks!

The adapt_image probably has a bit of correlated noise and other features, but it seems to be good enough to do it this way.

Eventually, when I get to my lower S/N data this might come up again, but I hope for now this will be the case throughout.

This sounds correct. image_mesh size certainly has nothing to do with interferometer, its just how we distribute our source pixels in the image plane.

This is great! Because then we can just push down shape_native down to the FOV of the antenna at that frequency for a pixel_scale with 5 pixels/beam, use the smallest circular mask as possible that captures the emission and then adapt the image_mesh for our science and run time.

You should probably try a smaller mask for the 0.1686" lens, feels like you're making the calculation slower than necessary.

Oops, my mistake this is just the other lens galaxy in the same image. There are two with different einstein radii in the same image with some offset. The other einstein rings for the other data sets are all either smaller or minimally larger.

Yep, this sounds realistic. We have done a lot of careful crafting of our settings in the past, and this took us from months to days... so its defo doable! Apologies its such a nightmare though, and the docs are stlil a bit lacking.

Thank you so much for your support! This is so much better than all the other lensing codes I tried before. Your code has already good documentation and I really appreciate your support! Doing this on my own would take many more months without your input.

from pyautolens.

Thanks for the issue :).

I did a lot of updates to adapt_images in the past week, including a new release, and certain example scripts may not have been properly updated!

It sounds like you landed on the fix yourself, but post here if you have any other issues or the search which started appears to be stuck.

from pyautolens.

Thank you for this fantastic tool ;)

The interferometer run eventually failed after 9 bounds, possibly due to an out-of-memory issue leading to a core dump, and no further report in its log file. Since I am working on a hpc, I cannot track exactly if this was the case.

Re-running the slam script should've started the process from where it left off, it did but immediately gave the following error:

This might be a consequence of some import issues from the previous run but is a bit too deep for me to understand where the problem comes from.

I renamed the output folder defined in the source_pix.py equivalent (I use a slightly modified version) in the output folder in my file system (not in the script) and am running the pipeline again. It appears to be running the search using the results of the parametric run.

from pyautolens.

Yeah Im unsure on the error above, but it seems reasonable that the results folder had a few incorrect files output there and rerunning it from scrtch (by renaming a few folders) would sort it.

If errors persist I can do an end-to-end interferometer run myself and check that it all works smoothly.

from pyautolens.

After a couple of weeks, the first of the two pixelation runs was finally completed. However, this time it resulted in the same error as before:

Which means it was not an incorrect file issue earlier in the post. This also means it should be a general issue, so I would appreciate it greatly if you could run the end-to-end interferometer run for yourself. I have not updated the code
since three weeks ago, but I couldn't find anything in the recent commits where I could see this issue fixed yet.
To me, it appears that some kind of plotting function cannot read or find the pixel_scales from the real_space_mask in the dataset, since we don't specify pixel_scales directly in the "dataset" object.

from pyautolens.

I have pushed a fix to the main branch of PyAutoArray.

However, you may wish to edit your source code with the fix manually, just to make sure the code does not lose backwards compatibility with your old results (I don't think it would, but better safe than sorry).

The fix requires you to go to this script in your virtual environment:

https://github.com/Jammy2211/PyAutoArray/blob/main/autoarray/inversion/plot/inversion_plotters.py

Go to line 282:

            self.mat_plot_2d.plot_array(
                array=self.inversion.data,
                visuals_2d=self.get_visuals_2d_for_data(),
                auto_labels=AutoLabels(title=f" Data"),
            )

And replace it with the following code:

        try:
            self.mat_plot_2d.plot_array(
                array=self.inversion.data,
                visuals_2d=self.get_visuals_2d_for_data(),
                auto_labels=AutoLabels(title=f" Data"),
            )
        except AttributeError:
            pass

from pyautolens.

Thank you for the fast fix! Works like a charm now.

from pyautolens.

Hello again,

I have noticed that actually my pixel run started from scratch instead of continuing/finishing the run and producing all the required files, i.e. ".completed", "model.results", "search.summary", etc. So I think that, in the step of updating the final results, where it produced the error last time, it automatically removed the checkpoint.hdf5 file which enabled it to start again from there. Now it always wants to start from the beginning, even though the sampling is supposed to be finished. The other files, see screenshots, are still there, and also "samples_info.json" shows that it has all 26200 samples.

Is there a way to either: create the checkpoint hdf5 file from what I have from the json and csv files?
Or, a way to have the search somehow understand that it is finished already and only has to create the final files from what is there?

The files I have of the results:

from pyautolens.

Can you try putting a .completed file in the unique identifier folder and running again?

Autolens checks if this file exists, and if it does it should bypass the search and load the result for doing stuff like visualization from samples.csv.

from pyautolens.

That's a great idea! It works. Out of curiosity, should the second analysis for the second run in source_pix.py have a positions_likelihood from the parametric run as the first pixel run did? Speaking about interferometer data, of course.

from pyautolens.

source_pix[2] does not need a positions_likelihood, because the mass model is fixed during this search (and we therefore know the positions threshold will always be met.

You are defo gonna wanna make sure your mass[1] search has a positions_likelihood though!

from pyautolens.

Perfect, thank you! I did it correctly then. Fingers crossed that there are no further issues in this pixel run.
It failed the very same moment I wrote the above:

Afterward, I set raise_inversion_positions_likelihood_exception=False explicitly in the analysisInterferometer object used in the second search.

from pyautolens.

Strange, are you able to input a positions_likelihood into source_pix[2] (something like analysis=analysis_1.positions_likelihood.

No idea why that exception is being raised, I'll note its a bug.

Alternatively, you can prob also disable it by changing the following setting in general.yaml:

test:
  disable_positions_lh_inversion_check: true

from pyautolens.

Two ways worked for me. One is:

    analysis = al.AnalysisInterferometer(
        dataset=analysis.dataset,
        adapt_images=result_1.adapt_images_from(use_model_images=True),
        **raise_inversion_positions_likelihood_exception=False,**
        settings_inversion=al.SettingsInversion(
            image_mesh_min_mesh_pixels_per_pixel=3,
            image_mesh_min_mesh_number=5,
            image_mesh_adapt_background_percent_threshold=0.1,
            image_mesh_adapt_background_percent_check=0.8, use_linear_operators=True,
        ),
    )

And alternatively setting:
analysis_1 = al.AnalysisInterferometer( dataset=analysis.dataset, positions_likelihood=analysis.positions_likelihood, ...
I did not try changing the setting in general.yaml since it is working atm.

While it is running now, it is only using a single core and only spawns one thread to generate the initial samples of the source_pix[2] and that is taking forever, i.e. more than 14h (unknown how long it actually needs to take). It usually takes only ~2h to do the first 300 samples with nautilus. Is there a way of allowing multithreading or some sort of parallelization in the autofit.non_linear.initializerstep?

2024-03-01 11:25:51,288 - autofit.non_linear.initializer - INFO - Generating initial samples of model, which are subject to prior limits and other constraints.

I already have set the number_of_cores setting in non_linear.yaml to 16 (max. I can use). The same is defined in settings_search.

Edit: Nvm, now it did spawn the other threads. However, it is only using a single core. Before, I've only seen this when the multiprocessing didn't work. I tested DynestyStatic on my laptop and the multiprocessing works fine (for "Imaging" data), all 4 cores I set up were being used.

from pyautolens.

Is there a way of allowing multithreading or some sort of parallelization in the autofit.non_linear.initializerstep?

Not currently, I think, I can raise an issue on PyAutoFit for us to implement this (its not often it really matters but clearly for your usecase it does lol).

I'm unsure why the threading isn't quite doing what you'd expect for source_pix[2]. Sounds like you're just about ok now, but I will keep an eye on it when I do my next round of profiling autolens.

from pyautolens.

Edit 2: This might be related to the code running on an hpc. I did the same test that I did on the laptop on an hpc and still had the same issue. On the hpc, changing number_of_cores to =1 in nest.yaml fixed the issue for the no_lens_light.py slam pipeline, which I used as a test.

It's a bit strange that you can set number_of_cores=16 in the nest.yaml for nautilus and it works fine, while for DynestyStatic I need to set it to =1. I don't get this behavior on my laptop for some reason, where I had number_of_cores=4 for DynestyStatic in nest.yaml and it worked. Also, on hpc it now spawns 52 threads for some reason and is still only using a single core for the interferometer dataset.

Not currently, I think, I can raise an issue on PyAutoFit for us to implement this (its not often it really matters but clearly for your usecase it does lol).

I think that would be great for the future, for using ALMA datasets with millions of visibilities. But compared to total computation time it's probably only a minor time saver, great for testing though.

from pyautolens.

Can I just check have you set up all this stuff correctly:

              """
                    The non-linear search is using multiprocessing (number_of_cores>1). 

                    However, the following environment variables have not been set to 1:

                    OPENBLAS_NUM_THREADS
                    MKL_NUM_THREADS
                    OMP_NUM_THREADS
                    VECLIB_MAXIMUM_THREADS
                    NUMEXPR_NUM_THREADS

                    This can lead to performance issues, because both the non-linear search and libraries that may be
                    used in your `log_likelihood_function` evaluation (e.g. NumPy, SciPy, scikit-learn) may attempt to
                    parallelize over all cores available.

                    This will lead to slow-down, due to overallocation of tasks over the CPUs.

                    To mitigate this, set the environment variables to 1 via the following command on your
                    bash terminal / command line:

                    export OPENBLAS_NUM_THREADS=1
                    export MKL_NUM_THREADS=1
                    export OMP_NUM_THREADS=1
                    export VECLIB_MAXIMUM_THREADS=1
                    export NUMEXPR_NUM_THREADS=1

                    This means only the non-linear search is parallelized over multiple cores.

                    If you "know what you are doing" and do not want these environment variables to be set to one, you 
                    can disable this warning by changing the following entry in the config files:

                    `config -> general.yaml -> parallel: -> warn_environment_variable=False`
                    """

from pyautolens.

I will raise an issue on PyAutoFit to have more detailed logging of what the parallelization is doing.

I think part of the problem is that the code is not displaying its set up to you and you're having to do annoying things to check if number_of_cores is working.

from pyautolens.

Small update:
The initialiser only used one core and took ~24h to do that iteration. Thankfully, with the current settings, the actual sampling after this uses the 16 cores as specified in the search settings in the script and 1000 iterations take 14h. It would be great if it would be possible to have the autofit.non_linear.initializer step use multiprocessing. Or should this be the case already (since that step does spawn ~32 unused threads)?

from pyautolens.

I am still running the second to last step in the slam pixelation pipeline and I was wondering if there are ways to speed up the process that are not obvious, i.e. Averaging the visibilities more, using fewer Hilbert image mesh pixels (currently using default 1000pix).

See the search.summary file:

Is there a way to speed up the "Log Likelihood Function Evaluation Time"? This is using DynestyStatic, where I get pretty long evaluation times. I know that Time per sample is ~2min because of the large number of visibilities and it seems to me to only be optimizable by reducing the number of visibilities (by averaging). Still, I would like to avoid doing that for now.

I will put this in multiple questions:
Should the value for nlive (DynestyStatic) be customized in the slam pipeline?
I.e. be tuned to be nlive=integer*ncores for hpcs (and sufficiently high number of cores)?
Currently, it is set to 100.

Should the Hilbert image mesh pixels be changed from the default value in the second pixelization search?

Am I missing some setting in non_linear.yaml to speed up the evaluation time?

from pyautolens.

How many visibilities are you modeling?

Can you send the full information of the run_time_dict?

With a time per sample of 2 mins there probably isn't much you can do to speed things up by changing Dynesty.

We are working on a faster visibilities fitting method but... its not an easy problem lol.

from pyautolens.

This dataset has ~5mil visibilites.

I created a separate run of the slam pipeline and inserted the run_time code before searching the second time in the pixelization run:

run_time_dict, info_dict = analysis_1.profile_log_likelihood_function(
        instance=model_2.random_instance()
)
print(run_time_dict)
search_2=af.DynestyStatic(...)

So this uses the model that is also input into the second pixelation run.
Right now I get the error:

I'll look into why this happens.

So falling back to just the functions, as defined in the interferometer/data_preparation/run_times.py script, everything else is the same:

Transformer Time = 0.8292598724365234

15158617.926531568
Fit Time = 1.4524788856506348

We are working on a faster visibilities fitting method but... its not an easy problem lol.

I fully understand, and I am already very grateful for all the help you already provided. I was just wondering, whether my settings were wrong, as all previous sampling was done with Nautilus and it took me some time to get everything running smoothly.

(Sorry for the longer wait. I first did this just using the run time and fit time from the interferometer data preparation script by mistake and now there is also another error.)

from pyautolens.

I'll look into why this happens.

If you input an inaccurate lens model into a fit, it can often lead to a ilnear algebra error.

This wouldn't happen if you input an "accurate" lens model.

The input thing at the start obviously relies on the random model you put in working... I forgot this could happen when I wrote the docs.

So falling back to just the functions, as defined in the interferometer/data_preparation/run_times.py script, everything else is the same:

Transformer Time = 0.8292598724365234

15158617.926531568
Fit Time = 1.4524788856506348

I'm a bit confused, why is this giving a fit time of 1.45 if your actual runs are > 2 mins? Are you inputting the same data?

Can you send the search.summary files of source_lp[1], source_pix[1] and source_pix[2]?

from pyautolens.

The input thing at the start obviously relies on the random model you put in working... I forgot this could happen when I wrote the docs.

I just realised this today and was able to run the script to get the run_time_dict. After running it a couple times and it stopping at the inversion check check_mesh_pixels_per_image_pixels in image_mesh/abstract and giving the InversionException error, I removed the

settings_inversion = al.SettingsInversion(
	image_mesh_min_mesh_pixels_per_pixel=3,
    image_mesh_min_mesh_number=5,
    image_mesh_adapt_background_percent_threshold=0.1,
    image_mesh_adapt_background_percent_check=0.8,)

part of the settings for the analysis object that is input in the run_time calculation. And since it doesn't check now, I got the run_time_dict:

{'fit_time': 512.0117816925049, 'image_2d_from_0': 0.009083032608032227, 'visibilities_from_0': 0.03899717330932617, 'relocated_grid_from_0': 0.09293508529663086, 'relocated_mesh_grid_from_0': 0.26378417378417015075684, 'mesh_grid_from_0': 2.6226043701171875e-05, 'image_2d_from_1': 0.008669853210449219, 'visibilities_from_1': 0.039293527603149414, 'regularization_matrix_0': 0.1901397705078125, 'precondreconditioner_matrix_0': 1.760345220565796, 'preconditioner_matrix_inverse_0': 0.08354568481445312, 'reconstruction_0': 477.9175887107849, 'mapped_reconstructed_image_dict_0': 0.17968225479125977, 'mapped_reconstructed_data_dict_0': 0.855323314666748, 'mapped_reconstructed_data_0': 0.04035782814025879, 'reconstruction_reduced_0': 6.67572021484375e-06, 'regularization_matrix_reduced_0': 3.0994415283203125e-06, 'regularization_term_0': 0.0006058216094970703, 'log_det_curvature_reg_matrix_term_0': 0.014694690704345703, 'log_det_regularization_matrix_term_0': 0.013769149780273438, 'figure_ofure_of_merit_0': 30.28476905822754}
But I am leaving the above setting_inversion for my normal slam run.

I'm a bit confused, why is this giving a fit time of 1.45 if your actual runs are > 2 mins? Are you inputting the same data?

I was inputting the same data, however I used

print_transformer_time_from(
    dataset=dataset, transformer_class=al.TransformerNUFFT, repeats=1
)
print_fit_time_from(
    dataset=dataset,
    transformer_class=al.TransformerNUFFT,
    use_linear_operators=True,
)

with the tracer using the output from source_pix[1] but without the af.Model wrappers.

Can you send the search.summary files of source_lp[1], source_pix[1] and source_pix[2]?

source_lp[1] search.summary :

Total Samples = 42300
Total Accepted Samples = 42300
Acceptance Ratio = 1.0
Time To Run = 5:43:20.742212
Time Per Sample (seconds) = 0.48701518232377145
Log Likelihood Function Evaluation Time (seconds) = 2.7229983806610107
Expected Time To Run (seconds) = 1 day, 7:59:42.831502
Speed Up Factor (e.g. due to parallelization) = 5.591198138155102

source_pix[1] search.summary :
unfortunately was somehow removed when the save_state file was removed, since the the run died during plotting after finishing.
But here is the run_time_dict for it, should that help:
{'fit_time': 192.92181277275085, 'image_2d_from_0': 0.00858449935913086, 'visibilities_from_0': 0.038361310958862305, 'relocated_grid_from_0': 0.18001127243041992, 'relocated_mesh_grid_from_0': 0.2624642848968506, 'mesh_grid_from_0': 2.5033950805664062e-05, 'image_2d_from_1': 0.008388042449951172, 'visibilities_from_1': 0.038111209869384766, 'regularization_matrix_0': 0.18108797073364258, 'preconditioner_matrix_0': 1.4672739505767822, 'preconditioner_matrix_inverse_0': 0.06569290161132812, 'reconstruction_0': 188.47503900527954, 'mapped_reconstructed_image_dict_0': 0.16329526901245117, 'mapped_reconstructed_data_dict_0': 0.8326575756072998, 'mapped_reconstructed_data_0': 0.03913140296936035, 'reconstruction_reduced_0': 9.775161743164062e-06, 'regularization_matrix_reduced_0': 3.0994415283203125e-06, 'regularization_term_0': 0.0005021095275878906, 'log_det_curvature_reg_matrix_term_0': 0.011707544326782227, 'log_det_regularization_matrix_term_0': 0.011070489883422852, 'figure_of_merit_0': 0.972273588180542}

and

Log Likelihood Evaluation Time (second) = 192.92181277275085
Estimated Run Time Upper Limit (seconds) =  1326337.4628126621

From the samples_info.json:

   "log_evidence": 15868836.349132173,
    "total_samples": 26200,
    "time": "1279071.1350502968",
    "number_live_points": 150

source_pix[2] search.summary :

Total Samples = 13172
Total Accepted Samples = 2533
Acceptance Ratio = 0.19230185241421197
Time To Run = 11 days, 0:02:56.990573
Time Per Sample (seconds) = 72.16648880755905
Log Likelihood Function Evaluation Time (seconds) = 502.7288188934326
Expected Time To Run (seconds) = 76 days, 15:25:44.002464
Speed Up Factor (e.g. due to parallelization) = 6.966236368156963

from pyautolens.

Can you send the SLaM script you are running plz

from pyautolens.

I basically want to check if your SLaM run is using these settings, which seemingly sped up your code 1.45 seconds:

    transformer_class=al.TransformerNUFFT,
    use_linear_operators=True,

from pyautolens.

Main script:

"""
SLaM (Source, Light and Mass): Mass Total + Source Inversion
============================================================

SLaM pipelines break the analysis of 'galaxy-scale' strong lenses down into multiple pipelines which focus on modeling
a specific aspect of the strong lens, first the Source, then the (lens) Light and finally the Mass. Each of these
pipelines has it own inputs which customize the model and analysis in that pipeline.

The models fitted in earlier pipelines determine the model used in later pipelines. For example, if the SOURCE PIPELINE
uses a parametric `Sersic` profile for the bulge, this will be used in the subsequent MASS TOTAL PIPELINE.

Using a SOURCE LP PIPELINE, SOURCE PIX PIPELINE and a MASS TOTAL PIPELINE this SLaM script fits `Interferometer` 
of a strong lens system, where in the final model:

 - The lens galaxy's light is omitted from the data and model.
 - The lens galaxy's total mass distribution is an `PowerLaw`.
 - The source galaxy is an `Inversion`.

This uses the SLaM pipelines:

 `source_lp`
 `source_pix`
 `mass_total`

Check them out for a detailed description of the analysis!

__Run Times and Settings__

The run times of an interferometer `Inversion` depend significantly on the following settings:

 - `transformer_class`: whether a discrete Fourier transform (`TransformerDFT`) or non-uniform fast Fourier Transform
 (`TransformerNUFFT) is used to map the inversion's image from real-space to Fourier space.

 - `use_linear_operators`: whether the linear operator formalism or matrix formalism is used for the linear algebra.

The optimal settings depend on the number of visibilities in the dataset:

 - For N_visibilities < 1000: `transformer_class=TransformerDFT` and `use_linear_operators=False` gives the fastest
 run-times.
 - For  N_visibilities > ~10000: use `transformer_class=TransformerNUFFT`  and `use_linear_operators=True`.

The dataset modeled by default in this script has just 200 visibilties, therefore `transformer_class=TransformerDFT`
and `use_linear_operators=False`.

The script `autolens_workspace/*/interferometer/run_times.py` allows you to compute the run-time of an inversion
for your interferometer dataset. It does this for all possible combinations of settings and therefore can tell you
which settings give the fastest run times for your dataset.

__Start Here Notebook__

If any code in this script is unclear, refer to the `autolens_workspace/imaging/advanced/chaining/slam/start_here.ipynb`
notebook.
"""
# %matplotlib inline
# from pyprojroot import here
# workspace_path = str(here())
# %cd $workspace_path
# print(f"Working Directory has been set to `{workspace_path}`")

import os
import sys
from os import path
import autofit as af
import autolens as al
import autolens.plot as aplt

sys.path.insert(0, os.getcwd())
import slam

"""
__Dataset + Masking__ 

Load the `Interferometer` data, define the visibility and real-space masks and plot them.
"""
real_space_mask = al.Mask2D.circular(shape_native=(1440, 1440), pixel_scales=0.023, radius=5.,sub_size=1)

dataset_name = "hers1_a_07_cont"

dataset_path = path.join("dataset", "interferometer", dataset_name)

settings_dataset = al.SettingsInterferometer(transformer_class=al.TransformerNUFFT)

dataset = al.Interferometer.from_fits(
    data_path=path.join(dataset_path, "data.fits"),
    noise_map_path=path.join(dataset_path, "noise_map.fits"),
    uv_wavelengths_path=path.join(dataset_path, "uv_wavelengths.fits"),
    real_space_mask=real_space_mask,
)


dataset = dataset.apply_settings(settings=al.SettingsInterferometer(transformer_class=al.TransformerNUFFT))

positions = al.Grid2DIrregular(values=[(1.794, 2.116), (-1.127, -1.495), (0.0942,3.11), (-0.59,3.24)])

visuals = aplt.Visuals2D(positions=positions)


"""
__Inversion Settings (Run Times)__

The run times of an interferometer `Inversion` depend significantly on the following settings:

 - `transformer_class`: whether a discrete Fourier transform (`TransformerDFT`) or non-uniform fast Fourier Transform
 (`TransformerNUFFT) is used to map the inversion's image from real-space to Fourier space.

 - `use_linear_operators`: whether the linear operator formalism or matrix formalism is used for the linear algebra.

The optimal settings depend on the number of visibilities in the dataset:

 - For N_visibilities < 1000: `transformer_class=TransformerDFT` and `use_linear_operators=False` gives the fastest
 run-times.
 - For  N_visibilities > ~10000: use `transformer_class=TransformerNUFFT`  and `use_linear_operators=True`.

The dataset modeled by default in this script has just 200 visibilties, therefore `transformer_class=TransformerDFT`
and `use_linear_operators=False`. If you are using this script to model your own dataset with a different number of
visibilities, you should update the options below accordingly.

The script `autolens_workspace/*/interferometer/run_times.py` allows you to compute the run-time of an inversion
for your interferometer dataset. It does this for all possible combinations of settings and therefore can tell you
which settings give the fastest run times for your dataset.
"""


settings_inversion = al.SettingsInversion(use_linear_operators=True,
	image_mesh_min_mesh_pixels_per_pixel=3,
    image_mesh_min_mesh_number=5,
    image_mesh_adapt_background_percent_threshold=0.1,
    image_mesh_adapt_background_percent_check=0.8,)

"""
We now create the `Interferometer` object which is used to fit the lens model.

This includes a `SettingsInterferometer`, which includes the method used to Fourier transform the real-space 
image of the strong lens to the uv-plane and compare directly to the visiblities. We use a non-uniform fast Fourier 
transform, which is the most efficient method for interferometer datasets containing ~1-10 million visibilities.
"""

dataset = dataset.apply_settings(settings=settings_dataset)
dataset_plotter = aplt.InterferometerPlotter(dataset=dataset, visuals_2d=visuals)
dataset_plotter.subplot_dataset()
#dataset_plotter.subplot_dirty_images()

"""
__Settings AutoFit__

The settings of autofit, which controls the output paths, parallelization, database use, etc.
"""
settings_search = af.SettingsSearch(
    path_prefix=path.join("interferometer", "slam"),
    unique_tag=dataset_name,
    info=None,
    number_of_cores=16,
    session=None,
)

"""
__Redshifts__

The redshifts of the lens and source galaxies, which are used to perform unit converions of the model and data (e.g. 
from arc-seconds to kiloparsecs, masses to solar masses, etc.).
"""
redshift_lens_0 = 0.201854
redshift_lens_1 = 0.201854
redshift_source = 2.55293


"""
__SOURCE LP PIPELINE__

The SOURCE LP PIPELINE uses one search to initialize a robust model for the source galaxy's light, which in
this example:

 - Uses a parametric `Sersic` bulge for the source's light (omitting a disk / envelope).
 - Uses an `Isothermal` model for the lens's total mass distribution with an `ExternalShear`.
 
__Settings__:
 
 - Mass Centre: Fix the mass profile centre to (0.0, 0.0) (this assumption will be relaxed in the SOURCE INVERSION 
 PIPELINE).
"""

positions_likelihood = al.PositionsLHPenalty(positions=positions, threshold=1.0)

analysis = al.AnalysisInterferometer(dataset=dataset,positions_likelihood=positions_likelihood, settings_inversion=al.SettingsInversion(use_linear_operators=True))

source_lp_results = slam.source_lp2Lenses.run(
    settings_search=settings_search,
    analysis=analysis,
    lens_bulge_0=None,
    lens_disk_0=None,
    lens_bulge_1=None,
    lens_disk_1=None,
    mass_0=af.Model(al.mp.Isothermal),
    mass_1=af.Model(al.mp.IsothermalSph),
    einstein_radius_0 =  2.4868,
    einstein_radius_1 = 0.1696,
    shear=af.Model(al.mp.ExternalShear),
    source_bulge=af.Model(al.lp.Sersic),
    mass_centre_0=(-0.532, 0.287),
    mass_centre_1=(-2.83, 0.8957),
    redshift_lens_0=redshift_lens_0,
    redshift_lens_1=redshift_lens_1,
    redshift_source=redshift_source,
)

"""
__SOURCE PIX PIPELINE__

The SOURCE PIX PIPELINE uses two searches to initialize a robust model for the `Pixelization` that
reconstructs the source galaxy's light. It begins by fitting an `Overlay` image-mesh, `Delaunay` mesh and `Constant` 
regularization, to set up the model and hyper images, and then:

- Uses a `Hilbert` image-mesh. 
- Uses a `Delaunay` mesh.
 - Uses an `AdaptiveBrightness` regularization.
 - Carries the lens redshift, source redshift and `ExternalShear` of the SOURCE LP PIPELINE through to the
 SOURCE PIX PIPELINE.

__Settings__:

 - Positions: We update the positions and positions threshold using the previous model-fitting result (as described 
 in `chaining/examples/parametric_to_pixelization.py`) to remove unphysical solutions from the `Inversion` model-fitting.
"""

settings_inversion = al.SettingsInversion(use_linear_operators=True,
	image_mesh_min_mesh_pixels_per_pixel=3,
    image_mesh_min_mesh_number=5,
    image_mesh_adapt_background_percent_threshold=0.1,
    image_mesh_adapt_background_percent_check=0.8,)
    
    
analysis = al.AnalysisInterferometer(
    dataset=dataset,
    positions_likelihood=source_lp_results.last.positions_likelihood_from(
        factor=3.0, minimum_threshold=0.2
    ),
    settings_inversion=settings_inversion,
)

source_pix_results = slam.source_pix2Lenses.run(
    settings_search=settings_search,
    analysis=analysis,
    source_lp_results=source_lp_results,
    image_mesh=al.image_mesh.Hilbert,
    mesh=al.mesh.VoronoiNN,
    regularization=al.reg.AdaptiveBrightnessSplit,
)

"""
__MASS TOTAL PIPELINE__

The MASS TOTAL PIPELINE uses one search to fits a complex lens mass model to a high level of accuracy, 
using the lens mass model and source model of the SOURCE PIX PIPELINE to initialize the model priors. In this 
example it:

 - Uses an `PowerLaw` model for the lens's total mass distribution [The centre if unfixed from (0.0, 0.0)].
 - Uses a `Pixelization` for the source's light [fixed from SOURCE PIX PIPELINE].
 - Carries the lens redshift, source redshift and `ExternalShear` of the SOURCE PIPELINES through to the MASS 
 PIPELINE.
 
__Settings__:

 - adapt: We may be using adapt features and therefore pass the result of the SOURCE PIX PIPELINE to use as the
 hyper dataset if required.

 - Positions: We update the positions and positions threshold using the previous model-fitting result (as described 
 in `chaining/examples/parametric_to_pixelization.py`) to remove unphysical solutions from the `Inversion` model-fitting.
"""
analysis = al.AnalysisInterferometer(
    dataset=dataset,
    adapt_images=source_pix_results[0].adapt_images_from(use_model_image=True),
    positions_likelihood=source_pix_results.last.positions_likelihood_from(
        factor=3.0, minimum_threshold=0.2
    ),
    settings_inversion=settings_inversion,
)

mass_results = slam.mass_total2Lenses.run(
    settings_search=settings_search,
    analysis=analysis,
    source_results=source_pix_results,
    light_results=None,
    mass_0=af.Model(al.mp.PowerLaw),
    mass_1=af.Model(al.mp.PowerLaw),
)

"""
Finish.
"""

from pyautolens.

source_lp2Lenses script:

import autofit as af
import autolens as al


from typing import Callable, Union, Optional, Tuple


def run(
    settings_search: af.SettingsSearch,
    analysis: Union[al.AnalysisImaging, al.AnalysisInterferometer],
    lens_bulge_0: Optional[af.Model] = af.Model(al.lp.Sersic),
    lens_disk_0: Optional[af.Model] = af.Model(al.lp.Exponential),
    lens_bulge_1: Optional[af.Model] = af.Model(al.lp.Sersic),
    lens_disk_1: Optional[af.Model] = af.Model(al.lp.Exponential),
    mass_0: af.Model = af.Model(al.mp.Isothermal),
    mass_1: af.Model = af.Model(al.mp.IsothermalSph),
    einstein_radius_0: Optional[float] = None,
    einstein_radius_1: Optional[float] = None,
    shear: af.Model(al.mp.ExternalShear) = af.Model(al.mp.ExternalShear),
    source_bulge: Optional[af.Model] = af.Model(al.lp.Sersic),
    source_disk: Optional[af.Model] = None,
    redshift_lens_0: float = 0.5,
    redshift_lens_1: float = 0.5,
    redshift_source: float = 1.0,
    mass_centre_0: Optional[Tuple[float, float]] = None,
    mass_centre_1: Optional[Tuple[float, float]] = None,
    clump_model: Union[al.ClumpModel, al.ClumpModelDisabled] = al.ClumpModelDisabled(),
) -> af.ResultsCollection:
    """
    The SlaM SOURCE LP PIPELINE, which provides an initial model for the lens's light, mass and source using a
    parametric source model (e.g. Sersics, an MGE).

    Parameters
    ----------
    analysis
        The analysis class which includes the `log_likelihood_function` and can be customized for the SLaM model-fit.
    lens_bulge
        The model used to represent the light distribution of the lens galaxy's bulge (set to
        None to omit a bulge).
    lens_disk
        The model used to represent the light distribution of the lens galaxy's disk (set to
        None to omit a disk).
    mass
        The `MassProfile` fitted by this pipeline.
    shear
        The model used to represent the external shear in the mass model (set to None to turn off shear).
    source_bulge
        The model used to represent the light distribution of the source galaxy's bulge (set to
        None to omit a bulge).
    source_disk
        The model used to represent the light distribution of the source galaxy's disk (set to
        None to omit a disk).
    redshift_lens
        The redshift of the lens galaxy fitted, used by the pipeline for converting arc-seconds to kpc, masses to
        solMass, etc.
    redshift_source
        The redshift of the source galaxy fitted, used by the pipeline for converting arc-seconds to kpc, masses to
        solMass, etc.
    mass_centre
       If input, a fixed (y,x) centre of the mass profile is used which is not treated as a free parameter by the
       non-linear search.
    clump_model
        Add additional clumps containing light and mass profiles to the lens model. These have a known input centre and
        are used to model nearby line of sight galaxies.
    """

    """
    __Model + Search + Analysis + Model-Fit (Search 1)__

    Search 1 of the SOURCE LP PIPELINE fits a lens model where:

     - The lens galaxy light is modeled using a light profiles [no prior initialization].
     - The lens galaxy mass is modeled using a total mass distribution [no prior initialization].
     - The source galaxy's light is a light profiles [no prior initialization].

    This search aims to accurately estimate an initial lens light model, mass model and source model.
    """

    if mass_centre_0 is not None:
        mass_0.centre = mass_centre_0
        
    if mass_centre_1 is not None:
        mass_1.centre = mass_centre_1
        
    if einstein_radius_0 is not None:
        mass_0.einstein_radius = einstein_radius_0
        
    if einstein_radius_1 is not None:
        mass_1.einstein_radius = einstein_radius_1

    model_1 = af.Collection(
        galaxies=af.Collection(
            lens_0=af.Model(
                al.Galaxy,
                redshift=redshift_lens_0,
                bulge=lens_bulge_0,
                disk=lens_disk_0,
                mass=mass_0,
                shear=shear,
            ),
            lens_1=af.Model(
                al.Galaxy,
                redshift=redshift_lens_1,
                bulge=lens_bulge_1,
                disk=lens_disk_1,
                mass=mass_1,
            ),
            source=af.Model(
                al.Galaxy,
                redshift=redshift_source,
                bulge=source_bulge,
                disk=source_disk,
            ),
        ),
        clumps=clump_model.clumps,
    )
    

    search_1 = af.Nautilus(
        name="source_lp[1]_light[lp]_mass[total]_source[lp]",
        **settings_search.search_dict,
        n_live=200,
    )

    result_1 = search_1.fit(
        model=model_1, analysis=analysis, **settings_search.fit_dict
    )

    return af.ResultsCollection([result_1])

from pyautolens.

And lastly source_pix2Lenses script:

import autofit as af
import autolens as al

from typing import Optional, Tuple, Union


def run(
    settings_search: af.SettingsSearch,
    analysis: Union[al.AnalysisImaging, al.AnalysisInterferometer],
    source_lp_results: af.ResultsCollection,
    image_mesh_init: af.Model(al.AbstractImageMesh) = af.Model(al.image_mesh.Overlay),
    mesh_init: af.Model(al.AbstractMesh) = af.Model(al.mesh.VoronoiNN),
    image_mesh_init_shape: Tuple[int, int] = (34, 34),
    regularization_init: af.Model(al.AbstractRegularization) = af.Model(
        al.reg.AdaptiveBrightnessSplit
    ),
    image_mesh: af.Model(al.AbstractImageMesh) = af.Model(al.image_mesh.Hilbert),
    mesh: af.Model(al.AbstractMesh) = af.Model(al.mesh.VoronoiNN),
    regularization: af.Model(al.AbstractRegularization) = af.Model(
        al.reg.AdaptiveBrightnessSplit
    ),
    image_mesh_pixels_fixed: Optional[int] = 1000,
) -> af.ResultsCollection:
    """
    The SLaM SOURCE PIX PIPELINE, which initializes a lens model which uses a pixelized source for the source
    analysis.

    Parameters
    ----------
    analysis
        The analysis class which includes the `log_likelihood_function` and can be customized for the SLaM model-fit.
    source_lp_results
        The results of the SLaM SOURCE LP PIPELINE which ran before this pipeline.
    image_mesh_init
        The image mesh, which defines how the mesh centres are computed in the image-plane, used by the pixelization
        in the first search which initializes the source.
    image_mesh_init_shape
        The shape (e.g. resolution) of the image-mesh used in the initialization search (`search[1]`). This is only
        used if the image-mesh has a `shape` parameter (e.g. `Overlay`).
    mesh_init
        The mesh, which defines how the source is reconstruction in the source-plane, used by the pixelization
        in the first search which initializes the source.
    regularization_init
        The regularization, which places a smoothness prior on the source reconstruction, used by the pixelization
        which fits the source light in the initialization search (`search[1]`).
    image_mesh
        The image mesh, which defines how the mesh centres are computed in the image-plane, used by the pixelization
        in the final search which improves the source adaption.
    mesh
        The mesh, which defines how the source is reconstruction in the source-plane, used by the pixelization
        in the final search which improves the source adaption.
    regularization
        The regularization, which places a smoothness prior on the source reconstruction, used by the pixelization
        in the final search which improves the source adaption.
    image_mesh_pixels_fixed
        The fixed number of pixels in the image-mesh, if an image-mesh with an input number of pixels is used
        (e.g. `Hilbert`).
    """

    """
    __Model + Search + Analysis + Model-Fit (Search 1)__

    Search 1 of the SOURCE PIX PIPELINE fits a lens model where:

    - The lens galaxy light is modeled using a light profiles [parameters fixed to result of SOURCE LP PIPELINE].

     - The lens galaxy mass is modeled using a total mass distribution [parameters initialized from the results of the 
     SOURCE LP PIPELINE].

     - The source galaxy's light is the input initialization imagemesh, mesh and regularization scheme [parameters of 
     regularization free to vary].

    This search improves the lens mass model by modeling the source using a `Pixelization` and computes the adapt
    images that are used in search 2.
    """

    analysis.adapt_images = source_lp_results.last.adapt_images_from(use_model_images=True)

    mass_0 = al.util.chaining.mass_from(
        mass=source_lp_results.last.model.galaxies.lens_0.mass,
        mass_result=source_lp_results.last.model.galaxies.lens_0.mass,
        unfix_mass_centre=True,
    )
    
    mass_1 = al.util.chaining.mass_from(
        mass=source_lp_results.last.model.galaxies.lens_1.mass,
        mass_result=source_lp_results.last.model.galaxies.lens_1.mass,
        unfix_mass_centre=True,
    )
    image_mesh_init.shape = image_mesh_init_shape

    model_1 = af.Collection(
        galaxies=af.Collection(
            lens_0=af.Model(
                al.Galaxy,
                redshift=source_lp_results.last.instance.galaxies.lens_0.redshift,
                bulge=source_lp_results.last.instance.galaxies.lens_0.bulge,
                disk=source_lp_results.last.instance.galaxies.lens_0.disk,
                mass=mass_0,
                shear=source_lp_results.last.model.galaxies.lens_0.shear,
            ),
            lens_1=af.Model(
                al.Galaxy,
                redshift=source_lp_results.last.instance.galaxies.lens_1.redshift,
                bulge=source_lp_results.last.instance.galaxies.lens_1.bulge,
                disk=source_lp_results.last.instance.galaxies.lens_1.disk,
                mass=mass_1,
            ),
            source=af.Model(
                al.Galaxy,
                redshift=source_lp_results.last.instance.galaxies.source.redshift,
                pixelization=af.Model(
                    al.Pixelization,
                    image_mesh=image_mesh_init,
                    mesh=mesh_init,
                    regularization=regularization_init,
                ),
            ),
        ),
        clumps=al.util.chaining.clumps_from(result=source_lp_results.last),
    )

    search_1 = af.Nautilus(
        name="source_pix[1]_light[fixed]_mass[init]_source[pix_init_mag]",
        **settings_search.search_dict,
        n_live=150,
    )

    result_1 = search_1.fit(
        model=model_1, analysis=analysis, **settings_search.fit_dict
    )

    """
    __Model + Search + Analysis + Model-Fit (Search 2)__

    Search 2 of the SOURCE PIX PIPELINE fits a lens model where:

    - The lens galaxy light is modeled using a light profiles [parameters fixed to result of SOURCE LP PIPELINE].
    - The lens galaxy mass is modeled using a total mass distribution [parameters fixed to result of search 2].
    - The source galaxy's light is the input final mesh and regularization.

    This search initializes the pixelization's mesh and regularization.
    """
    analysis_1 = al.AnalysisInterferometer(
        dataset=analysis.dataset,
	positions_likelihood=analysis.positions_likelihood,
        adapt_images=result_1.adapt_images_from(use_model_images=True),
        settings_inversion=al.SettingsInversion(
            image_mesh_min_mesh_pixels_per_pixel=3,
            image_mesh_min_mesh_number=5,
            image_mesh_adapt_background_percent_threshold=0.1,
            image_mesh_adapt_background_percent_check=0.8, use_linear_operators=True,
        ),
    )

    model_2 = af.Collection(
        galaxies=af.Collection(
            lens_0=af.Model(
                al.Galaxy,
                redshift=source_lp_results.last.instance.galaxies.lens_0.redshift,
                bulge=source_lp_results.last.instance.galaxies.lens_0.bulge,
                disk=source_lp_results.last.instance.galaxies.lens_0.disk,
                mass=result_1.instance.galaxies.lens_0.mass,
                shear=result_1.instance.galaxies.lens_0.shear,
            ),
            lens_1=af.Model(
                al.Galaxy,
                redshift=source_lp_results.last.instance.galaxies.lens_1.redshift,
                bulge=source_lp_results.last.instance.galaxies.lens_1.bulge,
                disk=source_lp_results.last.instance.galaxies.lens_1.disk,
                mass=result_1.instance.galaxies.lens_1.mass,
            ),
            source=af.Model(
                al.Galaxy,
                redshift=source_lp_results.last.instance.galaxies.source.redshift,
                pixelization=af.Model(
                    al.Pixelization,
                    image_mesh=image_mesh,
                    mesh=mesh,
                    regularization=regularization,
                ),
            ),
        ),
        clumps=al.util.chaining.clumps_from(result=source_lp_results.last),
    )

    if image_mesh_pixels_fixed is not None:
        if hasattr(model_2.galaxies.source.pixelization.image_mesh, "pixels"):
            model_2.galaxies.source.pixelization.image_mesh.pixels = (
                image_mesh_pixels_fixed
            )

    search_2 = af.DynestyStatic(name="source_pix[2]_light[fixed]_mass[fixed]_source[pix]",
        **settings_search.search_dict,
        nlive=100,
    )

    result_2 = search_2.fit(
        model=model_2, analysis=analysis_1, **settings_search.fit_dict
    )

    return af.ResultsCollection([result_1, result_2])

from pyautolens.

Finally, can you send me the script where you get a Fit Time of 1.45?

from pyautolens.

I have overwritten it by now but I should be able to recreate it relatively quickly.
The edit will come in ~30min.

from pyautolens.

Great.

Your scripts look correct (e.g. settings_inversion) and TransformNUFFT are being supplied in the way I expected.

I am now wondering if when you did the profiling test to get 1.45 seconds, your real_space_mask used very different setings to the one in your script currently.

Your real space mask is very high resolution (e.g. it has a lot of image pixels), so the slow run times are not that surprising to me.

from pyautolens.

Ok, i was able to recreate it:

Transformer Time = 0.8331937789916992

15158617.926531568
Fit Time = 1.4778831005096436

From the following code source_pix2Lenses:

import autofit as af
import autolens as al
import numpy as np
from os import path
import time
from typing import Optional, Tuple, Union
import autofit as af
import autolens as al
import autolens.plot as aplt

def run(
    settings_search: af.SettingsSearch,
    analysis: Union[al.AnalysisImaging, al.AnalysisInterferometer],
    source_lp_results: af.ResultsCollection,
    image_mesh_init: af.Model(al.AbstractImageMesh) = af.Model(al.image_mesh.Overlay),
    mesh_init: af.Model(al.AbstractMesh) = af.Model(al.mesh.VoronoiNN),
    image_mesh_init_shape: Tuple[int, int] = (34, 34),
    regularization_init: af.Model(al.AbstractRegularization) = af.Model(
        al.reg.AdaptiveBrightnessSplit
    ),
    image_mesh: af.Model(al.AbstractImageMesh) = af.Model(al.image_mesh.Hilbert),
    mesh: af.Model(al.AbstractMesh) = af.Model(al.mesh.VoronoiNN),
    regularization: af.Model(al.AbstractRegularization) = af.Model(
        al.reg.AdaptiveBrightnessSplit
    ),
    image_mesh_pixels_fixed: Optional[int] = 1000,
) -> af.ResultsCollection:
    """
    The SLaM SOURCE PIX PIPELINE, which initializes a lens model which uses a pixelized source for the source
    analysis.

    Parameters
    ----------
    analysis
        The analysis class which includes the `log_likelihood_function` and can be customized for the SLaM model-fit.
    source_lp_results
        The results of the SLaM SOURCE LP PIPELINE which ran before this pipeline.
    image_mesh_init
        The image mesh, which defines how the mesh centres are computed in the image-plane, used by the pixelization
        in the first search which initializes the source.
    image_mesh_init_shape
        The shape (e.g. resolution) of the image-mesh used in the initialization search (`search[1]`). This is only
        used if the image-mesh has a `shape` parameter (e.g. `Overlay`).
    mesh_init
        The mesh, which defines how the source is reconstruction in the source-plane, used by the pixelization
        in the first search which initializes the source.
    regularization_init
        The regularization, which places a smoothness prior on the source reconstruction, used by the pixelization
        which fits the source light in the initialization search (`search[1]`).
    image_mesh
        The image mesh, which defines how the mesh centres are computed in the image-plane, used by the pixelization
        in the final search which improves the source adaption.
    mesh
        The mesh, which defines how the source is reconstruction in the source-plane, used by the pixelization
        in the final search which improves the source adaption.
    regularization
        The regularization, which places a smoothness prior on the source reconstruction, used by the pixelization
        in the final search which improves the source adaption.
    image_mesh_pixels_fixed
        The fixed number of pixels in the image-mesh, if an image-mesh with an input number of pixels is used
        (e.g. `Hilbert`).
    """

    """
    __Model + Search + Analysis + Model-Fit (Search 1)__

    Search 1 of the SOURCE PIX PIPELINE fits a lens model where:

    - The lens galaxy light is modeled using a light profiles [parameters fixed to result of SOURCE LP PIPELINE].

     - The lens galaxy mass is modeled using a total mass distribution [parameters initialized from the results of the 
     SOURCE LP PIPELINE].

     - The source galaxy's light is the input initialization imagemesh, mesh and regularization scheme [parameters of 
     regularization free to vary].

    This search improves the lens mass model by modeling the source using a `Pixelization` and computes the adapt
    images that are used in search 2.
    """

    analysis.adapt_images = source_lp_results.last.adapt_images_from(use_model_images=True)

    mass_0 = al.util.chaining.mass_from(
        mass=source_lp_results.last.model.galaxies.lens_0.mass,
        mass_result=source_lp_results.last.model.galaxies.lens_0.mass,
        unfix_mass_centre=True,
    )
    
    mass_1 = al.util.chaining.mass_from(
        mass=source_lp_results.last.model.galaxies.lens_1.mass,
        mass_result=source_lp_results.last.model.galaxies.lens_1.mass,
        unfix_mass_centre=True,
    )
    image_mesh_init.shape = image_mesh_init_shape

    model_1 = af.Collection(
        galaxies=af.Collection(
            lens_0=af.Model(
                al.Galaxy,
                redshift=source_lp_results.last.instance.galaxies.lens_0.redshift,
                bulge=source_lp_results.last.instance.galaxies.lens_0.bulge,
                disk=source_lp_results.last.instance.galaxies.lens_0.disk,
                mass=mass_0,
                shear=source_lp_results.last.model.galaxies.lens_0.shear,
            ),
            lens_1=af.Model(
                al.Galaxy,
                redshift=source_lp_results.last.instance.galaxies.lens_1.redshift,
                bulge=source_lp_results.last.instance.galaxies.lens_1.bulge,
                disk=source_lp_results.last.instance.galaxies.lens_1.disk,
                mass=mass_1,
            ),
            source=af.Model(
                al.Galaxy,
                redshift=source_lp_results.last.instance.galaxies.source.redshift,
                pixelization=af.Model(
                    al.Pixelization,
                    image_mesh=image_mesh_init,
                    mesh=mesh_init,
                    regularization=regularization_init,
                ),
            ),
        ),
        clumps=al.util.chaining.clumps_from(result=source_lp_results.last),
    )

    search_1 = af.Nautilus(
        name="source_pix[1]_light[fixed]_mass[init]_source[pix_init_mag]",
        **settings_search.search_dict,
        n_live=150,
    )

    result_1 = search_1.fit(
        model=model_1, analysis=analysis, **settings_search.fit_dict
    )
    


    """
    __Model + Search + Analysis + Model-Fit (Search 2)__

    Search 2 of the SOURCE PIX PIPELINE fits a lens model where:

    - The lens galaxy light is modeled using a light profiles [parameters fixed to result of SOURCE LP PIPELINE].
    - The lens galaxy mass is modeled using a total mass distribution [parameters fixed to result of search 2].
    - The source galaxy's light is the input final mesh and regularization.

    This search initializes the pixelization's mesh and regularization.
    """
    real_space_mask = al.Mask2D.circular(shape_native=(1440, 1440), pixel_scales=0.023, radius=5.,sub_size=1)

    dataset_name = "hers1_a_07_cont"

    dataset_path = path.join("dataset", "interferometer", dataset_name)
    
    settings_dataset = al.SettingsInterferometer(transformer_class=al.TransformerNUFFT)
    
    dataset = al.Interferometer.from_fits(
        data_path=path.join(dataset_path, "data.fits"),
        noise_map_path=path.join(dataset_path, "noise_map.fits"),
        uv_wavelengths_path=path.join(dataset_path, "uv_wavelengths.fits"),
        real_space_mask=real_space_mask,
    )
    
    
    dataset = dataset.apply_settings(settings=al.SettingsInterferometer(transformer_class=al.TransformerNUFFT)) 
    
    settings_inversion = al.SettingsInversion(use_linear_operators=True,
            image_mesh_min_mesh_pixels_per_pixel=3,
        image_mesh_min_mesh_number=5,
        image_mesh_adapt_background_percent_threshold=0.1,
        image_mesh_adapt_background_percent_check=0.8,)

    analysis_1 = al.AnalysisInterferometer(
        dataset=dataset,
	positions_likelihood=analysis.positions_likelihood,
        adapt_images=result_1.adapt_images_from(use_model_images=True),
        settings_inversion=settings_inversion
    )

    model_2 = af.Collection(
        galaxies=af.Collection(
            lens_0=af.Model(
                al.Galaxy,
                redshift=source_lp_results.last.instance.galaxies.lens_0.redshift,
                bulge=source_lp_results.last.instance.galaxies.lens_0.bulge,
                disk=source_lp_results.last.instance.galaxies.lens_0.disk,
                mass=result_1.instance.galaxies.lens_0.mass,
                shear=result_1.instance.galaxies.lens_0.shear,
            ),
            lens_1=af.Model(
                al.Galaxy,
                redshift=source_lp_results.last.instance.galaxies.lens_1.redshift,
                bulge=source_lp_results.last.instance.galaxies.lens_1.bulge,
                disk=source_lp_results.last.instance.galaxies.lens_1.disk,
                mass=result_1.instance.galaxies.lens_1.mass,
            ),
            source=af.Model(
                al.Galaxy,
                redshift=source_lp_results.last.instance.galaxies.source.redshift,
                pixelization=af.Model(
                    al.Pixelization,
                    image_mesh=image_mesh,
                    mesh=mesh,
                    regularization=regularization,
                ),
            ),
        ),
        clumps=al.util.chaining.clumps_from(result=source_lp_results.last),
    )

    if image_mesh_pixels_fixed is not None:
        if hasattr(model_2.galaxies.source.pixelization.image_mesh, "pixels"):
            model_2.galaxies.source.pixelization.image_mesh.pixels = (
                image_mesh_pixels_fixed
            )
            
    lens_0=al.Galaxy(
        redshift=source_lp_results.last.instance.galaxies.lens_0.redshift,
        bulge=source_lp_results.last.instance.galaxies.lens_0.bulge,
        disk=source_lp_results.last.instance.galaxies.lens_0.disk,
        mass=result_1.instance.galaxies.lens_0.mass,
        shear=result_1.instance.galaxies.lens_0.shear)
    
    lens_1=al.Galaxy(
        redshift=source_lp_results.last.instance.galaxies.lens_1.redshift,
        bulge=source_lp_results.last.instance.galaxies.lens_1.bulge,
        disk=source_lp_results.last.instance.galaxies.lens_1.disk,
        mass=result_1.instance.galaxies.lens_1.mass)
    source=al.Galaxy(
        redshift=source_lp_results.last.instance.galaxies.source.redshift,
        pixelization=af.Model(
            al.Pixelization,
            image_mesh=image_mesh,
            mesh=mesh,
            regularization=regularization,))
        
    tracer=al.Tracer.from_galaxies(galaxies=[lens_0,lens_1, source])
        
    
    def print_transformer_time_from(dataset, transformer_class, repeats=1):
        settings_dataset = al.SettingsInterferometer(transformer_class=transformer_class)
        dataset = dataset.apply_settings(settings=settings_dataset)
    
        """
        __Numba Caching__
    
        Perform a single transformer call to ensure all numba functions are initialized.
        """
        image = tracer.image_2d_from(grid=dataset.grid)
    
        dataset.transformer.visibilities_from(image=image)
    
        """
        __Fit Time__
    
        Now profile the overall run-time of the transformer.
        """
        start = time.time()
        dataset.transformer.visibilities_from(image=image)
    
        transformer_time = (time.time() - start) / repeats
        print(f"Transformer Time = {transformer_time} \n")
    
    
    """
    __Fit Time__
    
    This function is used throughout this script to time how long a fit takes for each combination of settings.
    """
    
    
    def print_fit_time_from(dataset, transformer_class, use_linear_operators, repeats=1):
        settings_dataset = al.SettingsInterferometer(transformer_class=transformer_class)
        dataset = dataset.apply_settings(settings=settings_dataset)
    
        """
        __Numba Caching__
    
        Call FitImaging once to get all numba functions initialized.
        """
        fit = al.FitInterferometer(
            dataset=dataset,
            tracer=tracer,
            settings_inversion=al.SettingsInversion(
                use_linear_operators=use_linear_operators
            ),
        )
        print(fit.figure_of_merit)
    
        """
        __Fit Time__
    
        Time FitImaging by itself, to compare to run_times dict call.
        """
        start = time.time()
        for i in range(repeats):
            fit = al.FitInterferometer(
                dataset=dataset,
                tracer=tracer,
                settings_inversion=al.SettingsInversion(
                    use_linear_operators=use_linear_operators
                ),
            )
            fit.figure_of_merit
    
        fit_time = (time.time() - start) / repeats
        print(f"Fit Time = {fit_time} \n")
        
    print_transformer_time_from(
        dataset=dataset, transformer_class=al.TransformerNUFFT, repeats=1
    )
        
    print_fit_time_from(
        dataset=dataset,
        transformer_class=al.TransformerNUFFT,
        use_linear_operators=True,
    )
        
    search_2 = af.DynestyStatic(name="source_pix[2]_light[fixed]_mass[fixed]_source[pix]",
        **settings_search.search_dict,
        nlive=100,
    )
   
    result_2 = search_2.fit(
        model=model_2, analysis=analysis_1, **settings_search.fit_dict
    )

    return af.ResultsCollection([result_1, result_2])

from pyautolens.

I suspect the fit time being 1.45 seconds is due to an input error, albeit tracing the exact root cause would be a bit of a nightmare.

Are you able to simply input your data and mask into this function with your settings to get the fit time? It does not matter if things like the lens_galaxy are empty for the sake of profiling:

https://github.com/Jammy2211/autolens_workspace/blob/release/scripts/interferometer/data_preparation/examples/run_times.py

Obviously if you get a linear algebra error then this is no good.

from pyautolens.

Not sure if I can directly translate this. For the Hilber image mesh, I need an adapt image that I don't have when just checking the run time with this script.

I can give you the results for an overlay with shape (34,34) for now.

from pyautolens.

Yep, Overlay is fine.

from pyautolens.

Results from the run_times.py script are in:

Transformer Time = 0.8737094402313232

15849865.005819395
Fit Time = 193.9684615135193

Fit time is long, as should be, i suppose.

from pyautolens.

Ok, at least we have confirmed things are behaving as expected now.

The fit times of 1.45 were due to an input bug, I sort of know the problem (its to do with af.Model()) but probably not worth digging deeper into.

I anticipate your run times basically scale directly with number of visilibites (albeit you may as well quickly check that changing real_space_mask doesnt drive the fit time.

At this point, I think your options are to either reduce the number of visilibites to do the fit, or bite the bullet and accept these extremely long run times.

(I dont know what your science case is, but you can probably devise a strategy which fits with a reduced number of visibilities and then reconstructs the data and source with the full amount, using the max LH lens model).

from pyautolens.

Thank you for the support and explanations.

The fit times of 1.45 were due to an input bug, I sort of know the problem (its to do with af.Model()) but probably not worth digging deeper into.

Good to know, should I encounter this again sometime.

I'll play around and see where I can optimize.

(I dont know what your science case is, but you can probably devise a strategy which fits with a reduced number of visibilities and then reconstructs the data and source with the full amount, using the max LH lens model).

This sounds like a possible way forward, yes. That would mean doing a parametric and the first pixelation run with the reduced number of visibilities and then using that lens model for the second pixelization with the "full" number of visibilities, correct? Basically, I need to resolve as much as possible and follow the differential magnification of the image back to the spatial distribution of the source. Most likely, I will be able to use up to 32 cores and do a bunch of runs in parallel on different nodes on the hpc that I am using. This will be necessary since I want to do this for four galaxies with multiple continuum datasets and for a couple of emission/absorption line bins. Though for most of them, I can use the max LH lens model from the previous fits.

from pyautolens.

This sounds like a possible way forward, yes. That would mean doing a parametric and the first pixelation run with the reduced number of visibilities and then using that lens model for the second pixelization with the "full" number of visibilities, correct?

I think you need things to run fast all the way through to mass[1], so I would at least start by running SLaM end-to-end with reduced visibilitlies.

I would then, at the end, output the quantities you need using the full visibilties dataset, using the max LH lens model. You can then get on with figuring out your source analysis scripts using these reconstructions, and make sure everything is looking like it'll work fine.

In the mean time, whilst doing that, you can set off some source_pix[2] and mass[1] runs using the full visibilities and I guess after a month they'll have completed lol.

from pyautolens.

Thank you for the advice! Seems like a good plan to move forward with for now.

from pyautolens.

I anticipate your run times basically scale directly with number of visilibites (albeit you may as well quickly check that changing real_space_mask doesnt drive the fit time.

So this comment has occupied me for a little while. The shape_native parameter in the realspace mask certainly does impact the fit time a bit. But it's definitely the visibilities:
Testing sets with the same circular mask:
shape_native=[1024,1024]
with ~5mil visibilities

Transformer Time = 0.6205770969390869 

15852338.853299808
Fit Time = 187.65951585769653 

And the polarization averaged version ~2.5mil visibilities:

Transformer Time = 0.39135289192199707 

9649472.928159023
Fit Time = 106.84008979797363 

~5mil visibilities
shape_native=[1440,1440]

Transformer Time = 0.6800346374511719 

15852338.888308465
Fit Time = 193.51543259620667 

~2.5mil visibilities

Transformer Time = 0.47258710861206055 

9649472.963167684
Fit Time = 121.33718395233154 

Another question I had, was related to the PyNUFFT paper as they use very detailed masking for their image.
See Figure 6 on the left of the PyNUFFT paper

I can already implement quite complex masks,

but my question is, can such a complex mask help pixelation and the run time in PyAutoLens?
As far as I understood it, reducing the number of unmasked pixels from ~10⁵ to 10⁴ should decrease the overall fitting time too.
At the moment it even increases the fit time slightly. shape_native=[1440,1440], ~2.5mil visibilities and advanced mask as in the screenshot above.

Transformer Time = 0.459104061126709 

9658926.659100045
Fit Time = 133.0904037952423 

Is this the correct behavior, or is the mask ignored somehow?

from pyautolens.

Yeah, so now that I've had a chance to think over it your profiling all looks correct, at high numbers of visibilities (e.g. > 1 million) everything is fast compared to the NUFFT, which is called many times for a pixelized source. We are working on a solution to this (see JAX lol) but its a whilst off.

but my question is, can such a complex mask help pixelation and the run time in PyAutoLens?
As far as I understood it, reducing the number of unmasked pixels from ~10⁵ to 10⁴ should decrease the overall fitting time too.
At the moment it even increases the fit time slightly. shape_native=[1440,1440], ~2.5mil visibilities and advanced mask as in the screenshot above.

I would expect that it does speed things up (surprising it didnt for you) but as you've shown above, in terms of run times its secondary to calculations which scale with the number of visilbities. So I'm not sure you gain much using a tigher mask.

Visibilities associated calculations are performed N times, where N is the number of source pixels. So reducing this might speed things up, but again I doubt it'll move you from the world of "really really slow".

from pyautolens.

Thank you for your answer.
Great to hear that JAX will help in that regard in the future.

I am still testing how much quicker it runs on the HPC when reducing the visibilities. With quite a bit of testing later, I was able to comfortably reduce the dataset to ~675000 visibilities without losing very important information. Amazingly, with the JAX fix I can finally just run the new reduced dataset with and without the advanced masking in parallel on different nodes. Having reduced the visibilities again, I got down to fit times of ~90s from the run_times.py test. So at the very least the fit time was cut in half.

Visibilities associated calculations are performed N times, where N is the number of source pixels. So reducing this might speed things up, but again I doubt it'll move you from the world of "really really slow".

So from additional testing, I have found that most of the time saved can be achieved through reducing the image mesh pixels. I wasn't sure you meant this when talking about the number of source pixels. Going from a shape of 32x32 to 16x16 did cut the total time in half. I will see how this applies also to the different datasets.

from pyautolens.

It looks like the better the mask is the lower I can make the image mesh shape and still reconstruct the source nicely.

from pyautolens.

Are you using the Hilbert image mesh from source_pix[2] onwards?

Ah, this requires a circular mask, so would require you to change to that.

The flip side is the Hilbert can cluster image pixels over the source in a much better way -- you could prob get to 750 source pixels and have great source coverage. I guess the problem is whether a circular mask becomes too slow.

from pyautolens.

Yes, I am using the Hilbert image mesh for source_pix[2] and then using those results for the total_mass.
Is it then possible to use a mask (i.e. done through the GUI) for everything up to source_pix[1] and then for the Hilbert image mesh apply the circular mask? For source_pix[1] circular mask vs GUI mask I save ~30s fit time if using the GUI mask.

from pyautolens.

I would stick to a single mask, whilst it may be possible thigns like the regularization parameters depend on the mask and trying to swap masks during a fit sounds error prone.

from pyautolens.

Between source_pix[1] and source_pix[2], only the lens model is transferred. In source_pix[2] we start with a completely new pixelization and only the adapt_images from the result of source_pix[1] are used, which are in the image plane anyway, if I understood correctly. So regularization should be free of the previous fit, no?

from pyautolens.

To summarize:
Averaging by a factor of 8 in the visibilities, I could improve the fit time in the source_pix[1] run (using the same settings) by a factor of 3.

Now, since I cannot further average the visibilities without losing important information for my science, we find that the only other factors significantly affecting run times are the real_space_mask and the image_mesh pixels.

For Hilbert image meshes, I have to use a circular mask. The shape_native parameter for the mask should be reduced while still sampling uv-space well enough. The number of pixels in the overlay mesh and the Hilbert mesh both affect fit times +~100-200s.

These three, I have to balance somehow and still get good source reconstructions and run times, correct?

Also fit time from run_times.py appears to be extremely inconsistent, sometimes showing differences of up to 100s.

from pyautolens.

Between source_pix[1] and source_pix[2], only the lens model is transferred. In source_pix[2] we start with a completely new pixelization and only the adapt_images from the result of source_pix[1] are used, which are in the image plane anyway, if I understood correctly. So regularization should be free of the previous fit, no?

I think what you are saying is correct.

The one problem I can imagine happening is that you compute an adapt_image in a small mask, meaning that when you switch to a larger mask all values further out are zeros, which could cause issues with the adaptive features. I'll be honest I'm not certain if this would be an issue or not, but as I said above it all sounds very error prone and worth taking the hit on CPU run time to avoid strange results.

To summarize:
Averaging by a factor of 8 in the visibilities, I could improve the fit time in the source_pix[1] run (using the same settings) by a factor of 3.

Now, since I cannot further average the visibilities without losing important information for my science, we find that the only other factors significantly affecting run times are the real_space_mask and the image_mesh pixels.

For Hilbert image meshes, I have to use a circular mask. The shape_native parameter for the mask should be reduced while still sampling uv-space well enough. The number of pixels in the overlay mesh and the Hilbert mesh both affect fit times +~100-200s.

These three, I have to balance somehow and still get good source reconstructions and run times, correct?

Yes, this sounds like a good summary. We really need to speed up the interferometer code 🤣

I would be inclined to suggest you go down the Hilbert route, as I think your source reconstructions will come out a lot better in the end and once you've polished it I would guess the reduction in source pixels will give a speed up which males up for the increase in image pixels.

Also fit time from run_times.py appears to be extremely inconsistent, sometimes showing differences of up to 100s.

Is this for different mass models? I would guess that run time can change a lot depending on what the mass model does.

from pyautolens.

The one problem I can imagine happening is that you compute an adapt_image in a small mask, meaning that when you switch to a larger mask all values further out are zeros, which could cause issues with the adaptive features. I'll be honest I'm not certain if this would be an issue or not, but as I said above it all sounds very error prone and worth taking the hit on CPU run time to avoid strange results.

Can you help me understand how adapt_image works for interferometric datasets? I thought that the adapt_image is the image plane best fit of the previous search that is then inverted again to visibilities, which are then used in the analysis object for the current run. So it should produce visibilities that we want to fit similarly to. However, once you invert these visibilities again we should automatically lose all information about the mask, since the resulting image has to be continuous (due to the overlapping fringes), no?
Connected question: How are masks realized for visibilities? From what I got, they constrain the available area for the image which is then transformed to the visibilities that are fit to the dataset. So the larger the mask, the more visibilities you can fit to the image, which can turn into a problem through overfitting the noise.

meaning that when you switch to a larger mask all values further out are zeros, which could cause issues with the adaptive features
Does the Hilbert curve algorithm work on the visibilities or the image plane? If it works on the image plane then is there a way to set all the zeros a larger circular mask would produce to the weight_floor value?

Yes, this sounds like a good summary. We really need to speed up the interferometer code 🤣

I wish it wasn't so, because this will be very difficult to balance. Me hoping that there is a solution for this somewhere 🙏 lol

I would be inclined to suggest you go down the Hilbert route, as I think your source reconstructions will come out a lot better in the end and once you've polished it I would guess the reduction in source pixels will give a speed up which makes up for the increase in image pixels.

Unfortunately, i won't get enough speed out of just reducing the Hilbert pixels. source_lp was quick at ~5h, source_pix[1] took ~20 days and source_pix[2] 22 days with the standard settings, with mass_total still running. Even if this is sped up by a factor of 3, we're still talking about more than a month. And there will be ~4-12 datasets that will have to go through this full pipeline. So at least also the real_space_mask needs adjusting (in shape and pixel_size).
P.S. I haven't forgotten about running a very reduced dataset and applying it to the full one. I'll will need to benchmark this technique along with just going the pure imaging route and the slam pipeline somehow.

Is this for different mass models? I would guess that run time can change a lot depending on what the mass model does.

Honestly, It might have been human error lol. It might have been related to using a fixed mass model and a free one, but I did so many tests and didn't document them in enough detail for this. Let's disregard this for now.

from pyautolens.

I finally bit the bullet and reduced the circular real_space_mask's shape_native and pixel_scale. Since those values are directly related to the uv grinding they also lead to significant runtime improvements. So now I am going with shape_native=[960,960] and pixel_shape=0.033 instead of what is written above. This is sampling the restoring beam by ~5 pixels/minor beam axis. Before it was ~7.5. If we compare how much of the FOV of the antennas we get when using the image size and pixel size, we still take shape_native to be ~1.5*FOV at this resolution. The reason for this is simply the ALMA utilities provide these values "au.pickCellSize()" with imsize=True for the full baselines. Possibly, I could even go down to shape_native=670,670 and I should still get the full FOV but this is for further testing.
Naturally, the fit time is still related to the radius of the circular real_space_mask and the pixelation image_mesh size. Keeping the image_mesh the same (32x32) but varying the radius appears to show a non linear effect on the fit time. I.e. 5" mask (73s) takes longer than an 8" mask (52s) but an 11" mask also takes longer (109s). On the other hand, with fixed mask radius (5.5") and varying mesh size, we find that the fit_time increases when increasing the mesh_size. Now, I am running the whole slam run with the 5.5" radius, the new shape_native, pixel_scales, image_mesh of 25x25 and Hilbert pixels set to 750. Let's see how this shapes up.

Did you already find a way to set the image_mesh that optimizes run_time and has the source be well enough sampled in the pixelated grid?

from pyautolens.

Can you help me understand how adapt_image works for interferometric datasets? I thought that the adapt_image is the image plane best fit of the previous search that is then inverted again to visibilities, which are then used in the analysis object for the current run. So it should produce visibilities that we want to fit similarly to. However, once you invert these visibilities again we should automatically lose all information about the mask, since the resulting image has to be continuous (due to the overlapping fringes), no?

The adapt_image for interferometer data is an image and not visibilitlies (its the `dirty_model_image). It is used to control the source pixelization and regularization -- both things which can be performed in real space. Thus, the adapt image is never converted to visilibiltiles.

Therefore, thinking about it, the adapt_image is tied to the real_space_mask, so thats the thing you cannot change safelty between fits.

Connected question: How are masks realized for visibilities? From what I got, they constrain the available area for the image which is then transformed to the visibilities that are fit to the dataset. So the larger the mask, the more visibilities you can fit to the image, which can turn into a problem through overfitting the noise.

I think this sounds correct, I'm not the interferometry expert so can't give you a full answer! I dont think larger masks will mean more noise over fitting, the mask essentially defines in real space where you are modeling the lens.

Unfortunately, i won't get enough speed out of just reducing the Hilbert pixels. source_lp was quick at ~5h, source_pix[1] took ~20 days and source_pix[2] 22 days with the standard settings, with mass_total still running. Even if this is sped up by a factor of 3, we're still talking about more than a month. And there will be ~4-12 datasets that will have to go through this full pipeline. So at least also the real_space_mask needs adjusting (in shape and pixel_size).

Yep, its not ideal, but interferometry is just very slow :(. At least you have a clear assessment of the run time options now, I will keep trying to think of a way for you to speed it up.

On the other hand, with fixed mask radius (5.5") and varying mesh size, we find that the fit_time increases when increasing the mesh_size. Now, I am running the whole slam run with the 5.5" radius, the new shape_native, pixel_scales, image_mesh of 25x25 and Hilbert pixels set to 750. Let's see how this shapes up.

Sounds good, it depends on your lens but 5.5" still sounds very large to me, I'm used to ALMA lenses before < 3.0" and therefore masks of that size being appropriate. What are the Einstein Radius of your lenses?

Did you already find a way to set the image_mesh that optimizes run_time and has the source be well enough sampled in the pixelated grid?

This isn't really possible, the image_mesh that has the source be well enough sampled depends on the lens being fitted (e.g. more compact sources with need higher resolution image meshs).

The Hilbert does this for you (e.g. it distributes your 750 source pixels in the best way possible to resolve the source). Even if your image_mesh in source_pix[1] is a bit rubbish, it'll do a good enough job for the Hilbert in source_pix[2] onwards.

Make sure you are using these settings in your SLaM pipeline, they ensure the Hilbert mesh adapts better:

    analysis = al.AnalysisImaging(
        ...
        settings_inversion=al.SettingsInversion(
            image_mesh_min_mesh_pixels_per_pixel=3,
            image_mesh_min_mesh_number=5,
            image_mesh_adapt_background_percent_threshold=0.1,
            image_mesh_adapt_background_percent_check=0.8,
        ),
    )

from pyautolens.

Thank you for extensive answers!

The adapt_image for interferometer data is an image and not visibilitlies (its the `dirty_model_image). It is used to control the source pixelization and regularization -- both things which can be performed in real space. Thus, the adapt image is never converted to visilibiltiles.

Thank you for the explanation. Is the adapt image then only used for creating the pixelation (grid?), i.e. to not pick up any spurious signal or noise we might get from using the previous dirty image? Wouldn't it be better to do this in uv space then, since the dirty image, model or data, always have correlated noise? Or is there a S/N cut you apply to not pick up the noise somehow?

Therefore, thinking about it, the adapt_image is tied to the real_space_mask, so thats the thing you cannot change safelty between fits.

Good to know!

I think this sounds correct, I'm not the interferometry expert so can't give you a full answer! I dont think larger masks will mean more noise over fitting, the mask essentially defines in real space where you are modeling the lens.

No worries, I am also still learning about it. Usually, when going from uv space to image space both image size and pixel size are related to the size of the uv grid and its cell size, i.e. how the data is sampled. So it it correct, that shape_native and pixel_scales define how we sample the uv space and radius only defines the area we are modeling the lens in? If this is true, then radius and the image_mesh size are two related quantities which can be reduced to get significant improvements in run time.

Sounds good, it depends on your lens but 5.5" still sounds very large to me, I'm used to ALMA lenses before < 3.0" and therefore masks of that size being appropriate. What are the Einstein Radius of your lenses?

The einstein radii are 2.4868" and 0.1686". And the image is a little bit off center, with a 4.0" circular mask at the center being sufficient in picking up all important flux. However, it was not good in run time to use this mask. Somehow a 5" or 5.5" mask worked better than 4 or 4.5" mask. Even 8" was pretty fast. In the end, I guess the run time then depends on the image mesh size.

This isn't really possible, the image_mesh that has the source be well enough sampled depends on the lens being fitted (e.g. more compact sources with need higher resolution image meshs).
The Hilbert does this for you (e.g. it distributes your 750 source pixels in the best way possible to resolve the source). Even if your image_mesh in source_pix[1] is a bit rubbish, it'll do a good enough job for the Hilbert in source_pix[2] onwards.

Good! This is what I assumed was the case. With your quoted settings, I should at least get the factor ?5? potentially more pixels for the bright regions, which should be plenty for a magnification of factor ~20.

Since I didn't mention it yet (didn't want to tempt fate lol), with the above changes I am looking at factor 5 improvement for the source_lp, and even a factor ~10 improvement in run time for source_pix[1]. That's already fantastic, as 2 days run time is so much better than 20 lol. But it isn't finished yet.

from pyautolens.

Thank you for the explanation. Is the adapt image then only used for creating the pixelation (grid?), i.e. to not pick up any spurious signal or noise we might get from using the previous dirty image? Wouldn't it be better to do this in uv space then, since the dirty image, model or data, always have correlated noise? Or is there a S/N cut you apply to not pick up the noise somehow?

The adapt image is used to:

• Weight regularization by source pixel brightness in the AdaptiveImageBrightness scheme (we will soon be switching the default regularization to one which does not require an adapt_image).
• Cluster points in the source's brightest regions in the Hilbert pixelization.

Both the tasks above are controlled by the parameters fitted for in source_pix[2], which means that the adapt_image doesn't need to be perfect and these parameters will adapt to mitigate any imperfections.

The adapt_image probably has a bit of correlated noise and other features, but it seems to be good enough to do it this way.

radius only defines the area we are modeling the lens in? If this is true, then radius and the image_mesh size are two related quantities which can be reduced to get significant improvements in run time.

This sounds correct. image_mesh size certainly has nothing to do with interferometer, its just how we distribute our source pixels in the image plane.

The einstein radii are 2.4868" and 0.1686". And the image is a little bit off center, with a 4.0" circular mask at the center being sufficient in picking up all important flux. However, it was not good in run time to use this mask. Somehow a 5" or 5.5" mask worked better than 4 or 4.5" mask. Even 8" was pretty fast. In the end, I guess the run time then depends on the image mesh size.

Strange, I guess FFT's run times depend on crazy things like whether certain arrays are powers of 2 lol.

You should probably try a smaller mask for the 0.1686" lens, feels like you're making the calculation slower than necessary.

Since I didn't mention it yet (didn't want to tempt fate lol), with the above changes I am looking at factor 5 improvement for the source_lp, and even a factor ~10 improvement in run time for source_pix[1]. That's already fantastic, as 2 days run time is so much better than 20 lol. But it isn't finished yet.

Yep, this sounds realistic. We have done a lot of careful crafting of our settings in the past, and this took us from months to days... so its defo doable! Apologies its such a nightmare though, and the docs are stlil a bit lacking.

from pyautolens.