SIS Amplitudes

This test follows up on problems arising from Custom Amplitudes Files and tests calculation of roulette amplitudes from python.

Preparation¶

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
import toml
from copy import deepcopy

import CosmoSim as cs
from CosmoSim import Parameters
from CosmoSim.datagen import SimImage
from CosmoSim.roulettegen import Resim
import CosmoSim.Image as csimg

print( "CosmoSim version:", cs.__version__ )

CosmoSim version: 3.2.2b2

We use the same configuration as we have used before.

cfg = { 'simulator' : {
             "model" : "Raytrace",
             "centred" : False,
             "xireference" : True,
             "nterms" : 5, 
             "cropsize" : 256 }
      , 'lens': { 
            'mode' : "SIS",
            "amplitudefile" : "sis50.txt",
            'einsteinradius': 46 }
      , 'source': {
            'mode': 'SersicSphere',
            'sigma': 20,
            'luminosity' : 20,
            'position': 'cartesian'}
      , 'position': {'x': 10, 'y': 5 }
      }
param = Parameters(cfg)

We set up two simulators. We show only the raytrace simulation for now. This will be used as a reference to assess the fidelity of other simulations. Even though the simulator is set up for raytrace, it will also compute the roulette amplitudes for us. We centre the image, as this is the standard mode of operation in practice, and thus best tested.

sim = SimImage( param, verbose=0 )
ray = sim.getImage()
csimg.imshow( ray, 'Raytrace')

Roulette Amplitudes¶

As we have shown in previous demos, we can retrieve all the amplitudes as a pandas Series object.

df01 = sim.getData()

This calculation is made in double precision in the C++ library. The python code also offers a calculation in arbidtrary precision. The default is 64 digits at present. This computation is slow, so do not worry if you do not see the results immediately.

df02 = sim.getData(precision=64,verbose=0)

[['0', '1', '-g*x/sqrt(x**2 + y**2)', '-g*y/sqrt(x**2 + y**2)'], ['1', '0', '-0.5*g*(1.0*x**2 + 1.0*y**2)/(x**2 + y**2)**(3/2)', '0']]

We can compare the two, but we need a little trick to disregard non-numeric entries.

df = pd.DataFrame( [ df01, df02, df02.drop(["source","filename"])-df01 ],
                   index=[ "C++", "Python", "Difference" ] ).transpose()
display(df)

We can quickly check the worst discrepancy:

dfnumeric = df.drop( [ "source", "filename" ] )
display(dfnumeric.head())

diff = dfnumeric["Difference"].abs()
print( "Absolute error:", diff.max() )

Absolute error: 41.14365078599613041392879550465548273210737781685207332658450932

py = dfnumeric["Python"].abs()
py = py.mask( py == 0 )
df["Relative difference"] = diff/py
display( df )

We ignore the last amplitudes file that we have used, because it depends on the redundant orientation parameter, and the library does not currently handle this.

Resimulation¶

param1 = deepcopy( param )
resim1 = Resim( df01, param=param1, verbose=0 )
im01 = resim1.getImage()
csimg.drawAxes( im01 )
csimg.drawAxes( im01 )
csimg.imageCompare( im01, ray, "Resimulation", "Raytrace" )

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[getSource] mode=SourceSpec.Sphere, ltprf=LightProfileSpec.Sersic
[SphericalSource] constructor done

param2 = deepcopy( param )
resim2 = Resim( df02, param=param1, verbose=0 )
im02 = resim2.getImage()
csimg.imageCompare( im02, ray, "Resimulation", "Raytrace" )

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[getSource] mode=SourceSpec.Sphere, ltprf=LightProfileSpec.Sersic
[SphericalSource] constructor done

Further inspection¶

xi = sim.sim.getNu()
print(np.array(xi))

[51.14365079 25.57182539]

Preparation¶

Roulette Amplitudes¶

Resimulation¶

Further inspection¶

Conclusion¶