Interpolation of RF predictions with cosZD, for homogeneous performance #1320

moralejo · 2024-11-27T09:55:58Z

This should solve the issue of RF "performance jumps" at high zenith angles, when the pointing goes through the middle points between training nodes. See #1317

NOTE: this is a different implementation of the interpolation approach proposed by @gabemery, see branch https://github.com/cta-observatory/cta-lstchain/tree/dl2_RF_interpolate)

If the interpolation option is activated we call the RF predictors twice, for each of the two closest MC training nodes, then interpolate (or extrapolate) the values to the actual telescope pointing for each event.

Currently the training sample pointings are obtained from the path to the training sample(provided via config file). A better solution would be to add to the .sav files an array with the zenith and azimuth values of the MC training nodes.

to avoid jumps when using RFs trained on a discrete set of pointings

codecov · 2024-11-27T10:10:21Z

Codecov Report

Attention: Patch coverage is 43.85965% with 64 lines in your changes missing coverage. Please review.

Project coverage is 73.27%. Comparing base (9aaa78b) to head (95632d0).
Report is 9 commits behind head on main.

Files with missing lines	Patch %	Lines
lstchain/reco/dl1_to_dl2.py	20.33%	47 Missing ⚠️
lstchain/scripts/lstchain_dl1_to_dl2.py	69.09%	17 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1320      +/-   ##
==========================================
- Coverage   73.52%   73.27%   -0.25%     
==========================================
  Files         134      134              
  Lines       14215    14321     +106     
==========================================
+ Hits        10451    10494      +43     
- Misses       3764     3827      +63

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

command-line options The interpolation can be switched on for the three reconstructions independently (energy, gammaness, direction) By default the correction is applied. If no MC DL1 training directory is provided, the path will be built from the path of the RF models. If directory does not exist the option is simply deactivated

vuillaut · 2024-11-28T12:16:45Z

Hi @moralejo
Thank you for making this implementation.
Early feedback:
would it be possible to add the training and testing nodes of the current production as an ascii file in the lstchain data files to allow for testing without necessarily being on the cluster?
(astropy tables can be dumped to ascii easily)
I completely agree these pointing nodes should be saved with the models during training and then loaded from there in the future.

moralejo · 2024-11-28T13:37:08Z

Hi @moralejo Thank you for making this implementation. Early feedback: would it be possible to add the training and testing nodes of the current production as an ascii file in the lstchain data files to allow for testing without necessarily being on the cluster? (astropy tables can be dumped to ascii easily)

Will do. What is better, just the zenith and azimuth values, or a list of the directory names?

vuillaut · 2024-11-28T14:06:03Z

Hi @moralejo Thank you for making this implementation. Early feedback: would it be possible to add the training and testing nodes of the current production as an ascii file in the lstchain data files to allow for testing without necessarily being on the cluster? (astropy tables can be dumped to ascii easily)

Will do. What is better, just the zenith and azimuth values, or a list of the directory names?

I think the table of values that you can use directly at inference would be great (with appropriate metadata such as the date of export).

moralejo · 2024-11-28T14:19:05Z

Will do. What is better, just the zenith and azimuth values, or a list of the directory names?

I think the table of values that you can use directly at inference would be great (with appropriate metadata such as the date of export).

Can you clarify exactly what metadata (and how? Just lines starting with # or whatever?)

moralejo · 2024-11-28T14:24:15Z

What if we just write in the config the two arrays of zeniths and azimuths?
(considering it is just for local testing, not needed for normal execution)

I did not implement it yet, it would be like this:

"random_forest_zd_interpolation": {

"interpolate_energy": true,
"interpolate_gammaness": true,
"interpolate_direction": true,
"DL1_training_dir": null,
"DL1_training_pointings": [
  [ 69.19,   27.055],
  [290.81,   27.055],
  [ 63.197,  63.108],
  [ 69.537,  34.367],
  [ 66.706,  19.807],
  [293.294,  19.807],
  [ 67.397,  48.911],
  [296.803,  63.108],
  [292.603,  48.911],
  [ 65.502,  56.069],
  [ 58.737,  12.829],
  [290.463,  34.367],
  [ 60.505,  69.992],
  [328.91,    7.093],
  [ 68.802,  41.666],
  [291.198,  41.666],
  [294.498,  56.069],
  [301.263,  12.829],
  [ 31.09,    7.093],
  [299.495,  69.992]
]

}

(for RF interpolation)

Instead of closest angular distance we now take the closest nodes in alt_tel on the same side of the culmination (same sign of sin_az_tel)

Copilot reviewed 2 out of 3 changed files in this pull request and generated no suggestions.

Files not reviewed (1)

lstchain/data/lstchain_standard_config.json: Language not supported

Comments skipped due to low confidence (1)

lstchain/reco/dl1_to_dl2.py:215

The word 'dictionnary' is misspelled. It should be 'dictionary'.

config: dictionnary containing configuration

Copilot reviewed 2 out of 3 changed files in this pull request and generated no comments.

Files not reviewed (1)

lstchain/data/lstchain_standard_config.json: Language not supported

Comments suppressed due to low confidence (1)

lstchain/scripts/lstchain_dl1_to_dl2.py:164

Ensure that the training_pointings array is correctly populated and used. If the directory does not exist, the interpolation should be switched off.

if dl1_training_dir.is_dir():

vuillaut

Hi Abelardo,
Sorry for the delay (and I still need to finish this later...).
I have made a couple of comments for improvement in the code/readability but overall the core logic looks good to me.

vuillaut · 2025-01-13T15:23:34Z

lstchain/reco/dl1_to_dl2.py

+    DL2 pandas dataframe with additional columns alt0, alt1, w0, w1
+
+    """
+    pd.options.mode.copy_on_write = True


This will change the behaviour globally.
I'd suggest either to move it at the module level if that is wanted. Otherwise use context management to set it locally:

with pd.option_context('mode.copy_on_write', True): ...

vuillaut · 2025-01-13T15:24:53Z

lstchain/reco/dl1_to_dl2.py

@@ -51,6 +53,159 @@
    'update_disp_with_effective_focal_length'
 ]

+def add_zd_interpolation_info(dl2table, training_zd_deg, training_az_deg):


wouldn't using astropy units for quantities be better here?

vuillaut · 2025-01-13T15:25:45Z

lstchain/reco/dl1_to_dl2.py

+    alt_tel = dl2table['alt_tel']
+    az_tel  = dl2table['az_tel']
+
+    training_alt_rad = np.pi / 2 - np.deg2rad(training_zd_deg)


continuing comment about astropy quantities: it would avoid manual conversion

vuillaut · 2025-01-13T15:38:45Z

lstchain/reco/dl1_to_dl2.py

+    training_az_rad = np.deg2rad(training_az_deg)
+
+    tiled_az = np.tile(az_tel,
+                       len(training_az_rad)).reshape(len(training_az_rad),


tiled_az = np.broadcast_to(az_tel[:, np.newaxis], (len(dl2table), len(training_az_rad)))

would be more memory efficient as it would not create a copy of the array

vuillaut · 2025-01-13T16:07:20Z

lstchain/reco/dl1_to_dl2.py

+                  0.5, (cos_tel_zd - c0) / (c1 - c0))
+    w0 = 1 - w1
+
+    # Update the dataframe:


I don't see why converting to Series is necessary?

vuillaut · 2025-01-13T16:10:41Z

lstchain/reco/dl1_to_dl2.py

+
+    # Update the dataframe:
+    dl2table = dl2table.assign(alt0=pd.Series(closest_alt).values,
+                               alt1=pd.Series(second_closest_alt).values,


please use more explicit names for alt0, alt1, c0 and c1 in the dl2table that will end up being written

vuillaut · 2025-01-13T16:25:18Z

lstchain/reco/dl1_to_dl2.py

+    is_classifier = isinstance(rf, RandomForestClassifier)
+
+    # keep original alt_tel values:
+    param_array.rename(columns={"alt_tel": "original_alt_tel"}, inplace=True)


I think it would be easier to change the features rather than this game of name-swapping:

features_copy = features.copy() # avoid modifying features outside the function alt_index_in_features = features_copy.index('alt_tel') features_copy[alt_index_in_features] = 'alt0' predict_0 = rf.predict(param_array[features_copy]) features_copy[alt_index_in_features] = 'alt1' predict_1 = ...

that's it

vuillaut · 2025-01-13T16:34:00Z

lstchain/scripts/lstchain_dl1_to_dl2.py

+        c1 = dir.find('_theta_') + 7
+        c2 = dir.find('_az_', c1) + 4
+        c3 = dir.find('_', c2)


# Regular expression to extract theta and az values pattern = r"theta_(\d+\.\d+)_az_(\d+\.\d+)_" # Search for matches match = re.search(pattern, name) if match: theta = float(match.group(1)) # First captured group is theta az = float(match.group(2)) # Second captured group is az print(f"Theta: {theta}, Az: {az}") else: print("Pattern not found!")

would be more robust to find digital numbers

This function could even go to lstchain.paths module

vuillaut · 2025-01-13T16:38:12Z

lstchain/scripts/lstchain_dl1_to_dl2.py

+    # by the discrete set of pointings in the RF training sample.
+    if 'random_forest_zd_interpolation' in config.keys():
+        zdinter = config['random_forest_zd_interpolation']
+        if 'interpolate_energy' in zdinter.keys():


slightly more readable:

interpolate_energy = zdinter.get('interpolate_energy', False)

and then you don't need to set it first, it is set to False here if interpolate_energy not in zdinter

vuillaut · 2025-01-14T09:04:22Z

ely agree these pointing nodes should

Hi

I think we should go beyond testing actually.

Ideally, the training pointing directions should be saved with the models, and this PR should prepare for that imho.
However, in the meantime, to allow for using current models, the current training directions should be known.
I see two possibilities:

the training directions are saved as ascii files, produced automatically with the models. At inference this ascii file is loaded with the model to get training direction. But we could actually produce these ascii files manually for the current productions and add them to the models directories in the cluster to allow for interpolated inference already.
the training directions are saved within the pickle files with the trained models (Advantage: the information can't be separated from the model. Drawback: model loading is broken for past models) and loaded for inference. In the meantime, we should rely on loading the pointings from elsewhere (I'd say a config file within lstchain) - Thinking about it I think this 2nd solution actually has more drawbacks. Solution 1 could be implemented already and would avoid all the name guessing game at inference, thus simplifying the config.

To save the pointing directions as ascii files, I'd recommend starting from astropy tables and save to ecsv

vuillaut · 2025-01-14T08:23:21Z

lstchain/data/lstchain_standard_config.json

+    "interpolate_energy": true,
+    "interpolate_gammaness": true,
+    "interpolate_direction": true,
+    "DL1_training_dir": null


see my comment in the main discussion about saving training directions with the models

vuillaut · 2025-01-14T08:25:05Z

lstchain/scripts/lstchain_dl1_to_dl2.py

+        if 'random_forest_zd_interpolation' in config.keys():
+            zdinter = config['random_forest_zd_interpolation']
+            if 'dl1_training_dir' in zdinter.keys():
+                dl1_training_dir = zdinter['dl1_training_dir']


dl1_training_dir = zdinter.get('dl1_trainin_dir', None)

morcuended

Left some minor comments.

morcuended · 2025-01-14T11:25:24Z

lstchain/scripts/lstchain_dl1_to_dl2.py

+        c1 = dir.find('_theta_') + 7
+        c2 = dir.find('_az_', c1) + 4
+        c3 = dir.find('_', c2)


This function could even go to lstchain.paths module

morcuended · 2025-01-14T11:29:20Z

lstchain/scripts/lstchain_dl1_to_dl2.py

+            logger.warning('RF training pointings (az_deg, zd_deg):')
+            logger.warning(training_pointings)


log.info() here, it is the expected behavior

morcuended · 2025-01-14T12:06:33Z

lstchain/scripts/lstchain_dl1_to_dl2.py

+            # Build name of DL1 MC training files assuming standard pattern:
+            dummy = models_path.as_posix().replace('/data/models', '/data/mc/DL1')
+            dl1_training_dir = (Path(dummy[:dummy.rfind('/dec')] +
+                                     '/TrainingDataset/GammaDiffuse/' +
+                                     dummy[dummy.rfind('/dec') + 1:]))


I think this should go to lstchain.paths module.

morcuended · 2025-01-14T12:07:08Z

lstchain/scripts/lstchain_dl1_to_dl2.py

+        if interpolate_direction:
+            logger.warning('   direction reconstruction Random Forest')
+
+        if 'random_forest_zd_interpolation' in config.keys():


Suggested change

if 'random_forest_zd_interpolation' in config.keys():

if 'random_forest_zd_interpolation' in config:

morcuended · 2025-01-14T12:09:22Z

lstchain/scripts/lstchain_dl1_to_dl2.py

+    dl1_training_dir = None
+    training_pointings = None
+    if True in interpolate_rf.values():
+        logger.warning('Cos(zenith) interpolation will be used in:')


morcuended · 2025-01-14T13:42:37Z

Can you clarify exactly what metadata (and how? Just lines starting with # or whatever?)

Regarding metadata, if you're using astropy Tables to later write ECSV files, you can define the metadata as:

pointing_table = Table(...)
pointing_table.meta["export_date"] = datetime.today().strftime('%Y-%m-%d')
pointing_table.write("filepath_pointings.ecsv", format="ascii.ecsv")

moralejo added 2 commits November 27, 2024 10:47

New feature to perform interpolation of RF predictions vs. cos(zenith)

a2ef2f0

to avoid jumps when using RFs trained on a discrete set of pointings

Improved help message

908e6cf

moralejo mentioned this pull request Nov 27, 2024

Make RF response more homogeneous vs. zenith #1317

Open

moralejo added 11 commits November 27, 2024 16:04

Fixed various errors...

000f091

Add is_proba option for RF classifiers

19c6480

Proper treatment of predict_proba output

dba62a2

minor fix of matrix operations

d45ef90

Check RF type automatically - remove argument

d4a28db

Improved description

4b1c960

Change in treatment of missing dl1_training_dir case

309ab49

Fix building of dl1_training_dir from path to RF models

98a90ae

Typo

2acdf9f

Protection against MC training nodes with identical zenith

588a248

moralejo added 5 commits November 28, 2024 15:53

Preparations for better configuration of training pointings

2969a25

(for RF interpolation)

Moved function

32620d9

import glob

a65112b

Wrong variable scope!

9e37b7a

Better log message

80c0ffc

moralejo marked this pull request as ready for review November 29, 2024 14:58

Changed selection of pointing nodes for interpolation.

b0afdac

Instead of closest angular distance we now take the closest nodes in alt_tel on the same side of the culmination (same sign of sin_az_tel)

moralejo requested a review from vuillaut December 3, 2024 14:16

vuillaut requested a review from Copilot December 11, 2024 14:26

Copilot AI reviewed Dec 11, 2024

View reviewed changes

typo

95632d0

vuillaut requested a review from Copilot December 19, 2024 10:05

Copilot AI reviewed Dec 19, 2024

View reviewed changes

vuillaut requested changes Jan 13, 2025

View reviewed changes

vuillaut requested changes Jan 14, 2025

View reviewed changes

morcuended requested changes Jan 14, 2025

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Interpolation of RF predictions with cosZD, for homogeneous performance #1320

Interpolation of RF predictions with cosZD, for homogeneous performance #1320

moralejo commented Nov 27, 2024 •

edited

Loading

codecov bot commented Nov 27, 2024 •

edited

Loading

vuillaut commented Nov 28, 2024

moralejo commented Nov 28, 2024

vuillaut commented Nov 28, 2024

moralejo commented Nov 28, 2024

moralejo commented Nov 28, 2024 •

edited

Loading

vuillaut left a comment

vuillaut Jan 13, 2025

vuillaut Jan 13, 2025

vuillaut Jan 13, 2025

vuillaut Jan 13, 2025

vuillaut Jan 13, 2025

vuillaut Jan 13, 2025

vuillaut Jan 13, 2025

vuillaut Jan 13, 2025

morcuended Jan 14, 2025

vuillaut Jan 13, 2025

vuillaut commented Jan 14, 2025 •

edited

Loading

vuillaut Jan 14, 2025

vuillaut Jan 14, 2025

morcuended left a comment

morcuended Jan 14, 2025

morcuended Jan 14, 2025

morcuended Jan 14, 2025

morcuended Jan 14, 2025

morcuended Jan 14, 2025

morcuended commented Jan 14, 2025

		logger.warning('RF training pointings (az_deg, zd_deg):')
		logger.warning(training_pointings)

	if 'random_forest_zd_interpolation' in config.keys():
	if 'random_forest_zd_interpolation' in config:

Interpolation of RF predictions with cosZD, for homogeneous performance #1320

Are you sure you want to change the base?

Interpolation of RF predictions with cosZD, for homogeneous performance #1320

Conversation

moralejo commented Nov 27, 2024 • edited Loading

codecov bot commented Nov 27, 2024 • edited Loading

Codecov Report

vuillaut commented Nov 28, 2024

moralejo commented Nov 28, 2024

vuillaut commented Nov 28, 2024

moralejo commented Nov 28, 2024

moralejo commented Nov 28, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

vuillaut left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

vuillaut commented Jan 14, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

morcuended left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

morcuended commented Jan 14, 2025

moralejo commented Nov 27, 2024 •

edited

Loading

codecov bot commented Nov 27, 2024 •

edited

Loading

moralejo commented Nov 28, 2024 •

edited

Loading

vuillaut commented Jan 14, 2025 •

edited

Loading