border_pixels is a list of image pixels in the border, used to define the source-plane border. However, we don't necessarily want to trace image pixels to the source-plane, instead only using sub-pixels (and sparse pixels). So, we need a border_sub_pixels, which gives a list of the sub-pixel in each image pixel which is closest to the border.

I have made a visualize module, and am in the process of making individual functions for plotting every type of image.

Instead of manually passing the variables with dictate the matplotlib settings (e.g. figsize, colorscheme etc) we should make it so that each image can be customized using a config file. Im guessing we'll have a config file that looks like this:

[Image]
figsize=(10,15)
colorbar=linear
output_type=png

[Residuals]
figsize=(10,15)
colorbars=symmetric_log
output_type=fits

and so no - of course one should be able to skip these configs and use their manually input values

We are used to using weighted KMeans clustering to adapt the source-pixelization in the source plane. A consequence of this is we sample a unique discrete source-plane grid for every iteration of the method - something I believe helps remove systematic biases, but makes the method slower and more cumbersome. Its a shame we currently have no option to use a (nearly) fixed discreization, as this will speed things up a lot at the expense of introducing biases. However, we can switch to source-plane clustering at the end of the analysis, such that a long pipeline gains speed without bias.

To achieve this, we can perform weighted KMeans clutering on the source-weights in the image-plane (during the hyper-parameter optimization where we set up the source-plane pixelizlation, regularization etc). We then fix the clusters to the 'best-fit' points in the image plane for lens modeling, thereby mimicking the behaviour of a fixed source-plane discretization. This will reuse all the functions the source-plane clustering method uses and requires no extra methods, so should be simple to implement :).

This could also form the basis of a light-weighted Delaunay pixelization.

I have moved grids and masked_image (now called lensing_image) to a new module, called the lensing module, which also contains ray_tracing, galaxy, galaxy prior. Deciding where to put the grids module is a pain, but I opted for the lensing module as it is the grids that we use to perform lensing, and the choice of grid will now depend on the tracer we use.

For lens plane + source plane systems, we have a lot of freedom to customize our grid, as (unlike multi-planes) we are not held-back by the need to compute many intensities on non-uniform lensed grids. Therefore, I believe we should:

• Use an iterative grid when computing the image-plane (i.e. lens galaxy) intensities, which on-th-fly increase the sub-gridding level until a threshold numerical accuracy is reached for each light profile. This will ensure divergent Sersic profiles at low radii are resolved.

• A fixed sub-grid for deflection angle calculations and computation of the source-light profile / pixelized reconstruction. For pixelized reconstructions (the way we should be modelling sources anyway) we don't need an iterative grid to properly resolve the light profile of each pixel (our source light profile module will be a bit dodgy but oh well).

Once we have traced a grid to a source-plane, it is non-uniform. An iterative grid is no longer possible, as it relies on the regularity of the square-pixels. One could iteratively compute the deflection angles in the image-plane and then intensities in the source-plane, but this is a computation nightmare. Therefore, for multi-plane ray tracing I recommend we just stick to a regular sub-grid and take the hit on light-profile accuracy. I can't think of a better approach that doesn't get horribly complicated.

For simulated images, we need to ensure our images do reach numerical accuracy and thus need to perform iterative deflection + intensity calculations until a threshold accuracy is met. Luckily, for simulated images, we can easily do this for the entire tracer-instance at once for one sub-grid size, and then perform another tracer at higher sub-grid and compare them. This approach could work for actual modeling - lets see how it plays out for simulated images.

There are several different places that config is implemented for testing. Only one test config folder should exist.

Currently, one has to always manually pass the priors between two phases. This is good in that it provides full customization of the modeling, but has to major downsides:

1. Writing pipelines is cumbersome, because you have to specify how every prior links together.
2. It is not possible to link two models with different parametrizations (i.e. going from a SIE with einstein_radius as its normalization to an NFW which has kappa_s as its normalization).

I want us to aim for an automated prior linking scheme, whereby the 'pass_priors' function becomes as 'overwrite_priors' function such that users have some sense of how priors have automatically linked. Its bold, and a little bit crazy, but I believe it can work.

The first issue we need to address, is how do we tie two profiles together with two different parametrizations? This can be performed by pairing them based on their integrated quantities, for example, their luminosity and mass within an aperture.

For example, lets say we have an SIE mass profile, and we compute its mass within 1" is 10 mass units. We can go to the NFW mass profile we're linking it to, and ask, 'what value of kappa_s gives us 10 mass units in an aperture of 1"'. If we have two parameters that control mass, we can do this on a 2D grid. If we have 5 parameters, we can set-off a quick downhill simplex or MultiNest search to find the combination of 5 parameters that give 10 mass units in 1". In Implementation, lets always just use a non-linear search like MultiNest, it'll be easier.

The parameters we infer thus link two mass models. They can be used as Gaussian priors in the next phase. The sigma's on these priors can use the priors/width config files.

For light profiles, we can play the exact same game with luminosity.

Geometric parameters should not be linked in the way above. Here, we can simply pair each profile's centre, axis ratio and phi values.

This pairing should be done on a galaxy-by-galaxy basis. Provided the order of the galaxies supplied to a phase does not change (i.e the list of galaxi_models in lens_galaxies=[] and souce_galaxies=[]) then automatic prior linking works fine.

In fact, this may motivate a framework where the number of galaxy's does not change in a pipeline, unless an add_galaxy function is called. It would be worth discussing whether the downsides of this would conflict with model_mapper.

If two galaxy models are linked, where the previous galaxy model had aligned parameters or a parameter set to a constant, we need to define the default behaviour. Does the constant parameter become free? Do the align parameters stay aligned? In the pass priors, the current behaviour is that constants stay constant, linked parameters stay linked, etc, so I guess we would adopt this as default. Maybe we could have a boolean flag that switches all parameters to variable?

We've got a lot to work on, but with you currently working a lot on model mapper this may be something that you can start on (albeit, it well could not). I can write the function that pairs two models based on their mass / luminsoity.

I'm now in the processing of profiling the code, and comparing run-speeds to FORTRAN. I'll note in this issue aspects of the code which are too slow, and how we can speed them up.

Convolution via frame convolver:

ORIGINAL:

psf = 21x21

lsst_solution: 0.7838027477264404
euclid_solution: 2.7468605041503906
hst_solution: 9.063872337341309
hst_up_solution: 20.911508083343506

psf = 41x41

lsst_solution: 3.622267246246338
euclid_solution: 10.063962459564209
hst_solution: 36.11542248725891
hst_up_solution: 121.993239402771

JITTED:

psf = 21x21

lsst_solution: 0.03985404968261719
euclid_solution: 0.09821200370788574
hst_solution: 0.30489659309387207
hst_up_solution: 0.7619879245758057

psf = 41x41

lsst_solution: 0.07627749443054199
euclid_solution: 0.1665787696838379
hst_solution: 0.45037412643432617 ( FORTRAN = 0.04)
hst_up_solution: 1.0230822563171387 (FORTRAN = 0.096)

NOTES:

• Numba is key to speeding convolution up to reasonable run-times.
• The frame convolver is x10 slower than FORTRAN - and is giving run-speeds that will be problematic for a realistic lens analysis.

The frame convolver should not be slower than FORTRAN after jitting - as the FORTRAN code includes the process of mapping the image to 2D for convolution and reducing it back to 1D.

The reason the frame convolver is slower is because of its large number of -1 entries, which represents pixel in the frame that are outside of the masked region. The code reads every -1 and then continues onto the next index. Below is the fraction of values for each frame array that are -1s:

LSST:

frame_array: 0.4806965585227189
buring_frame_array: 0.8761812348231439

Euclid

frame_array: 0.24736280724621756
buring_frame_array: 0.8446561570493754

HST:

frame_array: 0.12451287733922718
buring_frame_array: 0.826131601559918

HST_up:

frame_array: 0.07482219308446714
buring_frame_array: 0.8173658286983688

Therefore, we need a new data representation that removes all -1's, albeit we need to think carefully about memory.

I suggest that we remove all of the -1's and store a new list of arrays which represent, for every element in frame_array / blurring_frame_array, the index of the kernel they map too. Then, instead of iterating over the length of the kernel in the frame convolver, we would iterative over the length of each frame_array / blurring_frame_array and grab the kernel index using the new array.

Currently frame_array and blurring_frame_array are lists of numpy arrays. This new scheme would allow us to make them one giant NumPy array. Numba might be happier if we did this, but it could produce a hit on memory?

This problem will also effect the pixelization matrix convolution and be solved in the same manner!

For a HyperPhase, we need an unmasked image of every galaxy's model image and a masked image of every inversion's reconstructed model image(s). Extracting these images is not necessarily easy though, as the number of images that we need as well as where they 'come from' depends on the number of lens / source galaxies and whether a pixelization was used.

We should be able to easily put the logic necessary to extract these images in the fitting module, which has access to the tracer to determine the number of galaxies in each plane. For example, we already use the tracer to generate the unmasked image of each galaxy - I recommend we therefore make the 'hyper_galaxy_images' array be a list of planes of galaxies, e.g.:

hyper_galxy_images[0][0] - galaxy 0 in plane 0 (the image plane)
hyper_galxy_images[0][1] - galaxy 1 in plane 0 (the image plane)
hyper_galxy_images[0][2] - galaxy 2 in plane 0 (the image plane)
hyper_galxy_images[1][0] - galaxy 0 in plane 1 (the source plane - this could be the image of a profile or an inversion).

The HyperPhase can then just iterate over every image, optimizing the HyperGalaxy's associated parameters. It would be nice if the hyper_galaxies in the pipeline could follow the same structure (as opposed to being a 1D list) - lets think about this.

Extracting this array from the internal logic of the fitting classes should be straightforward..

The data structure above generalizes fairly well to multi-plane lensing.

We also need to optimize the pixelization / regularization parameters of every inversion. As we cannot have multiple inversions per plane, we can simply extract a list of inversions, where None corresponds to not inversion, e.g.:

hyper_inversioin[0] = None (No inversiion in image plane)
hyper_inversion[1] = Rectangular, RegConstant (source plane has pixelization and regularization)

Again, once we can extract these quantities running their optimization in a HyperPhase is straightforward.

Many of the objects used in the analysis are described mathematically as matrices.

A simplistic representation of a matrix in python is a list of lists:

[[]]

If this object is sparse then a lot of memory is allocated for empty entries. Furthermore, some operations require iteration over each element or each row/column of a matrix. If zero values do not contribute to the result of the operation then this time is wasted.

One way to represent the neighbour relation between source place pixels is with such a sparse matrix:

[[1, 1, 0, 0], 
 [1, 1, 1, 0],
 [0, 1, 1, 1], 
 [0, 0, 1, 1]]

In this representation a pixel address is a column or a row and a 1 indicates that two pixels are related. So pixel 0 has pixel 1 as a neighbour, pixel 1 has pixels 0 and 2 as a neighbour and so on.

This representation is very efficient for checking if two pixels are related. m[a, b] == 1 can be performed in constant time. However, if we wanted to perform some operation over all of the neighbours of pixel 0, for instance, we would have to check four entries in the matrix even though only two contain a non-zero value.

A more efficient representation is an list of lists of addresses:

[[1], [0, 2], [1, 3], [2]]

From sixteen elements we have reduced the size of the object to six. If we were to iterate over the neighbours of some pixel we could do that in O(N) time where N is the average number of neighbours.

If we wanted to check whether two pixels were neighbours that would now take O(N) time as well. We can recover a constant time lookup by using sets rather than lists:

[{1}, {0, 2}, {1, 3}, {2}]

These two objects still contain all the information that a simplistic matrix representation contained. We have all of the neighbour relationships and we can find the number of neighbours by checking the length of the list.

Another matrix representation of information in the analysis is the matrix that maps source pixels to image pixels:

[[0, 1, 0,    0], 
 [1, 0, 0.5, 0],
 [0, 0, 0.5, 0], 
 [0, 0, 0,    1]]

In this example, source pixel 0 maps to image pixel 1, source pixel 1 maps to image pixel 0 and so on. Source pixel 2 is an exception, as it maps half to source pixel 1 and half to source pixel 2.

This partial mapping means we can't use the representations used in the source-neighbour objects.

An appropriate representation could be a list of dictionaries:

[{1: 1}, {0: 1}, {1: 0.5, 2:0.5}, {3: 1}]

Here the same information is given but only non-zero entries need be accounted for.

Now convolution acts on an individual dictionary.

Let's imagine a simple 3x3 kernel:

[[0, 0.2, 0],
 [0.2, 0.4, 0.2],
 [0, 0.2, 0]]

We take the dictionary for source pixel 1:

{1: 1}

In image space this looks something like:

[[0, 1],
 [0, 0]]

Applying the kernel would give:

[[0.2, 0.4],
 [0, 0.2]]

However, in our efficient regime we can neglect image space entirely and simply perform a transformation on our dictionary representation:

{1: 1} -> {0: 0.2, 1: 0.4, 3: 0.2}

Our f matrix is now formed of a list of convolved dictionaries.

This representation makes computation of each element in the F matrix far more efficient.

Take F01. Ignoring noise, this is the sum over j of f[0, j] * f[1, j].

Using the classic matrix representation of f we would have to check the entries in f[0] and f[1] for each j. That's O(I) time complexity where I is the number of image pixels.

Given the list of dictionary representation, we only need to check for keys that actually exist in the dictionary.

Let d0 be f[0] and d1 be f[1]. We iterate through the keys in d1 and if a key exists in d2 we take the values, multiple them and add them to the sum.

sum = 0
for image_pixel_index in d0.keys():
    if image_pixel_index in d1:
        sum += d0[image_pixel_index] * d1[image_pixel_index]

Checking a key is in a dictionary is performed in constant time. In effect, the number of operations performed is equal to the number of image pixels that source pixel 1 maps to. This means that calculation of F matrix operations using this scheme is vastly more efficient.

Okay - if you have a spare moment, it'd be good to get your opinion on this before I push ahead.

The galaxy class is a bit weird, if you look at the code, you'll see that we setup each Galaxy, and then feed them to the RayTracingPlane. The RayTracingPlane class then updates a lot of the properties of the galaxies, as they have a lot of properties that require knowledge of the entire ray tracing plane (e.g. the redshifts of the galaxies before / after each galaxy) in order to calculate.

This gets quite nasty, as it means the Galaxy class has a lot of methods which won't work until they've been passed through the RayTracingPlane class.

I'm going to split the Galaxy class in two, i.e. the Galaxy class as we use it now, and a RayTracingGalaxy, which the RayTracingPlane uses to set up its galaxies. The RayTracingGalaxy will inherite from the Galaxy class, and include all the methods which require knowledge of the RayTracingPlane to perform.

There are two circumstances where we'd like to model multiple images of the same strong lens galaxy:

1. When we have multiple images of a lens taken during the same telescope observation (and therefore at the same wavelength). These images are most likely 'dithered', so that an under-sammpled PSF can be more accurately modeled. We may wish to model the individually dithered images (rather than their combined image) as this will reduce correlated noise (I need to grasp what trade off there is here with the PSF accuracy).

2. When we have multiple images taken at different wavelengths, potentially with completely different telescopes.

Multiple source-planes

The biggest change to the analysis when modeling multiple images is that one no-longer has one source-plane, but instead a single source-plane for each image. Regardless of whether the images share the same resolution, we can always reconstruct these multiple source's using the same pixelization (i.e. we derive the pixelization on one image and use it across the remaining images).

How do we represent data when there are multiple images

If we have one image, we'd currently mask it so that it is a list of image pixel indexes (lets say, 0-1000).

To represent a mapping from an image to source pixel (e.g. the f matrix), we use a list of size source pixels containing dictionaries of each image pixel mapping. So, if source pixel 2 mapped to image pixels 3, 5 and 8, with values of 0.5, 1.5 and 2.5:

mapping_matrix[2] = {3 : 0.5, 5 : 1.5, 8 : 2.5}

Likewise, for the F covariance matrix we have a dictionary, containing a tuple of every source pixel pair and the value of their covariance. So, if source pixels 4 and 8 were covariant and the value of this was 80.5:

F[(4, 8)] = 80.5

The first question that arises is how do we want to number multiple images and their source plane's?

We have two options:

1. All image pixels and source pixels continue their numbering from the previous image. So, if image 1 had image pixels 0-1000 and source pixels 0-100, image 2 would have image pixels 1001-2000 and source pixels 101-200. This does not fit in with how we coded up the analysis module thus far, which will often implicitly assume image / source pixels are ordered from 0.

2. We build 'super' structures for f and F, which handle the multiple images individually. So, following the mapping matrix example above (e.g. the same values), if that example were for image 7 we would build its mapping matrix and then store it:

multi_image_mapping_matrices[7][2] = {3 : 0.5, 5 : 1.5, 8 : 2.5}

F we can do the same:

multi_image_F[7][(4, 8)] = 80.5

Clearly, option 2 is the best - it reuses all our code and keeps image pixels / source pixels defined in a unified way that makes it easier to apply other functions. Our multi_image variables are now just lists of the structures of the mapping_matrix and F (and of course also D and H). There is one (obvious) downside to this option, detailed below.

Regularizing the source-plane across images

To reconstruct each image, we invert its F and D matrix to compute S. However, it is desirable that we reconstruct all source's simultaneously - that is, we compute S for all source-planes by fitting for all of the F and D matrices simultaneously. Say we've two images, each with an F matrix of size 200 x 200 and D matrix of size 200. To invert both simultaneously we'd build an F_all matrix of size 400 x 400 and D_all matrix of size 400. The F_all matrix would be formed by making its top-left corner F of the first matrix and bottom-right corner F_all of the second matrix. D_all would be a stack of the two individual D matrices. We could then compute S_all by inverting F_all and D_all exactly as we would an individual F and D matrix. However, the bottom-left and top-right corners of F_all are currently all zeros, which means there is no covariance between source-planes and that inverting F_all and D_all in this way gives an identical result to inverting F and D for each image individually.

it is when we introduce source-pixel covariances across the different source planes invert F_all and D_all is important, via the regularization matrix H_all. H_all is initially built the same way as F_all, by adding the individual H matrix of each source-planes to the top-left and bottom-right corners of H_all (Which are identical, as we are using an identical source-pixelization for each image). However, we also want to regularize the same source pixel in each source-pixelization (as they are fitting the same region of the image). This involves adding non-zero terms to the H_all matrix in the bottom-left and top-right corners (as, if the source pixels are numbered 1-400, we can regularize between say source pixel 5 and 205, by adding a non-zero entry here. We'd do this for every source pixel pair in each source plane eg. 1 and 201, 2 and 202, 3 and 203 etc. as well as 203 and 3, 202 and 2, 201 and 1 etc.). Simple, and we have our covariant source reconstruction :).

Herein lies the problem, the description above requires our image pixel indexing to go from 0->400, not be two lists going from 0->200. This doesn't effect our mapping_matrix, so the idea of
multi_image_mapping_matrices[7][2] = {3 : 0.5, 5 : 1.5, 8 : 2.5} is fine, its essentially multi_image_F[7][(4, 8)] = 80.5 we need a rethink, as this doesn't allow us to represent the covariances across different images (via their source planes).

The solution is that multi_image_F is not a list of dictionaries, but a dictionary of dictionaries, where the first dictionary stores the mapping between different F matrices. So, to represent the covariances in image 7 we would use multi_image_F[(7, 7)][(4, 8)] = 80.5, but if images 5 and 7 were covariant we'd have multi_image_F[(5, 5)][(1, 1)] = 1.0, multi_image_F[(5, 5)][(2, 2)] = 1.0, multi_image_F[(5, 5)][(3, 3)] = 1.0 etc. Note that in this framework, because source pixels are typically only covariant with their equivalent source pixel in another source plane, this makes all the off-diagonal F matrices identical and simple a line of values down the diagonal.

Excellent - now to just write the Cholesky Decomposition and inversion routine for all this (I'm sure beautiful symmetries will emerge)...

To speed up KMeans, we can track the positions of a set of image pixels and track how they move in terms of centres shifting, rotation, distance from the centre etc.

We can then move the initial source-pixel centers of the source-plane to align with the expected grid.

The grids module is now in a state where I intend to leave it, until Rich is back. This issue is to describe how I set it up and what needs to be done next. It is very possible it'll need a refactor - but its in a state where we can immediately figure this out when you're back!

In Scripts, check out the 'grid_demo.py' file to see how grids work with code ;).

Grids

The grids module essentially contains a set of structured NumPy arrays, which hold the (x,y) coordinates of an image. Coordinates *are held as either :

GridRegular - a NumPy array where the coordinates correspond to a uniform regular grid. This has shape ndarray[image_pixels, 2]

GridSub - a NumPy array where the coordinates correspond to the sub-grid of each pixel. This has shape ndarray[image_pixels, sub_pixels, 2]

Grids which correspond to regions of the data inherit from the above two classes:

GridImage - The regular grid of unmasked image pixels.

GridImageSub - The sub grid of unmasked image pixels.

GridBlurring - The regular grid of pixels which are outside the mask, but need to be included in the analysis due to PSF blurring effects.

Finally, there is a GridCollection, which stores all of the above grids and is the object passed to the ray-tracing module for lensing calculations.

I anticipate we will have other grids added in the future, for different types of pixelizations / calculations.

GridData and GridMapper

In addition to the grids above, which take care of coordinates, we'll also have:

• GridData - This grid follows the same structure as GridRegular e.g. it is of shape [image_pixels]. It does not contain coordinates, but data instead (e.g the image fluxes, noise-map, etc.). We'll ultimately be comparing the images we make on GridRegular to the actual data, thus its ideal that they share the same structure.

GridMapper - This maps pixel indexes between grids. Currently, the code has a GridMapperSparse, which maps each image pixel in GridImage to their nearest pixels defined on a sparse uniform grid (e.g. if grid_mapper.image_to_sparse[4] = 3, then the 5th image pixel maps to the 4th sparse pixel). This is used purely for optimization, which I suspect most other GridMappers will be.

Currently, we can't map any of our grids, whether they be GridData, GridRegular or GridSub back to their original 2d images. This is straight forward to add, and requires a GridMapper to be added to do this. I'll have probably added this by the time you read this!

Computing Deflections and Intensities

GridRegular has functions which:

1. Take in a list of galaxies.
2. Compute the deflections / intensities for all mass / light profiles in each galaxy.
3. Return the deflection angles / intensities in the same structure as the grid (noting the grid's last dimension is 2 as it is an (x,y) grid, so for a grid of computed intensities this dimension is dropped).

GridSub has functions which do the same, but:

• If it is computing deflection angles, the complete set of deflection angles on every sub-pixel are returned.
• If it is computing intensities, each set of sub-pixel intensities are averaged to compute the overall intensity in the pixel.

These choices reflect the fact that a sub-grid of deflection angles is required for setting up the next grid of sub-pixels and performing the source reconstruction, whereas there is no merit in keeping track of sub-pixel intensities as only the pixel intensities can be compared to the actual image data.

Deflection angles can be computed simultaneously across all grids at the GridCollection level.

Ray Tracing and Grids

The ray_tracing module receives an instance of GridCollection as its input, as well as a list of galaxies. It groups the galaxies into different planes in the ray-tracing calculation:

• ImagePlane - The first plane in the ray-tracing calculation, corresponding to the lowest redshift galaxy(s). Note that because this is the image-plane, its coordinates correspond to the regular uniform grid of the CCD imaging data.

• LensPlane - A plane in the ray-tracing calculation where the grid coordinates have been lensed by a foreground galaxy(s) and are thus not on a uniform grid. This is not the last plane, so the deflection angles at the plane's coordinates are computed such that they can be traced to the next plane.

• SourcePlane - The last plane in the ray-tracing calculation. Because it is the last plane, no deflection angles are computed as there is not another plane to trace the coordinates to.

To set up each plane, the deflection angles are computed in each plane using the GridCollection deflection angle functions. Each Plane therefore has two _GridCollection_s - a collection of (traced) image coordinate and collection of deflection angles at those coordinates.

To keep things simple, the ray_tracing module currently assumes just two planes, an image-plane and source-plane. In the future, we'll want to implement multi-plane lensing, where the distribution of galaxy redshifts dictates the number of planes, and the angular diameter distances separating these planes determines the deflection angles.

I have made sure the way I've designed this code won't hinder us setting up multi-plane lensing. My suspicion is it will make it so that we store the planes as a list, as opposed to explicitly calling them image_plane and source_plane as we do now.

Generating Images from a Plane

Once ray-tracing is compute, and we therefore have an instance of the RayTrace class, it is now possible to generate what the images of each plane of galaxies looks like. A given plane is comprised of: (i) grids of coordinates and; (ii) galaxies with light profiles. Generating its image requires us to evaluate all of the plane's galaxy light profiles at its grids of coordinates, which is performed by passing the galaxies into our Grids.

Computing the deflection angles on each grid (e.g. image, sub, blurring) is straightforward. For each grid, whether it be a GridImage, GridImageSub, or GridBlurring, all you do is compute the deflection angles due to each plane's galaxies at the coordinates defined by the grid.

For generating a galaxy image, via its light profile and intensities, things are not so straight-forward. This is because there are multiple ways that one might want to compute the intensities from our GridCollection:

• Using GridRegular, thus ignoring sub-pixel gridding.

• Using GridSub, such that each intensity is the average of the sub-pixel intensities over that pixel.

In light of this, I have added to the light_profiles module a class, LightProfileSettings. This allows one to customize, for each light profile instance, how its intensities are computed from an input gird. Thus, we could have a galaxy where the first light profile is computed using GridRegular and second using GridSub (this isn't an unrealistic occurrence).

Here is where it gets tricky, as when we go to each plane in the ray-tracing instance, we need to use its galaxies to go to their light profiles and correctly choose each grid which it uses. The code currently uses intensity calculations in the Galaxy class, thus this may require us to filter light profiles.

Whilst both methods of computing intensities above are implemented in the code, the code does not choose its method based on the LightProfileSettings. I figured that was the point where refactoring a bad solution could take time, so stopped there. Currently, the code defaults to always use GridRegular.

Iterative Sub-grid

There is a third approach to computing intensities (which is also controlled by LightProfileSettings). Here, GridRegular is used in-conjunction with and an iterative sub-grid, whereby intensity values are computed on iteratively higher resolution sub-grids until the change in value reaches a threshold value of precision. In the centre of a light profile, sub-grids of order (1000 x 1000) may be necessary, which are clearly too high resolution to be set up using GridSub.

In the geometry_profiles module, you've already coded up something to this effect, albeit not using the GridRegular data structure.

This approach is only valid when the coordinates are on a regular uniform grid, which is true only when they are in the ImagePlane. Thus, it should be turned off for any other grid, regardless of LightProfileSettings.

There are also circumstances where we may wish to evaluate a light profile in the image-plane, even though it traces to a higher redshift LensPlane or SourcePlane. In this case, we will still assign the galaxy to the higher redshift Plane, but will evaluate its light profile using the ImagePlane grid. In LightProfileSettings, I have added a variable 'image_plane_override' to control this behaviour. It is not yeet implemented in ray_tracing, but I'd guess the ImagePlane class could have an extra input in addition to galaxies called something like 'extra_light_profiles'.

MappingMatrix

Above, we're using light profiles to generate images. If a Plane has a Pixelization, we will instead use the plane coordinates to generate a MappingMatrix. This is the matrix we've discussed and used in the past, which maps every image pixel to source pixel using 0s and 1s (or fractional values for a sub grid). This isn't implemented yet, as the pixelization module has some work to do, but it should be fairly straightforward to implement as it doesn't have the confusing issues light profiles do above.

Analysis

So, once this is all set up and working, the point is we can use the ray_tracing module to:

• Generate model images of all galaxies, which we can analyse by comparing to the image data (e.g. blurring with the PSF, subtracting from the image data to compute residuals, evaluating a goodness of fit).

• Generate the mapping matrices of our pixelizations, which we can blur with the PSF and fit the image data with.

This will take place in some form of 'analysis' module, which I'll start writing on Friday :).

Everything is well tested and behaves as expected, but its not verbose and a bit of a mess to read. A lot of the deflection angle / potential routines come from a weird hierarchy too.

The problem is basically that, in order to exploit faster computations via symmetries in spherical / simpler profiles, we have to over-ride a lot of methods.

I don't think there's necessarily a simpler way to write a lot of these routines, so can we make things more clear by splitting mass_profile into different files and having clear documentation for each routine? For each profile, we should write explicitly where it is inherited routine from and add latex documentation with the equations used.

Once we have computed a model image of a galaxy and its blurring region, this region performs their PSF convolution by:

• converting the image to 2D
• converted the blurring region to 2D
• Performing the PSF convolution in 2D
• Converting the blurred image back to 1D

We should do this in 1D, by using an adapted version of the frame convolver. This will pretty work exactly the same as it does for the source-plane, but must also include the blurring region pixels surrounding the mask. Provided its set up correctly, this will be as efficient as possible.

We may want to use different sized PSF's for this convolution and that performed in the source plane. Therefore, we probably want to set this frame convolver up independently from the source-plane one. I'm a bit wary of how much memory this could take, but probably best not to worry for now.

The results of the unit tests will not need changing once we implement this.

When executing the line

profile_plotters.plot_intensities(light_profile=sersic_light_profile, grid=image_grids.image)

In the "howtolens" tutorial "2_profiles.ipynb" I receive the following error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
\<ipython-input-6-27d07a0d0848\> in \<module\>()
      1 # We can use a profile plotter to plot this intensity map (the image is mapped to 2D before plotting).
----> 2 profile_plotters.plot_intensities(light_profile=sersic_light_profile, grid=image_grids.image)

~/anaconda3/lib/python3.7/site-packages/autolens/plotting/profile_plotters.py in plot_intensities(light_profile, grid, output_path, output_filename, output_format)
     12         array=intensities, points=None, grid=None, as_subplot=False,
     13         units='arcsec', kpc_per_arcsec=None,
---> 14         xticks=grid.xticks, yticks=grid.yticks, xyticksize=16,
     15         norm='linear', norm_min=None, norm_max=None, linthresh=0.05, linscale=0.01,
     16         figsize=None, aspect='auto', cmap='jet', cb_ticksize=16,

AttributeError: 'ImageGrid' object has no attribute 'xticks'

The model.info files now double output the source and lens galaxies in the priors, e.g.:

EllipticalSersic

lens_centre_0 GaussianPrior, mean = 0.0, sigma = 1.0
lens_centre_1 GaussianPrior, mean = 0.0, sigma = 1.0
lens_axis_ratio UniformPrior, lower_limit = 0.2, upper_limit = 1.0
lens_phi UniformPrior, lower_limit = 0.0, upper_limit = 180.0
lens_intensity UniformPrior, lower_limit = 0.0, upper_limit = 1.0
lens_effective_radius UniformPrior, lower_limit = 0.0, upper_limit = 4.0
lens_sersic_index UniformPrior, lower_limit = 0.8, upper_limit = 8.0

EllipticalSersic

lens_galaxies_centre_0 GaussianPrior, mean = 0.0, sigma = 1.0
lens_galaxies_centre_1 GaussianPrior, mean = 0.0, sigma = 1.0
lens_galaxies_axis_ratio UniformPrior, lower_limit = 0.2, upper_limit = 1.0
lens_galaxies_phi UniformPrior, lower_limit = 0.0, upper_limit = 180.0
lens_galaxies_intensity UniformPrior, lower_limit = 0.0, upper_limit = 1.0
lens_galaxies_effective_radius UniformPrior, lower_limit = 0.0, upper_limit = 4.0
lens_galaxies_sersic_index UniformPrior, lower_limit = 0.8, upper_limit = 8.0

These are the same model, but being output twice.

I have fixed this by preventing 'lens_galaxies' and 'source_galaxies' from being output:

` @Property
def info(self):
"""
Use the priors that make up the model_mapper to generate information on each parameter of the overall model.

    This information is extracted from each priors *model_info* property.
    """
    info = []

    for prior_model_name, prior_model in self.flat_prior_model_tuples:

        # TODO : clean this up

        if 'lens_galaxies' not in prior_model_name and 'source_galaxies' not in prior_model_name:

            info.append(prior_model.cls.__name__ + '\n')

            prior_model_iterator = prior_model.prior_tuples + prior_model.constant_tuples

            for i, attribute_tuple in enumerate(prior_model_iterator):
                attribute = attribute_tuple[1]

                line = prior_model_name + '_' + attribute_tuple.name
                info.append(line + ' ' * (60 - len(line)) + attribute.info)

            info.append('')

    return '\n'.join(info)`

This now produces the following output in model.info

EllipticalSersic

lens_centre_0 GaussianPrior, mean = 0.0, sigma = 1.0
lens_centre_1 GaussianPrior, mean = 0.0, sigma = 1.0
lens_axis_ratio UniformPrior, lower_limit = 0.2, upper_limit = 1.0
lens_phi UniformPrior, lower_limit = 0.0, upper_limit = 180.0
lens_intensity UniformPrior, lower_limit = 0.0, upper_limit = 1.0
lens_effective_radius UniformPrior, lower_limit = 0.0, upper_limit = 4.0
lens_sersic_index UniformPrior, lower_limit = 0.8, upper_limit = 8.0

This is good - previousouly this looked like

EllipticalSersic

lens_galaxies_lens_centre_0 GaussianPrior, mean = 0.0, sigma = 1.0
lens_galaxies_lens_centre_1 GaussianPrior, mean = 0.0, sigma = 1.0
lens_galaxies_lens_axis_ratio UniformPrior, lower_limit = 0.2, upper_limit = 1.0
lens_galaxies_lens_phi UniformPrior, lower_limit = 0.0, upper_limit = 180.0
lens_galaxies_lens_intensity UniformPrior, lower_limit = 0.0, upper_limit = 1.0
lens_galaxies_lens_effective_radius UniformPrior, lower_limit = 0.0, upper_limit = 4.0
lens_galaxies_lens_sersic_index UniformPrior, lower_limit = 0.8, upper_limit = 8.0

The issue is because im guessing we can do a better fix than an if loop

These two functions don't pass the non-grid attrbutes (e.g. the mask) to new ImagingGrids setup during ray_tracing. These attributes must be passed for visualization.

In ray tracing... you can see I've made a prertty nasty hack to fix this for now. If you could fix these two functions then we'll have neat tidy code again!

Nan has a GUI lens modeling interface, which allows one to initialize a lens model using a mouse. This can be seen using the link below, and works brilliantly :).

(Nan attatch link below this post)

The plan is to attach this GUI code and interface to AutoLens, such that in AutoLens one can:

1. Load an image from a fits file for lens analysis.
2. Bring up that image in the GUI.
3. Setup the initial lens model parameters using a mouse, just like one can following the link above.
4. Begin an AutoLens analysis using this initialized mass model.

In terms of code, the GUI interface is a separate module, which uses only the AutoLens Profile module (e.g the mass profiles for deflection calculations, light profiles for galaxy light, etc.). The GUI interface does not need to feed anything 'back' to AutoLens. Rather, it will output a new config file (e.g. GUI.ini), such that the AutoLens analysis priors can be initialized using the GUI interface's SIE mass profile.

To begin implementing this, Nan needs PyAutoLens to be able to do a non-linear optimization fit of a lens model, specifically a mass (SIE) + parametric source (Gaussian) model. So, I need to get the ray_tracing module capable of generating model images, and the bones of an analysis module together for computing chi^2. We may also require a different non-linear optimizer than MultiNest.

Long-term functionality of GUI Interface

1. AutoLens uses 'positional image pixels' to speed up the analysis, by discarding mass models where these pixels do not trace within a threshold value to one another. Thus, these pixels should correspond to the same region of the source in the source plane, or what 'looks like' the multiple images of the image in the image-plane.

It'd be cool if the user could draw these pixels on in the GUI interface, by highlighting a circle around each image and locating the brightest pixel, or something.

2. Currently, the GUI loads an image which includes the lens's light, and has no means by which to subtract it. Obviously, for low Einstein radii lenses, this means the source is not visible. Thus, we should ruin a non-linear light profile fit to the image beforoe the GUI runs, to get an initial set of Sersic parameters, and the user should be able to adjust these so as to better reveal the source light. This therefore also initilize's the lens light profile, which is neat.

As discussed in #11 , we'll have bespoke analysis routines for computing light profiles / mass profiles during the lens analysis.

After an analysis is run, we want a lot of flexibility in how one can visualize the results For example, after we have a lens model, we may wish to plot the light profile of an individual component of a galaxy (e.g. LightProfile), all the light profiles of a galaxy (e.g. LightProfileCollection or Galaxy) or all the galaxies in the model (e.g. GalaxyCollection). The same is true of mass profiles.

Therefore, I suggest we include in all these classes plot routines which use an array of data to make the plot (you'd already started working on this, so I'll use the code that already exists). This code is different to that in the analysis routine, but we can tie them together via unit tests. To achieve this we need to do the following:

1. For light profiles, one should be able to retrieve either the light profile or the PSF convolved light profile.

However, the reality is the routines above won't handle the visualization of results after an analysis that is particularly computed, for the reasons below:

1. For a pixelized source reconstructions, this must run on the light-profile subtracted image, or the linear inversion will fit the lens galaxy light and look stupid. If there are multiple source-planes, the sources must be reconstructed simultaneously, information which again cannot be built into the Pixelization class.

2. When there are multiple galaxies, the light profile is computed not on the image coordinates but the coordinates already deflected by foreground galaxies. Equally, things like mass profiles should reflect the deflection angles of foreground deflections. Neither of these effects are reflected in the LightProfile or MassProfile class.

3. Some things we want to visualization are inherently tied to the analysis itself (e.g. the image, chi^2 map, residuals etc.). This includes a dependence on hyper-parameters.

Therefore, I suggest that we'll need a Visualization class, which takes the results of an analysis and computes all the quantities one could ever want to plot and stores them in memory, e.g.:

1. The light profile, surface density, deflection angles, potential etc. for every Galaxy, allowing one to plot what they look like after foreground lensing and in their true source plane.

2. All of the source reconstructions, as they look after the final analysis so including the complications above.

3. Everything else - the residuals, chi sq image, hyper-parameter associated images and so forth.

We could make a lot of flexible routines to allow one to do weird things that arn't necessarily obvious (e.g. subtracting specific components of a galaxy light profile from an image).

Input welcome, but I don't see a neater way to do this. It'd include all the parameter degeneracies images too.

So, i was a bit stupid - this is a fairly trivial task. We have the list of neighbor pairs, e.g.:

pairs = [[0,1],[0,2],[1,4],[4,6]] etc.

We just have to go to each pair and put them in a neighbors vector i.e.:

for pair in pairs:
neighbors[pair[0]] = pair[1]
neighbors[pair[1]] = pair[0]

Im struggling to set up the list so it adds points... or do this as a numpy array. But its late, but if you know a quick code then let us know!

We have an efficient data structure for our F matrix, where each value is stored in a dictionary. We can do the same for the H matrix trivially. Unfortunately, we currently need to convert this to a matrix to compute the determinant and solve the set of linear equations of these matrices.

Thus, we now want to use our data structure to perform the above tasks. We need the determinant of two matrices (H and F+H) and need to solve for one matrix (S, from D and F+H). For both tasks, we need to compute the Cholesky decomposition of our matrices, e.g. F+H = LL* and H = LL*, which represents them as the multiplication of a diagonal matrix (i.e. all zeros in the upper or lower diagonal) multiplied by its conjugate transpose. There are standard algorithms for doing this.

The determinant is now easy to compute, it is simply the sum of the natural log of all diagonal terms of L, multiplied by 2:

Det_F = 2.0* sum_i log (L(i,i))

The inversion requires more care, using forward and back substitution . I haven't yet wrapped my head around how exactly that works, but the c# code supplied in the link is only 30 lines long so it can't be too tricky (right?...).

scipy.spatial.Voronoi does not return a list giving the points that a voronoi vertex is adjacent with. i.e., if a Voronoi cell had 5 neighbours, whos points corresponded to the entries [1, 6, 10, 25, 39] in an input list of points, the list [1, 6, 10, 25, 39] is not returned from the algorithm.

It can be calculated from the results of the algorithm using the ridge_points array, which gives a list of every neighbouring pair of points (i.e. if points 3 and 7 share a Voronoi vertex, and points 4 and 12 do, Voronoi.ridge_points will have entries of [3,7] and [4,12]). Making the desired array could be slow though, as one would need to apply a lot of filtering w/ logic.

Alternatively, I think we can compute what we need using a Delaunay triangulation (scipy.spatial.Delaunay), which is probably pretty quick (given that Voronoi is so fast anyway). Still, its a shame to waste time and write code on something the voronoi algorithm should be able to compute.

The other option is we find another Voronoi algorithm.

We can use pass priors to set parameters to constants, e.g.:

 def pass_priors(self, previous_results):
         self.lens_galaxies[0].sersic.axis_ratio = 0.2
         self.lens_galaxies[0].sersic.phi = 90.0

    phase1 = MMPhase(lens_galaxies=[gp.GalaxyModel(sersic=lp.EllipticalSersic)],
                               optimizer_class=nl.MultiNest, phase_name="{}/phase1".format(pipeline_name))

We then pass priors to make parameters variable in the next phase:

def pass_priors(self, previous_results):
    self.lens_galaxies = previous_results[0].variable.lens_galaxies`

However, the constants are not being linked consistently. In the example above, phi remains a constant whereas axis_ratio becomes variable. Checkout the 'link_variable_with_constant_floats.py integration test.

There are also issues when linking tuples, e.g.:

  class LensPlaneGalaxy0Phase(ph.LensPlanePhase):
        def pass_priors(self, previous_results):
            self.lens_galaxies[0].elliptical_sersic.centre_0 = -1.0
            self.lens_galaxies[0].elliptical_sersic.centre_1 = -1.0

If in the next phase we write:

    class LensPlaneBothGalaxyPhase(ph.LensPlanePhase):
        def pass_priors(self, previous_results):
            self.lens_galaxies[0] = previous_results[0].variable.lens_galaxies[0]

This cases an error as the code doesnt know what to do with the centres anymore. We have to manually link the constants:

class LensPlaneBothGalaxyPhase(ph.LensPlanePhase):
    def pass_priors(self, previous_results):
           self.lens_galaxies[0] = previous_results[0].variable.lens_galaxies[0]
           self.lens_galaxies[1] = previous_results[1].variable.lens_galaxies[0]
           self.lens_galaxies[0].elliptical_sersic.centre_0 = -1.0
           self.lens_galaxies[0].elliptical_sersic.centre_1 = -1.0

Checkout the 'link_variable_with_constant_tuples.py' integration test.

The desired behaviour is that constants remain constant and variables remain variable.

It may be nice to have the options. For example, if we put:

   class LensPlaneBothGalaxyPhase(ph.LensPlanePhase):
        def pass_priors(self, previous_results):
            self.lens_galaxies[0] = previous_results[0].variable.lens_galaxies[0]

constants remain constant. However, if we put:

    class LensPlaneBothGalaxyPhase(ph.LensPlanePhase):
        def pass_priors(self, previous_results):
            self.lens_galaxies[0] = previous_results[0].all_variable.lens_galaxies[0]

The constants also become variable.

To adapt to the source's morphology, we need to cluster source pixels in its brightest regions but weighting the clustering to its brighter pixels. In Fortran, this used a weighted k-means algorithm, which Python doesn't appear to support.

We either need to find another algorithm, or come up with an approach that addresses this by removing faint pixels from the vector of points fed into the KMeans, thus making similiar behaviour. This approach is appealing in that it will greatly speed up the KMeans clustering, if it works as it should.

I'm going to change things slightly from fortran AutoLens, which tried hard to keep noise terms associated with standard data reduction procedures. In hindsight, this wasn't necessary and convoluted the code design.

For an image, one can easily define a 'baseline variance map'. This is essentially the variance in each pixel after accounting for all known noise sources (e.g. Poisson counting statistics from the galaxies / background sky, and read noise). If our image is background sky subtracted and in electrons per second, each variance (in EpS) is therefore:

1. sqrt( (sigma_BG + d )* exptime) / exptime

Where exptime could be a 2D image of the effective exposure time of each pixel.

The scaled variance map used in the Fortran AutoLens would add variances generating separate Poisson and Background variance maps - which represented separately the variances due to the galaxy counts and background. This is what I'm going to do away with, instead only generating a ScaledNoiseMap, which increases the variances relative to the BaselineNoiseMap in one of three ways:

1. It increases the in the lens galaxy, by generating a lens galaxy contribution map (see AutoLens paper) and increasing the variances. In Fortran AutoLens, this would explicitly convert each pixel from eps to counts, compute the Poison variance of each pixel by sqrting its value and then convert back to eps. I'm suggesting instead that we now simply apaply the contribution map to the baseline variances - yes these include the background but it'll make the code simpler. Ultimately, the hyper-parameters should just correct for this.

2. The same process in the source, using the source contribution map. My suggestion is the same as above.

3. The background noise map. I simply suggest the code can increase all background variances by a constant multiplication term, using the effective exposure time map.

Background Sky Subtraction

This can also now simply be included in AutoLens as a term which increases all image pixel values by a multiplication term. This is because the image is in electrons per second. The benefit of doing away with PoisonNoise and BackgroundNoise terms / classes is that this doesn't need to recompute each for a new sky background subtraction.

Weight Maps

Instead of using a variance map in the analysis, it may be better to use a weight map. The weight of each pixel is simply 1.0(/variance**2), thus it naturally enters the F and D matrices in the analysis. The benefit of this is that its more natural to 'mask' out a part of the image by setting its weight to 0, rather than its variance to infinity. I need to check how this interfaces with other parts of the code.

Okay, so the problem with paramnames is as follows:

A paramnames file is composed of two things, the paramnames_names and paramnames_labels.

• The paramnames_names are determined from the class constructor via model_mapper, thus when we change the number of parameters in a priors_pass the code knows that the number of paramnames_names has changed and adapts accordingly.
• The paramnames_labels are instead drawn from a property / list in the class. This list does not change / reduce its size when parameters become fixed by a priors pass.

Thus, we end up in a situation where the paramnames_names and paramnames_labels do not match, and often lead to different index sizes causing the error in the integration test.

I've split the create_paramnames function up a bit to make it easier to read / test, however to properly fix the issues I think we need some sort of class_priors_labels_dict in model mapper, which would be used in the same way paramnames_names uses class_priors_dict.

So, the tests are still all commented out for now but at least we know the issue now.

The interface will have the following three Phase's for the user to choose from:

LensPlanePhase - Only a lens plane, with lens_galaxies, is included - this is used purely for fitting the foreground objects with a light profile (i.e. no mass-profiles / deflection angles).

LensSourcePlanePhase - A lens and source plane phase, thus including lens and source galaxies. This is used for 99% of use cases where everything is treated on just two planes. The source-galaxy can have a pixelization or use light profiles.

MultiPlanePhase - An arbritrary number of planes, where each is composed on galaxies with light, mass and / or pixelizations. This phase will be integrated into AutoLens in a few months, once we're sure the above two are optimal.

However, the phase module itself will have many child phases, which are sub-sets of the above phases but specify whether the lens and source-plane use light-profiles or pixelizations to model the light. They also include hyper-parameters optimization phases and specify what is visualized. By my count, we need the following child phases:

• LensPhase - As LensPlanePhase above.

• LensHyperPhase - Use the lens galaxy light profiles fitted in LensPlanePhase to setup the HyperGalaxy of each galaxy.

• LensMassOnlySourceProfilePhase - The lens galaxy has no light profile (mass-profile only). The source uses light profiles.

• LensMassOnlySourceProfileHyperPhase - Use the source-galaxy light profile above to setup a HyperGalaxy.

• LensMassOnlySourcePixelizationPhase - The lens galaxy has no light profile (mass-profile only). The source uses a pixelization.

• LensMassOnlySourcePixelizationHyperPhase - Use the source-galaxy pixelization above to setup a HyperGalaxy and the pixelization's adaptive regularization / tesselation hyper-parameters.

• LensProfileSourceProfilePhase - The lens galaxy has light and mass profiles. The source uses light profiles.

• LensProfileSourceProfileHyperPhase - Use the source-galaxy light profile above to setup a HyperGalaxy.

• LensProfileSourcePixelizationPhase - The lens galaxy has light and mass profiles. The source uses a pixelization.

• LensProfileSourcePixelizationHyperPhase - Use the source-galaxy pixelization above to setup a HyperGalaxy and the pixelization's adaptive regularization / tesselation hyper-parameters.

Lets worry about child-phases of MultiPlanePhase later, once we have a clear idea of how the phase module looks above. Multi-planes will also care about whether profiles and pixelization are used.

Currently, the paramnames use lens_galaxies or source_galaxies, but without an index, e.g:

lens_galaxies_centre_0 x_{\mathrm{geometry_profile11}}
lens_galaxies_centre_1 y_{\mathrm{geometry_profile11}}
lens_galaxies_axis_ratio q_{\mathrm{geometry_profile11}}
lens_galaxies_phi phi_{\mathrm{geometry_profile11}}
lens_galaxies_intensity I_{\mathrm{geometry_profile11}}
lens_galaxies_effective_radius R_{\mathrm{geometry_profile11}}
lens_galaxies_sersic_index n_{\mathrm{geometry_profile11}}

The problem is that if we have two lens galaxies with repeated entries, we get duplicate paramnames, e.g.:

lens_galaxies_centre_0 x_{\mathrm{geometry_profile11}}
lens_galaxies_centre_1 y_{\mathrm{geometry_profile11}}
lens_galaxies_axis_ratio q_{\mathrm{geometry_profile11}}
lens_galaxies_phi phi_{\mathrm{geometry_profile11}}
lens_galaxies_intensity I_{\mathrm{geometry_profile11}}
lens_galaxies_effective_radius R_{\mathrm{geometry_profile11}}
lens_galaxies_sersic_index n_{\mathrm{geometry_profile11}}
lens_galaxies_centre_0 x_{\mathrm{geometry_profile11}}
lens_galaxies_centre_1 y_{\mathrm{geometry_profile11}}
lens_galaxies_axis_ratio q_{\mathrm{geometry_profile11}}
lens_galaxies_phi phi_{\mathrm{geometry_profile11}}
lens_galaxies_intensity I_{\mathrm{geometry_profile11}}
lens_galaxies_effective_radius R_{\mathrm{geometry_profile11}}
lens_galaxies_sersic_index n_{\mathrm{geometry_profile11}}

This causes problems with GetDist and leads to an error. So lets add an index to the paramnames:

lens_galaxies_0_centre_0 x_{\mathrm{geometry_profile11}}
lens_galaxies_0_centre_1 y_{\mathrm{geometry_profile11}}
lens_galaxies_0_axis_ratio q_{\mathrm{geometry_profile11}}
lens_galaxies_0_phi phi_{\mathrm{geometry_profile11}}
lens_galaxies_0_intensity I_{\mathrm{geometry_profile11}}
lens_galaxies_0_effective_radius R_{\mathrm{geometry_profile11}}
lens_galaxies_0_sersic_index n_{\mathrm{geometry_profile11}}
lens_galaxies_1_centre_0 x_{\mathrm{geometry_profile11}}
lens_galaxies_1_centre_1 y_{\mathrm{geometry_profile11}}
lens_galaxies_1_axis_ratio q_{\mathrm{geometry_profile11}}
lens_galaxies_1_phi phi_{\mathrm{geometry_profile11}}
lens_galaxies_1_intensity I_{\mathrm{geometry_profile11}}
lens_galaxies_1_effective_radius R_{\mathrm{geometry_profile11}}
lens_galaxies_1_sersic_index n_{\mathrm{geometry_profile11}}

Currently, to compute an intensity / surface density / potential / deflection angle value at a set of coordinates, we can only pass an individual set of coordinates into the code, e.g.:

sersic = light_profiles.EllipticalSersic(intensity=1.0, effective_radius=0.6, sersic_index=4.0)
intensity = sersic.intensity_at_coordinates(np.array([1.0, 0.0]))

We are also restricted to individual coordinates if we do this via a Galaxy class:

sersic = light_profiles.EllipticalSersic(intensity=1.0, effective_radius=0.6, sersic_index=4.0)
galaxy = galaxy.Galaxy(redshift=0.5, light_profiles=[sersic])
intensity = galaxy.intensity_at_coordinates(np.array([1.0, 0.0]))

However, in many calculations, we don't pass one set of coordinates to a function but a 'grid', which is a NumPy array with some defined structure. For example, an image.grid is a NumPy array with dimensions ( image_grid[image_pixels, 2] ), such that each entry is a coordinate on the image grid.

It'd be good if our _at_coordinates routines could take in an input numpy array of any shape, and output the results in a NumPy array following the same structure.

Currently, I have made a half-baked implementation of this by making separate routines for computing the values on a grid:

    def deflection_angles_grid(self, coordinates):

        deflection_angles = np.zeros(coordinates.shape)

        for defl_count, coordinate in enumerate(coordinates):
            deflection_angles[defl_count, :] = self.deflection_angles_at_coordinates(coordinates=coordinate)

        return deflection_angles

and on a sub grid:

    def deflection_angles_sub_grid(self, sub_coordinates):

        deflection_angles = np.zeros(sub_coordinates.shape)

        for sub_count, image_pixel in enumerate(sub_coordinates):
            for defl_count, sub_coordinate in enumerate(image_pixel):
                deflection_angles[sub_count, defl_count, :] = self.deflection_angles_at_coordinates(coordinates= sub_coordinate)

        return deflection_angles

It'd be good if we didn't need to set up separate functions for this. We could either make the functions generic for any NumPy array, or specific to the structures that our grids module uses (which at the moment is just the two above, image_grid[image_pixels, 2] and sub_grid[image_pixels, sub_grid_pixels, 2].

@rhayes777 When you're back, can your coding trickery easily set this up?

NOTE: For the intensity / surface_density / potential calculations, the input grid is in (x,y) coordinates, so always has size 2 in its last dimensions 2, e.g. image_grid[image_pixels, 2]. However, the output grid is not an (x,y) grid, so the dimensions output_grid[image_pixels]. In contrast, because deflection angles are (x,y) coordinates, the output grid is output_grid[image_pixels, 2]

I have written a very basic simulate module, which takes an image and adds PSF blurring / noise. However, the simulate module will be a key part of AutoLens, and will requite a lot more effort!

For simulating images, basically two things are needed:

• The code to simulate instrumental effects (e.g. read noise, poisson noise, PSF, background sky, etc.)
• An interface that allows the user to control these effects whilst not being bogged down in hard to understanding parameters.

On point 1 - The majority of effects we might want to simulate already have existing functions from libraries like AstroPy. In truth, other than the PSF, background noise and Poisson noise we probably won't want to worry about another effect for a while. I guess correlated noise due to image dithering would be the next, which goes hand-in-hand with multi image analysis and thus should be postponed until then.

On point 2 - To produce an image with appropriate signal-to-noise, I had to both change the background level of noise and the exposure time. Sometimes, the source brightness, lens light brightness and lensing magnification will also be parameters one is changing. For a user messing with parameters it becomes confusing whether the images they are producing are sensible or realistic.

When generating simulated images, users will want to interact with the code and generate images follwing two very different approaches:

• As physical as possible, using realistic magnitudes for the lens / source galaxies and realistic mass models, as well as physical parameters for the instrument. Thus, we probably want some classes which default to observed data-set results and instruments.

• Exploratory - e.g. where you ask questions like 'what does a broader PSF or lower S/N do to my results'. Here, you just want the images to come out 'looking right' but in a way you can change whatever you're interested in. For the problem above, we'd have want an automated routine that produces an image to a specified S/N.

Furthermore, it worth remember that we will want to tie the simulate module directly to the analysis modules, so one can swiftly simulate images and then see how fits behave.

The way we currently simulate images is a bit dodgy - you can't access stuff like the noise properly and don't have a great deal of customization. Furthermore, the generation of the image from a tracer is handled separately.

I believe we can make a much better simulation module by structuring it in the same way that Fitter is now structured. Specifically, we can perform all processing (in the correct order) in the class constructor (e.g. simulate image_plane_image -> blur with psf -> add sky -> add noise) and then access the different simulated quantities via properties which are automatically mapped to 2D. This will allow us to turn different features on and off without need copious amounts of if loops.

For the F matrix, the basic problem is we go to a row of our f mapping matrix, say row 0 for source pixel 0 and need to know which image pixels it maps to (after sub-gridding and PSF convolution).

However, we know from the clustering algorithm the set of image-pixels that source pixel 0 maps to. Thus, we can immediately find the image-pixels it directly maps to. The question is how do we efficiently locate the remaining image pixels (which it maps to only because of sub-gridding and PSF convolution).

Well, what we can do is we can go to the Voronoi neighbours of source pixel 0 and look for all non-zero image pixels. If we don't get them all, we can go to the neighbours of neighbours until we do. Image pixels will only correlate with neighbouring image pixels, which by definition are allocated to neighbouring source pixels!

The key questions for this scheme is when do we stop looking?

I think it is true that once you hit a source-pixel neighbour that maps to an image pixel whose value in its f matrix is 0 (i.e. they are not correlated due to subgridding or PSF convolution), you no longer need to look over its neighbours. its hard to imagine the next pixels outwards won't be all 0's (anything that isn't 0 is probably tiny anyway).

If the above isn't true, we could have a counting during or after the PSF convolution step, which tells us how many non-zero elements are in f for each source pixel. This could be computationally expensive.

One could implement a numerically exact method and compare to ensure accuracy.

Okay, so I've just spent the past 5 hours trying to sort out the imaging module, and failed miserably. I reckon this is something Rich is a lot better suited to, so I'll leave it to him. Nevertheless, I'll leave some notes on what I've done to the module and where I think it gets confusing:

• The imaging module basically consists of all the CCD data we use in the lens analysis. This is the Image from the CCD, as well as its Noise, BackgroundNoise, ExposureTime, PSF and potentially derived quantities like the SignalToNoise.

• Where things like an Image or Noise are alwyas 2D, the BackgroundNoise or ExposureTime could be a single float giving the value over the entire image, or a 2D image containing the noise / exposure time in each individual pixel.

• Previously, all unit conversions (e.g. x_pixel_to_x_arc_second) were module level functions, which require the pixel scale and image dimensions to be input into them. I have restructured this so that every class above (e.g. Image, Noise, etc.) inherit a DataGrid (e.g. units and conversions). This inheritance uses the pixel_scale and image dimensions for coordinate calculations internally.

• Previously, all classes inherited from np.ndarray via a new method. I am happy for you to go back to this, I changed it as I didnt know how to code with that style properly. From what I could tell, this made it difficult for the different classes to have different inputs and behave in different ways, which they will depending on the data at hand.

• Every Data class (e.g. the Image, Noise, ExposureTime) requires their own input pixel_scale and can be padded / trimmed independently. However, for a CCD, these different parts of the data all come from the same CCD, so they should all have the same 'shared' pixel scale and coordinates. There should be some way of grouping the different Data's into one class (e.g. CCDData) that does this.

• There are lots of different ways to set up each Data. For example, the Noise can be input, loaded from a fits file, estimated from an image, estimated from an Image + ExposureTime, estimated from an Image + Backgrond noise, etc. This may just require us to have lots of class methods for each circumstance.

• The specifics of handling each image in the point above is what makings the imaging module so hard to streamline. Ideally, one could just give the code a set of data and it runs, but each Data can be initialized in many different ways.

Be sure to read the post below after this one - most of what is in this post is useless

Once we have defined an image and masked it, we need to set up the sub-pixel grid on which deflection angles are computed and traced to the source-plane. I suggest we do this pre-analysis, once we have set up the mask. This grid will be defined by a list of lists, where each entry in the list contains tuples of the image-pixel sub-pixel coordinates, e.g. for a 2x2 sub-grid :

sub_grid[1][:] = [(1.575, 1.675), (1.675, 1.675), (1.675, 1.575), (1.575, 1.575))]

This means we do not need to recompute the image-pixel coordinates for the sub-grid for each iteration of the code. Map operations can then be used for tracing the image pixels to the source plane.

Source-plane sub-pixels

Thus far, the code has assumed that the source-plane sub-grid is a 1D list of size image_pixels*sub_pixels**2. Clearly, by adopting the list of list notation in the image-plane sub-grid, we can do so in the source-plane notation, and very much should. A few routines will require tweaking (e.g. the mapping of sub-pixels to source pixels and creation of the f matrix), but the effort required here is low. A benefit is that these routines and the code in general will require less inputs and provide a cleaner mapping of image and source coordinates.

Deflection angle interpolation

In the fortran code, I define 3 regions of sub-pixel deflection angle interpolations, depending on where in the mass profile we are (noting that close to the centre of the mass profile, the density is changing rapid and therefore needs to be computed more accurately):

1. The central regions (e.g. < 0.6") - here there is no interpolation, as the diverging mass profile would lead to inaccurate results.

2. The intermediate regions (e.g. 0.6" < 2.0" ) - here the the density is smoothing, so interpolation can be applied. To achieve this, I compute the deflection angle at the centre of each image pixel in this region, and interpolate with the sub-pixels of these imge pixels.

3. The outer regions (e.g. > 2.0") - Here the density is so flat, I compute the deflection angles on a sparse grid of image pixels (e.g. spaced 2 or 3 pixels apart) and interpolate the sub-pixels from this grid.

Adopting this scheme, we don't just store each sub-pixel image plane coordinate, but also the 4 image pixels we interpolate it with to get its value and the interpolation weights. This is quite memory intensive so I want to think of a better way to set this up. However, lets get AutoLens up and running before we develop this routine (lets just not bother with the interpolation for now).

Specific changes to existing routines

1. The routine sub_coordinates_to_source_pixels_via_sparse_pairs assumes the input list of sub_coordinates is a 1D list. We will change this to the structure above. It also assumes we have a list mapping sub pixels to sparse pixels (sub_pixel_to_sparse_pixel). We can change this to simply be the image_pixel_to_sparse_pixel, which is much easier to set up and less memory intensive (it wont be any slower).

2. The routine create_mapping_matrix. This uses two 1d lists sub_to_source and sub_to_imae to map sub-pixels to source pixels. sub_to_source is currently a 1D list, but will now become a list of lists. This will mean that we do not need to input the sub_grid_size into this routine, nor will we need to use sub_to_image.

The Sparse Grid and KMeans

The KMeans algorithm requires the input list of source-plane coordinates is 1D. Fortunately, this does not impact our sub-grid, as we define a sparse-grid of coordinates to compute this, which we will keep in 1D. This produces a 1D list of sparse_to_source_pixels, which is used to compute the sub-grid using the list structure above.

The Fortran code computed the sparse pixel deflection angles directly, via interpolation. We may be better served computing each sparse pixel source plane coordinate as the weighted sum of its four surrounding sub-pixels.

I setup a phase with two galaxies:

    phase1 = LensProfileGalaxy1Phase(lens_galaxies=[gp.GalaxyPrior(elliptical_sersic=lp.EllipticalSersicLP),
                                                                                        gp.GalaxyPrior(elliptical_sersic=lp.EllipticalSersicLP)])

I then over ride the first galaxy with constants in pass priors

    class LensProfileGalaxy1Phase(ph.LensProfilePhase):
        def pass_priors(self, previous_results):
            self.lens_galaxies[0] = previous_results[0].constant.lens_galaxies[0]

This means that the results from the analysis now only has one galaxy in lens_galaxies.variable, and two in lens_galaxies.constant. Furthermore, it means that what we call result.lens_galaxies[0].variable corresponds to result.lens_galaxies[1].constant.

Is it possible for the lens galaxy overirde to not completely remove the variable lens galaxy and keep the list structure in order? I can work around for now but its unintuitive.

The source PSF convolution routine currently goes to each pixel and computes the pixels within the mask that it needs to be blurred with during PSF convolution (it does this in 2D, but we agreed 1D would be more efficient). It does this by overlaying the kernel over the pixel and book-keeping all pixels within the mask and kernel.

However, a lot of the image-pixels in this mask will have essentially no source-light in them (the mask is big because of the lens galaxy's size). For these pixels, there is no need to use a big PSF kernel, as they really don't contain much information on the source. So, lets implement a routine which reduces the kernel size for image-pixels we know contain no source light. This should be straight-forward, as it should just come down to changing the kernel size input to each 'frames' you use to setup the convolution pre-analysis.

Determining image pixels without source light is simple - we just look at the model of a previous run of the code (we do this anyway for the adaptive source reconstruction features).

We need classes for storing the Lens and Source galaxy(s), for two reasons:

• The lens and source galaxies each have a set of observation properties that do not change depending on the mass / light models (e.g. their redshifts, cosmological distances, velocity dispersions etc.). So, we want a class (Something like LensGalaxy) which the mass and light profiles inherit from, so we can change their models without their overall properties. Will this require the LensGalaxy class to have constructors for the mass model?

• For simplicity, we may want to evaluate the light profile or deflection angles of all the lens galaxyies (imagine a weird system where the 'lens' is two merging galaxies). Thus, the LensGalaxy class should then have some functions 'compute all deflection angles' and 'compute all light profiles'.

Integration using functions such as quad raises the warning that the function argument should be a low level callable. e.g. in line 157 of mass_profile.py:

potential = quad(self.potential_func, a=0.0, b=1.0, args=(coordinates,))[0]

self.potential function should be a LowLevelCallable.

The docs describe how a compiled c function can be used in this respect.

However, SciPy supports using python functions written in a very particular as LowLevelCallables.

Ideally all deflection functions should be modified or duplicated to be used as LowLevelCallables.

This would make calling these functions confusing. It may make sense to define any function to be used in integration as an inner function of the function that performs the integration.

In the development branch, you should now be able to find a script 'grid_demo.py' in the scripts folder which explains how the profiles, galaxy, grids and ray tracing modules work. If things are unclear, just ask me in this GitHub issue.

The end of the script also has a discussion of how I would implement the potential grid, again just ask me for clarity (And I might be misunderstanding something important).

The code as I've shown it needs a bit of refactoring, some of the classes get a bit clunky with their long attribute names (the word grids comes up a lot). And try not to think too much about optimization at the moment, we'll worry about that once we can profile the code.

If we pair a two floats:

    class MMPhase(ph.LensPlanePhase):

        def pass_priors(self, previous_results):

            self.lens_galaxies[0].sersic.intensity = self.lens_galaxies[0].sersic.axis_ratio

The model.info file does not reflect this pairing:

EllipticalSersic

lens_galaxies_centre_0 GaussianPrior, mean = 0.0, sigma = 0.5
lens_galaxies_centre_1 GaussianPrior, mean = 0.0, sigma = 0.5
lens_galaxies_intensity UniformPrior, lower_limit = 0.5, upper_limit = 1.0
lens_galaxies_axis_ratio UniformPrior, lower_limit = 0.0, upper_limit = 180.0
lens_galaxies_phi UniformPrior, lower_limit = 0.5, upper_limit = 1.0
lens_galaxies_effective_radius UniformPrior, lower_limit = 0.0, upper_limit = 2.0
lens_galaxies_sersic_index UniformPrior, lower_limit = 0.8, upper_limit = 8.0

Checkout 'pair_floats.py' integration test.

If we pair a tuple parameter:

    class MMPhase(ph.LensPlanePhase):

        def pass_priors(self, previous_results):

            self.lens_galaxies[0].sersic.centre_0 = self.lens_galaxies[0].sersic.axis_ratio

model.info does not reflect this:

EllipticalSersic

> lens_galaxies_centre_1                  UniformPrior, lower_limit = 0.5, upper_limit = 1.0
> lens_galaxies_axis_ratio                GaussianPrior, mean = 0.0, sigma = 0.5
> lens_galaxies_centre_0                  UniformPrior, lower_limit = 0.5, upper_limit = 1.0
> lens_galaxies_phi                       UniformPrior, lower_limit = 0.0, upper_limit = 180.0
> lens_galaxies_intensity                 UniformPrior, lower_limit = 0.0, upper_limit = 5.0
> lens_galaxies_effective_radius          UniformPrior, lower_limit = 0.0, upper_limit = 2.0
> lens_galaxies_sersic_index              UniformPrior, lower_limit = 0.8, upper_limit = 8.0

And neither does mn.results:

Most likely model, Likelihood = -774.3340543069343

lens_galaxies_centre_1                            0.025323588465887274
lens_galaxies_centre_0                            0.7339913315711855
lens_galaxies_phi                                 156.04476328513027
lens_galaxies_intensity                           2.133283347234559
lens_galaxies_effective_radius                    1.9811613864886086
lens_galaxies_sersic_index                        0.8314012538231821

Most probable model (3 sigma limits)

lens_galaxies_centre_1                            0.025841115091876164 (-0.03571302847822744, 0.0934866142252781)
lens_galaxies_centre_0                            0.7243205070654029 (0.6521633755644127, 0.7553890092227087)
lens_galaxies_phi                                 156.7717683276137 (146.38954660304958, 172.1219018122385)
lens_galaxies_intensity                           2.11720631837675 (2.005912715932006, 2.2152704517064503)
lens_galaxies_effective_radius                    1.9794334007339063 (1.9398774306863937, 2.011765761578133)
lens_galaxies_sersic_index                        0.8373621306262105 (0.7776963481590332, 0.9287504513654512)

Most probable model (1 sigma limits)

lens_galaxies_centre_1                            0.025841115091876164 (0.0054955320593991774, 0.045710262533824375)
lens_galaxies_centre_0                            0.7243205070654029 (0.7184671038499029, 0.7397098989676293)
lens_galaxies_phi                                 156.7717683276137 (152.71842853013382, 159.51021196442056)
lens_galaxies_intensity                           2.11720631837675 (2.0859044139025413, 2.1516978255446717)
lens_galaxies_effective_radius                    1.9794334007339063 (1.9686753917960578, 1.9913354624475075)
lens_galaxies_sersic_index                        0.8373621306262105 (0.8142795187064701, 0.8559290797237339)

Constants

Checkout 'pair_tuples_x1.py' integration test.

However, if we pair both parameters in the tuple, the issue in mn.results goes away

There are also issues when both tuples are paired, checkout pair_tuples_x2.py'

dimensionless_mass_within_ellipse is not implemented for:

EllipticalGeneralizedNFW
SphericalGeneralizedNFW
EllipticalNFW
SphericalNFW

In computer imagery (0, 0) generally refers to the top left point in an image.

In the image module we are using the convention that the central image pixel has arcsecond coordinates of (0, 0) with the (0, 0) pixel coordinates being negative x and y. This implicitly assumes that images defined in arcseconds use the same convention as computer imagery. It is worth being aware of this assumption as it may cause results to be rotated from what we would expect them to be.

This follows on my post in #19, where I suggest we make RayTracingPlane are core model class on which our priors module operators.

In the priors module, I added two more routines, 'reconstruction_most_likely' and 'reconstructon_most_probable'. These essentially load the MultiNest output and setup a model using it. Currently, they assume the ClassMap has already been set up, so it uses an existing instance of class map to setup the model.

We actually want this to be a layer 'above' ClassMap e.g., given a directory of MultiNest results, the code automatically reads the MultiNest output files and knows how to set up the RayTracingPlane class (this will require us to make a few new output files containing the model information). It'd do this using the ClassMap module, thus setting up the model as one normally would without MultINest, This is neat, as one can go to a set of results and instantly set up the model as the results of an analysis.

So, I guess I'm proposing a 'MultiNest' layer which handles the interface between analysis results and model setup. This would allow us to easily integrate other analysis tools instead of MultiNest.

I have already started doing this when I made a module called MultiNestTools. I think actually this shouldn't be 'Tools' but a layer of code, with a direct interface with ClassMap as discussed above. I'll begin writing code that outputs the details of a model to the MultiNest dir, which the MultiNest layer can read to set up the RayTracingPlane.

An analysis is performed within a mask, and therefore produces model images only within this mask. The model image (and its galaxy images) are used to produce the hyper images. If the mask changes between phases, the model images must be mapped to 2D, and then back to 1D using the new mask. This is currently what phase does.

However, this still isn't enough, as if the mask grows between phases (which is probably the expected behaviour), there will be zeros in the regions where the analysis was masked previously. This will mess up the use of hyper images.

For profiles, we can simply compute the profile image on an unmasked 2D image, and use that to setup the 1D hyper galaxy images.

For pixelizations however, using an unmasked 2D image will mess up the reconstructed image and this approach will not be viable. Thus, we can employ the following two rules:

• If no pixelization is used, we are free to change mask and have hyper galaxies be set up correct.
• If a pixelization is used, the mask must remain fixed thereafter.

I have reverted the code to use only 1D images for now and ignore 2D mapping. I will build pipelines for now that don't change masks after a hyper phase. Hopefully I can get round to this soon, but it doesnt feel like a big enough issue atm.

Something is going wrong with getdit, whereby when it reads the multinest.paramnames file of the 2_x2_galaxies howtolens pipeline incorrectly.

I believe this is tied to whether getdist decides a parameter 'counts'. If it feels the samples of a parameter are too narrow or have too small a range it chooses to ignore it 0_0. This is a nightmare, and would mean we need to stop using getdist for parameter linking.

Ive put a raise error in non_linear.MultiNest.pdf, to check whether getdist has the same number of parameters as MultiNest, so that when this happens we know...

Lets say we are modeling a lens galaxy's light profile using a double Sersic profile:

sersic1 = EllipticalSersic(...)
sersic2 = EllipticalSersic(...)
light_profile = CombinedLightProfile(sersic1, sersic2)

If the two light profiles are completely independent, this means each is controlled by 7 nonlinear parameters, giving us a total of 14 non-linear parameters.

In many modeling situations, it is unclear if the two Sersic's should share the same geometry. For example, we might want the second Sersic to share the same centre as the first sersic. This must be reflected in the priors we give sersic2, as it no longer draws its values of the centre from a prior but simply acquires the values from sersic1. The same could happen for the rotation angle and ellipticity. Thus, we need the code to be able to know if the centers of two light profiles are shared, that there are now just 12 non-linear parameters.

This could feasibly happen when there are 3 Light profiles (e.g. a Sersic + Sersic + Sersic) and we need flexibility in choosing if sersic3 and sersic2 are aligned with sersic1, or not.

This also occurs for mass profiles. For example, we may have a dark matter halo given by an EllipticalNFWProfile, and may want to align it with the centre or rotation angle of our CombinedLightProfile (somehow picking which component we align it with). Thus, for any mass profile, we need the option too tie its geometry to given component of a CombinedLightProfile.

Mass profiles also have some level of shared parameters. For example, a SersicMassProfile has a mass to light ratio. If a CombinedMassProfile is formed of two SersicMassProfiles, the question arises whether they share the same mass_to_light ratio, or use indepedent values.

The model.info file displays the priors of an analysis, e.g.:

EllipticalSersic

lens_galaxies_centre_0 GaussianPrior, mean = 0.1936895385006247, sigma = 0.02
lens_galaxies_centre_1 GaussianPrior, mean = 0.19347809618045228, sigma = 0.02
lens_galaxies_axis_ratio GaussianPrior, mean = 0.5953876055651927, sigma = 0.1
lens_galaxies_phi GaussianPrior, mean = 48.674970226754766, sigma = 10.0
lens_galaxies_intensity GaussianPrior, mean = 4.944301653522131, sigma = 0.2249443899895267
lens_galaxies_effective_radius GaussianPrior, mean = 0.7043911318725903, sigma = 1.0
lens_galaxies_sersic_index GaussianPrior, mean = 1.2070190682963673, sigma = 1.2

However, if there are constant parameters, it does not list them. It'd be nice if they were listed as well, e.g.:

Variable:

EllipticalSersic

lens_galaxies_centre_0 GaussianPrior, mean = 0.1936895385006247, sigma = 0.02
lens_galaxies_centre_1 GaussianPrior, mean = 0.19347809618045228, sigma = 0.02
lens_galaxies_axis_ratio GaussianPrior, mean = 0.5953876055651927, sigma = 0.1
lens_galaxies_phi GaussianPrior, mean = 48.674970226754766, sigma = 10.0
lens_galaxies_intensity GaussianPrior, mean = 4.944301653522131, sigma = 0.2249443899895267
lens_galaxies_effective_radius GaussianPrior, mean = 0.7043911318725903, sigma = 1.0
lens_galaxies_sersic_index GaussianPrior, mean = 1.2070190682963673, sigma = 1.2

Constant:

SphericalIsothermal

lens_galaxies_centre_0 0.0
lens_galaxies_centre_1 1.0
lens_galaxies_einstein_radius 1.0

The reconstruction_for_vector in the priors module works perfectly, assigning physical values from a unit hypercube as expected to setup a reconstruction class. This includes for cases where the priors are updated after being loaded from the config file or set equal to one another - I've added tests showing this.

However, the value_vector_for_hypercube_vector still swaps the physical values that are output, if the priors are either updated manually after the initial set up or set equal to one another. This is still problematic, as the unit hypercube and physical hypercube must maintain the same order. I have added a unit test showing the problem, which I've pasted below:

        def test__order_maintained_with_prior_change(self):

            collection = prior.ClassMap(
                prior.Config(config_folder_path=data_path+"config_test"),
                profile_1=geometry_profiles.EllipticalProfile, profile_2=geometry_profiles.Profile,
                profile_3=geometry_profiles.EllipticalProfile)

            collection.profile_1.axis_ratio = prior.UniformPrior(100, 200)

            assert collection.value_vector_for_hypercube_vector([0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]) == \
                   [150.0, 0.8, 0.5, 0.25, 0.5, 0.25, 0.75, 0.8, 0.5, 0.25]

The error is that the 150.0 entry does not appear first in the output list, but last.

I am confident I have the tools I need to fix this myself - infact I have already hacked together a solution (in the CTI code) which addresses cases like the one above - by moving the physical values to their correct entries in the list. I haven't yet solved this for settings priors equal to one another (e.g. collection.profile_1.axis_ratio = collection.profile_2.axis_ratio )

All Spherical profiles inherit from their elliptical counterpart. This means that they don't use the SphericalProfile transform grid function, which is faster than the EllipitcalProfile equivalent.

We should change our inheritance structure to address this and add unit tests testing the inheritance