Evaluation on the test set

This document evaluates one trained machine learning model, see details below. We use both conventional metrics and comparison with a resimulated image.

We use three datafiles, which must be downloaded if this document is to be executed.

testing.csv is ground truth for model training
test.csv is the predicted amplitudes from machine learning.
dataset.csv is the original lens parameters used to generate the training, testing, and validation data, i.e. it has more rows than the other two sets.

The neural network used for machine learning is almost arbitrarily chosen and has not been tuned. The architecture is one of the best performing in Nicolò’s experiments on other datasets from CosmoSim, but the hyperparameters are arbitrarily chosen. The training set used is 16000 images.

The specification of the distribution is discussed in Sample Datasets and can be downloaded (dataset.toml).

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( "Neural network used:", cfg["settings"]["model"] )
print( json.dumps( cfg["hyperparameters"], indent=4 ) )

CosmoSim version is 3.2.2b2
Data from directory experiment001
Neural 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-012001.png    2.349062e-11
image-012002.png    7.091439e-11
image-012003.png    2.546456e-10
image-012004.png    1.730464e-10
image-012005.png    5.698352e-12
                        ...     
image-015996.png    1.359637e-11
image-015997.png    4.388787e-11
image-015998.png    5.255600e-11
image-015999.png    5.295995e-11
image-016000.png    1.193498e-11
Length: 4000, dtype: float64

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

sse.nlargest(3)

filename
image-013232.png    4.497526e-10
image-012946.png    4.085509e-10
image-013180.png    3.678154e-10
dtype: float64

sse.nsmallest(3)

filename
image-014934.png    8.154099e-14
image-012366.png    2.576096e-13
image-012399.png    3.433474e-13
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.nlargest(3).index)
worst = list(sse.nsmallest(3).index)
print( "Best:", best )
print( "Worst:", worst )

Best: ['image-013232.png', 'image-012946.png', 'image-013180.png']
Worst: ['image-014934.png', 'image-012366.png', 'image-012399.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
[getSource] mode=SourceSpec.Sphere, ltprf=LightProfileSpec.Sersic
[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
[getSource] mode=SourceSpec.Sphere, ltprf=LightProfileSpec.Sersic
[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, fn, "Ground Truth", axiscross=True )

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

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

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

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

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

[getSource] src=SersicSphere, ltprf0=None, verbose=1
[getSource] mode=SourceSpec.Sphere, ltprf=LightProfileSpec.Sersic
[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
[getSource] mode=SourceSpec.Sphere, ltprf=LightProfileSpec.Sersic
[SphericalSource] constructor done

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

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

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

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

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

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"]["model"] = "Raytrace"
cfg["simulator"]["centred"] = False
cfg["lens"] = { "mode" : "SIE" }
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
[getSource] mode=SourceSpec.Sphere, ltprf=LightProfileSpec.Sersic
[SphericalSource] constructor done

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

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

The match is not perfect, but we do see that the primary image is correctly placed. The shape is not accurate on the edges. It is significant that the shadow in the difference image around the primary image is either black or white, never both. This means that the image is smaller or larger, and not rotated or awfully misshaped in general. Clearly, this is a limitation of the roulette formalism, since there is no difference between reconstructed and ground truth 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
[getSource] mode=SourceSpec.Sphere, ltprf=LightProfileSpec.Sersic
[SphericalSource] constructor done

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

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

Interestingly, the best images do not perform any better than the worst in terms of visual comparison between roulettes and raytrace.

Conclusion¶

We see that this machine learning model make accurate prediction as far as optical perception goes, but there are limitations to the roulette representations.

This means that there is nothing to gain from further research on machine learning models at this stage.