Machine Learning on Cluster Lenses

This document evaluates one trained machine learning model, see details below. We use both conventional metrics and comparison with a resimulated image. There is one methodological flaw in the model training. The early stopping criterion is based on testing on the same test set as we use here.

Configuration¶

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import IPython.display as idp
import toml
import json
import CosmoSim.Image as csimg
import CosmoSim as cs
import CosmoSim.datagen as dg
import CosmoSim.roulettegen as rg
imshow = csimg.imshow

Let us check the parameters underlying this test. This gives the neural network and hyperparameters used.

cfg = toml.load( "ml.toml" )
print( "CosmoSim version is", cs.__version__ )
print( "Data from directory", cfg["output"]["directory"] )
print( "Nerual network used:", cfg["settings"]["model"] )
print( json.dumps( cfg["hyperparameters"], indent=4 ) )

CosmoSim version is 3.2.4b1
Data from directory experiment002
Nerual network used: regnet_y_32gf
{
    "learning-rate": 0.001,
    "batch-size": 32,
    "weight-decay": 0.0001
}

Let’s load the test set, both the ground truth (gt) and the predictions (df).

gt = pd.read_csv( "../testing.csv", index_col="filename" )
df = pd.read_csv( "test.csv", index_col="filename" )
display( gt.head() )
display( df.head() )

As we can see, we have columns for all the parameters in question, and rows for the individual objects. There is no superfluous data in the data frames.

We can also check the difference.

rawerrors = gt - df
display( rawerrors.head() )

That was straight forward. Now we come to the tricky part - making sense of the numbers. There are two questions that we want to consider.

How well does the model perform overall?
How do particular images behave, and what are the best and the worst samples?

Model evaluation¶

It is interesting to see how the model behaves on each column. To see this, we can aggregate the rows, computing the mean of absolute errors and mean of squared errors.

mse = (rawerrors**2).mean()
mae = rawerrors.abs().mean()
mse.name = "MSE"
mae.name = "MAE"
colerrors = pd.concat( [ mse, mae ], axis=1 )
display( colerrors )

In practice, relative errors may be more interesting. We can calculate that too.

colerrors["mean"] = gt.abs().mean()
colerrors["Relative MAE"] = colerrors["MAE"] / colerrors["mean"]
display( colerrors )

This looks good with mean relative errors less than 10^-7, and no column is particularly bad or good.

Object evaluation¶

To investigate individual objects, we start by calculating sum of squared errors for each one.

sse = (rawerrors**2).sum(axis=1)
display(sse)

filename
image-010001.png    5.349902e-11
image-010002.png    1.096119e-10
image-010003.png    1.376738e-10
image-010004.png    7.306417e-11
image-010005.png    1.553233e-10
                        ...     
image-013996.png    1.033530e-10
image-013997.png    1.716025e-10
image-013998.png    3.867597e-11
image-013999.png    3.155409e-10
image-014000.png    1.530504e-10
Length: 4000, dtype: float64

We can look at the three best and three worst images, thus,

sse.nlargest(3)

filename
image-011941.png    9.255588e-10
image-011530.png    8.878920e-10
image-010865.png    8.739316e-10
dtype: float64

sse.nsmallest(3)

filename
image-012974.png    1.132337e-12
image-011907.png    2.474771e-12
image-010461.png    3.606011e-12
dtype: float64

We note that the errors are small, but there is also a huge span between the best and the worst. For the purpose of this test, we do not assume that we have access to the original images, but we can resimulate them from the roulette amplitudes. First we record the filenames.

best = list(sse.nsmallest(3).index)
worst = list(sse.nlargest(3).index)
print( "Best:", best )
print( "Worst:", worst )

Best: ['image-012974.png', 'image-011907.png', 'image-010461.png']
Worst: ['image-011941.png', 'image-011530.png', 'image-010865.png']

Simulation from Roulette Amplitudes¶

Before we can simulate, we need to set up some basic parameters.

cfg = { "simulator" : { "imagesize" : 512, "cropsize" : 256, "xireference" : True }
      , "source" : { "mode" : "SersicSphere" } }
param = cs.Parameters( cfg )

Lets first consider the worst images and the ground truth. The data is extracted as

gtw = gt.loc[ worst ]
display( gtw )

To simulate, we instantiate the Resim class which is the simulator for roulette resimulation. For now, we test in just one image.

imsim = rg.Resim(gtw.iloc[0],param=param,verbose=0)
im = imsim.getImage()
imshow( im )

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

We can make the same simulation from the reconstructed amplitudes and compare.

dfw = df.loc[ worst ]
imsim2 = rg.Resim(dfw.iloc[0],param=param,verbose=0)
im2 = imsim2.getImage()
csimg.imageCompare( im, im2, "Ground Truth", "Reconstructed" )

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

This looks like perfect match.

We can do the same for all the top and bottom three.

for fn in worst:
    dfsim = rg.Resim(df.loc[fn],param=param,verbose=0)
    dfim = dfsim.getImage()
    gtsim = rg.Resim(gt.loc[fn],param=param,verbose=0)
    gtim = gtsim.getImage()
    csimg.imageCompare( dfim, gtim, "Reconstructed", "Ground Truth" )

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

No visible discrepancy. We can continue with the best images, obviously expecting perfect match again.

for fn in best:
    dfsim = rg.Resim(df.loc[fn],param=param,verbose=0)
    dfim = dfsim.getImage()
    gtsim = rg.Resim(gt.loc[fn],param=param,verbose=0)
    gtim = gtsim.getImage()
    csimg.imageCompare( dfim, gtim, "Reconstructed", "Ground Truth" )

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

Numerical errors¶

gtw = gt.loc[ worst ]
dfw = df.loc[ worst ]
display( gtw - dfw )

display( (gtw - dfw)/gtw )

There are some numbers that stick out as particularly large, but that’s still on the order of 10^-5 at worst.

gtb = gt.loc[ best ]
dfb = df.loc[ best ]
display( gtb - dfb )

display( (gtb - dfb)/gtb )

Simulations from Lens Parameters¶

We can also load the original lens parameters, from which the ground truth was computed.

orig = pd.read_csv( "../dataset.csv", index_col="filename" )
display( orig.head() )

This dataset is different from the others, giving physical parameters of the lens instead of roulette amplitudes in a point of observation. Thus, we need a different simulator. On a positive note, we can use raytrace simulation, which is accurate.

The SimImage simulator is parameterised in a slightly different way. We need to add the row data from the dataset to the Parameters object instead of passing it as a separate argument. To avoid interference, we make a copy of params.

cfg["simulator"]["config"] = "raysie"
cfg["simulator"]["centred"] = True
p2 = cs.Parameters( cfg )
for fn in worst:
    dfsim = rg.Resim(df.loc[fn],param=param,verbose=0)
    dfim = dfsim.getImage()
    p2.setRow( orig.loc[fn] )
    gtsim = dg.SimImage(param=p2,verbose=0)
    gtim = gtsim.getImage()
    csimg.imageCompare( dfim, gtim, fn, "Original raytrace simulation", axiscross=True )

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

This is not good, and with the small visible images, it is difficult to tell what the problem is.

Finally, we can also plot the best objects with raytrace simulation.

for fn in best:
    dfsim = rg.Resim(df.loc[fn],param=param,verbose=0)
    dfim = dfsim.getImage()
    p2.setRow( orig.loc[fn] )
    gtsim = dg.SimImage(param=p2,verbose=0)
    gtim = gtsim.getImage()
    csimg.imageCompare( dfim, gtim, fn, "Original raytrace simulation", axiscross=True )

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[SphericalSource] constructor done

Conclusion¶

We see that this machine learning model make accurate prediction as far as optical perception goes. Other datasets may prove harder, but this dataset gives no room for further tuning.