Merge pull request cms-sw#40779 from aloeliger/anomalyDetectionFromExternal

"CICADA"/calorimeter anomaly detection trigger emulator using external repositories

Author: cmsbuild

L1Trigger/L1TCaloLayer1/plugins/BuildFile.xml

@@ -5,4 +5,7 @@
 <use name="DataFormats/L1TCalorimeter"/>
 <use name="L1Trigger/L1TCalorimeter"/>
 <use name="L1Trigger/L1TCaloLayer1"/>
+<use name="hls"/>
+<use name="hls4mlEmulatorExtras"/>
+<use name="CICADA"/>
 <flags EDM_PLUGIN="1"/>

L1Trigger/L1TCaloLayer1/plugins/L1TCaloSummary.cc

@@ -57,6 +57,10 @@
 #include "L1Trigger/L1TCaloLayer1/src/UCTLogging.hh"
 #include <bitset>
 
+//Anomaly detection includes
+#include "ap_fixed.h"
+#include "hls4ml/emulator.h"
+
 using namespace l1tcalo;
 using namespace l1extra;
 using namespace std;
@@ -102,6 +106,9 @@ class L1TCaloSummary : public edm::stream::EDProducer<> {
   edm::EDGetTokenT<L1CaloRegionCollection> regionToken;
 
   UCTLayer1* layer1;
+
+  hls4mlEmulator::ModelLoader loader;
+  std::shared_ptr<hls4mlEmulator::Model> model;
 };
 
 //
@@ -127,7 +134,8 @@ L1TCaloSummary::L1TCaloSummary(const edm::ParameterSet& iConfig)
       boostedJetPtFactor(iConfig.getParameter<double>("boostedJetPtFactor")),
       verbose(iConfig.getParameter<bool>("verbose")),
       fwVersion(iConfig.getParameter<int>("firmwareVersion")),
-      regionToken(consumes<L1CaloRegionCollection>(edm::InputTag("simCaloStage2Layer1Digis"))) {
+      regionToken(consumes<L1CaloRegionCollection>(edm::InputTag("simCaloStage2Layer1Digis"))),
+      loader(hls4mlEmulator::ModelLoader(iConfig.getParameter<string>("CICADAModelVersion"))) {
   std::vector<double> pumLUTData;
   char pumLUTString[10];
   for (uint32_t pumBin = 0; pumBin < nPumBins; pumBin++) {
@@ -147,6 +155,10 @@ L1TCaloSummary::L1TCaloSummary(const edm::ParameterSet& iConfig)
     }
   }
   produces<L1JetParticleCollection>("Boosted");
+
+  //anomaly trigger loading
+  model = loader.load_model();
+  produces<float>("anomalyScore");
 }
 
 L1TCaloSummary::~L1TCaloSummary() {}
@@ -161,6 +173,8 @@ void L1TCaloSummary::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
 
   std::unique_ptr<L1JetParticleCollection> bJetCands(new L1JetParticleCollection);
 
+  std::unique_ptr<float> anomalyScore = std::make_unique<float>();
+
   UCTGeometry g;
 
   // Here we read region data from the region collection created by L1TCaloLayer1 instead of
@@ -175,6 +189,11 @@ void L1TCaloSummary::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
   if (!iEvent.getByToken(regionToken, regionCollection))
     edm::LogError("L1TCaloSummary") << "UCT: Failed to get regions from region collection!";
   iEvent.getByToken(regionToken, regionCollection);
+  //Model input
+  //This is done as a flat vector input, but future versions may involve 2D input
+  //This will have to be handled later
+  //Would also be good to be able to configure the precision of the ap_fixed type
+  ap_ufixed<10, 10> modelInput[252];
   for (const L1CaloRegion& i : *regionCollection) {
     UCTRegionIndex r = g.getUCTRegionIndexFromL1CaloRegion(i.gctEta(), i.gctPhi());
     UCTTowerIndex t = g.getUCTTowerIndexFromL1CaloRegion(r, i.raw());
@@ -189,7 +208,22 @@ void L1TCaloSummary::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
     UCTRegion* test = new UCTRegion(crate, card, negativeEta, region, fwVersion);
     test->setRegionSummary(i.raw());
     inputRegions.push_back(test);
+    //This *should* fill the tensor in the proper order to be fed to the anomaly model
+    //We take 4 off of the GCT eta/iEta.
+    //iEta taken from this ranges from 4-17, (I assume reserving lower and higher for forward regions)
+    //So our first index, index 0, is technically iEta=4, and so-on.
+    //CICADA v1 reads this as a flat vector
+    modelInput[14 * i.gctPhi() + (i.gctEta() - 4)] = i.et();
   }
+  //Extract model output
+  //Would be good to be able to configure the precision of the result
+  ap_fixed<11, 5> modelResult[1];
+  model->prepare_input(modelInput);
+  model->predict();
+  model->read_result(modelResult);
+
+  *anomalyScore = modelResult[0].to_float();
+
   summaryCard.setRegionData(inputRegions);
 
   if (!summaryCard.process()) {
@@ -253,6 +287,8 @@ void L1TCaloSummary::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
   }
 
   iEvent.put(std::move(bJetCands), "Boosted");
+  //Write out anomaly score
+  iEvent.put(std::move(anomalyScore), "anomalyScore");
 }
 
 void L1TCaloSummary::print() {}

L1Trigger/L1TCaloLayer1/python/uct2016EmulatorDigis_cfi.py

@@ -47,5 +47,6 @@
                                       boostedJetPtFactor = cms.double(1.5),
                                       verbose = cms.bool(False),
                                       # See UCTLayer1.hh for firmware version
-                                      firmwareVersion = cms.int32(1)
+                                      firmwareVersion = cms.int32(1),
+                                      CICADAModelVersion = cms.string("CICADAModel_v1")
                                       )

L1Trigger/L1TCaloLayer1/test/BuildFile.xml

@@ -7,3 +7,10 @@
 
 <bin name="testUCTLayer1HF" file="testUCTLayer1HF.cpp">
 </bin>
+
+<bin name="testCICADAEmulation" file="testCICADAEmulation.cppunit.cc">
+    <use name="hls"/>
+    <use name="hls4mlEmulatorExtras"/>
+    <use name="CICADA"/>
+    <use name="cppunit"/>
+</bin>

L1Trigger/L1TCaloLayer1/test/testCICADAEmulation.cppunit.cc

@@ -0,0 +1,52 @@
+//Test of the external CICADA model emulation model loading and model unloading
+//Developed by Andrew Loeliger, Princeton University, Feb 23, 2023
+
+//We can't test a load of a bad model here, since that is a segfault, not an exception, which is
+//OS level and cppunit cannot test against that in any way that qualifies as a success
+
+//TODO: However, it would be good in the future to assure that loading multiple CICADA models at the
+//same time have the correct function symbols assigned to each simultaneously
+//i.e. CICADA_v1's predict is not overwritten by CICADA_v2's predict if it is loaded later with a
+//CICADA_v1 still around
+
+//TODO: might also be nice to have a test for model integrity? Known test cases producing known outputs?
+//This may not be appropriate for unit testing however.
+
+#include "ap_fixed.h"
+#include "hls4ml/emulator.h"
+
+#include "cppunit/extensions/HelperMacros.h"
+#include <memory>
+#include "Utilities/Testing/interface/CppUnit_testdriver.icpp"
+
+class test_CICADA : public CppUnit::TestFixture {
+  CPPUNIT_TEST_SUITE(test_CICADA);
+  CPPUNIT_TEST(doModelV1Load);
+  CPPUNIT_TEST(doModelV2Load);
+  CPPUNIT_TEST(doMultiModelLoad);
+  CPPUNIT_TEST_SUITE_END();
+
+public:
+  void doModelV1Load();
+  void doModelV2Load();
+  void doMultiModelLoad();
+};
+
+CPPUNIT_TEST_SUITE_REGISTRATION(test_CICADA);
+
+void test_CICADA::doModelV1Load() {
+  auto loader = hls4mlEmulator::ModelLoader("CICADAModel_v1");
+  auto model = loader.load_model();
+}
+
+void test_CICADA::doModelV2Load() {
+  auto loader = hls4mlEmulator::ModelLoader("CICADAModel_v2");
+  auto model = loader.load_model();
+}
+
+void test_CICADA::doMultiModelLoad() {
+  auto loader_v1 = hls4mlEmulator::ModelLoader("CICADAModel_v1");
+  auto loader_v2 = hls4mlEmulator::ModelLoader("CICADAModel_v2");
+  auto model_v1 = loader_v1.load_model();
+  auto model_v2 = loader_v2.load_model();
+}