CosmoSim Demo III Roulette Resimulation

This Demo assumes that you are familiar with the principles from CosmoSim Demo I. We start with the same imports

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import tomllib as tl

from CosmoSim.datagen import SimImage
import CosmoSim.Image as csimg
from CosmoSim import Parameters

import CosmoSim as cs
print( "CosmoSim version", cs.__version__ )

CosmoSim version 3.3.0

Additionally, we require the Resim class for resimulation.

from CosmoSim.roulettegen import Resim

The Resim class shares the GenericSim superclass with SimImage, which we used in CosmoSim Demo I but simulates from precomputed roulette parameters and not from a regular lens model.

Basic Configuration¶

We start with parameters similar to the ones from Demo01.md, to make a baseline raytrace simulation using a SIE lens. We have just rounded some of the parameters to get neater numbers. We will use the simulation from SimImage to get roulette amplitudes for resimulation.

with open( "Demo03.toml", 'rb' ) as f:
            toml = tl.load(f)

param = Parameters( toml )
imsim = SimImage( param, verbose=0 )
im = imsim.getImage()
csimg.imshow( im, "Baseline simulation" )

Tip

For debugging, it may be useful to change the verbosity level on the fly. On this point, the C++ code functions differently from the python code. Where the python code often has per-class and per-function verbosity settings, the C++ code uses a single global variable which can be set thus:

from CosmoSim.CosmoSimPy import setDebug
setDebug(3)

This is set every time the Resim or SimImage classes are instantiated, which works well enough if only one instance is made, or all use the same verbosity level. For more complex situations, the above setDebug may be useful. We used it during debugging, but it has now been from the final demo.

Roulette parameters¶

We can retrieve the roulette amplitudes with the getData() method. The resulting data is a pandas Series.

row = imsim.getData()
display( row )

filename              test.png
source         SersicEllipsoid
x                           12
y                            6
sigma                       15
sigma2                      35
theta                       45
luminosity                  70
n_sersic                     4
lensX                        0
lensY                        0
centreX                      0
centreY                      0
reletaX                   12.0
reletaY                    6.0
offsetX               0.007619
offsetY              49.989834
xiX                  14.328906
xiY                  58.311121
alpha[0][1]          -2.321287
alpha[1][0]          -0.293492
alpha[1][2]          -0.260066
alpha[2][1]          -0.001387
alpha[2][3]           0.004043
alpha[3][0]          -0.000131
alpha[3][2]           0.000201
alpha[3][4]            0.00022
alpha[4][1]          -0.000001
alpha[4][3]           0.000003
alpha[4][5]          -0.000012
alpha[5][0]               -0.0
alpha[5][2]           0.000001
alpha[5][4]               -0.0
alpha[5][6]          -0.000001
beta[0][1]          -52.312994
beta[1][0]                 0.0
beta[1][2]            0.136027
beta[2][1]            0.007891
beta[2][3]             0.00621
beta[3][0]                 0.0
beta[3][2]            0.000057
beta[3][4]           -0.000191
beta[4][1]            0.000011
beta[4][3]           -0.000008
beta[4][5]            -0.00001
beta[5][0]                 0.0
beta[5][2]                 0.0
beta[5][4]                -0.0
beta[5][6]            0.000001
Name: test.png, dtype: object

Here we see that we have roulette amplitudes up to order 5, which is the maximum implemented for analytical SIE.

rsim = Resim( row, verbose=0 )
resimImage = rsim.getImage()
csimg.imshow( resimImage, "Resimulation" )

[getSource] src=SersicEllipsoid, ltprf0=None, verbose=1

For reference, we can do the original roulette simulation, using a SimImage object.

param["simulator"]["model"] = "Roulette"
roulette = SimImage( param, verbose=0 )
rouletteImage = roulette.getImage()
csimg.imageCompare( rouletteImage, resimImage, "Roulette simulation", "Resimulation" )

The two roulette simulations look similar but displaced compared to each other.

The Source¶

We can also have a look at the original source image.

actual = rsim.getActualImage()
csimg.imshow( actual, "Actual source image" )

Comparison¶

To compare the images, we can plot side by side We also add reference lines to each image, except the difference image. If the code is confusing, it is because we

imdiff = csimg.imageDiff(rouletteImage,resimImage)
rdiff = csimg.imageDiff(rouletteImage,im)

csimg.drawAxes(actual)
csimg.drawAxes(im)
csimg.drawAxes(rouletteImage)
csimg.drawAxes(resimImage)

fig = plt.figure(figsize=(15,10))
fig.tight_layout(pad=0.0)
plt.subplots_adjust(hspace=0.1, wspace=0.1) 

fig.add_subplot(2, 3, 1)
csimg.imshow( actual, "Actual source image" )
fig.add_subplot(2, 3, 2)
csimg.imshow( im, "Original raytrace simulation" )
fig.add_subplot(2, 3, 3)
csimg.imshow( rdiff, "Roulette/raytrace difference" )
fig.add_subplot(2, 3, 5)
csimg.imshow( imdiff, "Difference image" )
fig.add_subplot(2, 3, 4)
csimg.imshow( rouletteImage, "Roulette simulation" )
fig.add_subplot(2, 3, 6)
csimg.imshow( resimImage, "Resimulation" )

Centring the image¶

In the original simulation, centring is straight forward, by setting simulator.centred to True.

param["simulator"]["centred"] = True
param["simulator"]["model"] = "Raytrace"
imsim = SimImage( param, verbose=0 )
imC = imsim.getImage(verbose=1)

fig = plt.figure(figsize=(10,5))
fig.tight_layout(pad=0.0)
plt.subplots_adjust(hspace=0.1, wspace=0.1) 

fig.add_subplot(1, 2, 1)
csimg.imshow( imC, "Original simulation - centred" )
fig.add_subplot(1, 2, 2)
csimg.imshow( csimg.imageDiff(imC,im), "Difference" )

[getImage] centring from parameters: True

We have plotted the difference image to see that centring does make a difference. We can also note the difference it makes to the data for resimulation.

rowC = imsim.getData()
display( pd.DataFrame( [row, rowC] ).transpose() )

Loading...

The roulette amplitudes do not change. The are calculated in the same point. However the relative source position reletaX/reletaY does. It is the sum of the original position and the centrepoint centreX/centreY.

We can do the same with roulette; centred is still True.

param["simulator"]["model"] = "Roulette"
roulette = SimImage( param, verbose=0 )
rouletteImageC = roulette.getImage(verbose=0)

fig = plt.figure(figsize=(10,5))
fig.tight_layout(pad=0.0)
plt.subplots_adjust(hspace=0.1, wspace=0.1) 

fig.add_subplot(1, 2, 1)
csimg.imshow( rouletteImageC, "Roulette simulation - centred" )
fig.add_subplot(1, 2, 2)
csimg.imshow( csimg.imageDiff(rouletteImageC,rouletteImage), "Difference" )

Finally we resimulate from the data from the centred simulatior. We compare the result with both the original resimulation and with the centred raytrace.

rsim = Resim( rowC, verbose=0 )
resimImageC = rsim.getImage()

fig = plt.figure(figsize=(15,5))
fig.tight_layout(pad=0.0)
plt.subplots_adjust(hspace=0.1, wspace=0.1) 

fig.add_subplot(1, 3, 1)
csimg.imshow( resimImageC, "Resimulation from centred image" )
fig.add_subplot(1, 3, 2)
csimg.imshow( csimg.imageDiff(resimImageC,resimImage), "Difference with previous resimulation" )
fig.add_subplot(1, 3, 3)
csimg.imshow( csimg.imageDiff(resimImageC,imC), "Difference with centred raytrace" )

[getSource] src=SersicEllipsoid, ltprf0=None, verbose=1

We note that the image has shifted compared to the original resimulation.
The right hand image shows that the co-ordinate system is shifted to match the centred raytrace image.

Closure¶

Resmulation is going to be critical to evaluate machine learning models in practical scenarioes. The resimulated image can be compared to the observation without having a ground truth for the roulette amplitudes.

There is some work in progress here.
I hope we have managed to demonstrate some of the caveats, particularly with respect to centring. Documentation of Implementation and the definitions of the different co-ordinate systems is ongoing work.