Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329
/*
 * OnlineDQMDigiAD_cmssw.cpp
 *
 * Created on: Jun 10, 2023
 * Author: Mulugeta W.Asres, UiA, Norway
 *
 * The implementation follows https://github.com/cms-sw/cmssw/tree/master/PhysicsTools/ONNXRuntime
 */

// #include "FWCore/Utilities/interface/Exception.h"
// #include "FWCore/Utilities/interface/thread_safety_macros.h"
// #include "FWCore/Framework/interface/Event.h"
#include "PhysicsTools/ONNXRuntime/interface/ONNXRuntime.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"
#include "HeterogeneousCore/CUDAUtilities/interface/requireDevices.h"

#include <algorithm>
#include <cassert>
#include <functional>
#include <iostream>
#include <memory>
#include <numeric>
#include <algorithm>

#include "DQM/HcalTasks/plugins/OnlineDQMDigiAD_cmssw.h"

// using namespace std;
using namespace cms::Ort;

// Constructor
OnlineDQMDigiAD::OnlineDQMDigiAD(const std::string model_system_name,
                                 const std::string &modelFilepath,
                                 Backend backend) {
  std::string instanceName{"DESMOD Digioccupancy Map AD inference"};

  /**************** Initailize Model Memory States ******************/
  InitializeState();  // initailize model memory states to zero

  /**************** Create ORT session ******************/
  // Set up options for session
  auto session_options = ONNXRuntime::defaultSessionOptions(backend);
  // Create session by loading the onnx model
  model_path = edm::FileInPath(modelFilepath).fullPath();
  auto uOrtSession = std::make_unique<ONNXRuntime>(model_path, &session_options);
  ort_mSession = std::move(uOrtSession);

  // check model availability
  hcal_subsystem_name = model_system_name;

  IsModelExist(hcal_subsystem_name);  // assert model integration for the given hcal system name

  if (hcal_subsystem_name == "he") {
    std::vector<std::vector<int64_t>> input_shapes_ = {
        {batch_size, 64, 72, 7, 1},
        {batch_size, 1},
        {1, 1},
        {batch_size, model_state_inner_dim, model_state_layer_dims[0][0]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[0][0]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[0][1]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[0][1]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[1][0]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[1][0]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[1][1]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[1][1]}};  // input dims
    input_shapes = input_shapes_;
  }

  else if (hcal_subsystem_name == "hb") {
    std::vector<std::vector<int64_t>> input_shapes_ = {
        {batch_size, 64, 72, 4, 1},
        {batch_size, 1},
        {1, 1},
        {batch_size, model_state_inner_dim, model_state_layer_dims[0][0]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[0][0]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[0][1]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[0][1]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[1][0]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[1][0]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[1][1]},
        {batch_size, model_state_inner_dim, model_state_layer_dims[1][1]}};  // input dims
    input_shapes = input_shapes_;
  }
}

void OnlineDQMDigiAD::IsModelExist(std::string hcal_subsystem_name) {
  if (std::find(hcal_modeled_systems.begin(), hcal_modeled_systems.end(), hcal_subsystem_name) ==
      hcal_modeled_systems.end()) {
    std::string err =
        "ML for OnlineDQM is not currently supported for the selected " + hcal_subsystem_name + " system!\n";
    throw std::invalid_argument(err);
  }
}

void OnlineDQMDigiAD::InitializeState() {
  // model memory states vectors init, only when the runs starts or for the first LS
  std::fill(input_model_state_memory_e_0_0.begin(),
            input_model_state_memory_e_0_0.end(),
            float(0.0));  // init model memory states-encoder_layer_0_state_0 to zero
  std::fill(input_model_state_memory_e_0_1.begin(),
            input_model_state_memory_e_0_1.end(),
            float(0.0));  // init model memory states-encoder_layer_0_state_1 to zero
  std::fill(input_model_state_memory_e_1_0.begin(),
            input_model_state_memory_e_1_0.end(),
            float(0.0));  // init model memory states-encoder_layer_1_state_0 to zero
  std::fill(input_model_state_memory_e_1_1.begin(),
            input_model_state_memory_e_1_1.end(),
            float(0.0));  // init model memory states-encoder_layer_1_state_1 to zero
  std::fill(input_model_state_memory_d_0_0.begin(),
            input_model_state_memory_d_0_0.end(),
            float(0.0));  // init model memory states-decoder_layer_0_state_0 to zero
  std::fill(input_model_state_memory_d_0_1.begin(),
            input_model_state_memory_d_0_1.end(),
            float(0.0));  // init model memory states-decoder_layer_0_state_1 to zero
  std::fill(input_model_state_memory_d_1_0.begin(),
            input_model_state_memory_d_1_0.end(),
            float(0.0));  // init model memory states-decoder_layer_1_state_0 to zero
  std::fill(input_model_state_memory_d_1_1.begin(),
            input_model_state_memory_d_1_1.end(),
            float(0.0));  // init model memory states-decoder_layer_1_state_1 to zero

  // model_state_refresh_counter = 15; // counter set due to onnx double datatype handling limitation that might cause precision error to propagate.
  model_state_refresh_counter =
      1;  // DQM multithread returns non-sequential LS. Hence, the model will not keep states (experimental)
}

std::vector<float> OnlineDQMDigiAD::Serialize2DVector(const std::vector<std::vector<float>> &input_2d_vec) {
  std::vector<float> output;
  for (const auto &row : input_2d_vec) {
    for (const auto &element : row) {
      output.push_back(element);
    }
  }
  return output;
}

std::vector<std::vector<float>> OnlineDQMDigiAD::Map1DTo2DVector(const std::vector<float> &input_1d_vec,
                                                                 const int numSplits) {
  if (numSplits <= 0)
    throw std::invalid_argument("numSplits must be greater than 0.");

  std::size_t const splitted_size = input_1d_vec.size() / numSplits;

  if (splitted_size * numSplits != input_1d_vec.size())
    throw std::invalid_argument("Conversion is not allowed! The input vector length " +
                                std::to_string(input_1d_vec.size()) + " must be divisible by the numSplits " +
                                std::to_string(numSplits) + ".");

  std::vector<std::vector<float>> output_2d_vec;

  for (int i = 0; i < numSplits; i++) {
    std::vector<float> chunch_vec(input_1d_vec.begin() + i * splitted_size,
                                  input_1d_vec.begin() + (i + 1) * splitted_size);
    output_2d_vec.push_back(chunch_vec);
  }
  return output_2d_vec;
}

std::vector<float> OnlineDQMDigiAD::PrepareONNXDQMMapVectors(
    std::vector<std::vector<std::vector<float>>> &digiHcal2DHist_depth_all) {
  std::vector<float> digi3DHistVector_serialized;

  for (const std::vector<std::vector<float>> &digiHcal2DHist_depth : digiHcal2DHist_depth_all) {
    std::vector<float> digiHcalDHist_serialized_depth = Serialize2DVector(digiHcal2DHist_depth);
    digi3DHistVector_serialized.insert(digi3DHistVector_serialized.end(),
                                       digiHcalDHist_serialized_depth.begin(),
                                       digiHcalDHist_serialized_depth.end());
  }

  return digi3DHistVector_serialized;
}

std::vector<std::vector<std::vector<float>>> OnlineDQMDigiAD::ONNXOutputToDQMHistMap(
    const std::vector<std::vector<float>> &ad_model_output_vectors,
    const int numDepth,
    const int numDIeta,
    const int selOutputIdx) {
  // each output_vector is a serialized 3d hist map

  const std::vector<float> &output_vector = ad_model_output_vectors[selOutputIdx];
  std::vector<std::vector<float>> output_2d_vec = Map1DTo2DVector(output_vector, numDepth);

  std::vector<std::vector<std::vector<float>>> digiHcal3DHist;
  for (const std::vector<float> &output_vector_depth : output_2d_vec) {
    std::vector<std::vector<float>> digiHcal2DHist_depth = Map1DTo2DVector(output_vector_depth, numDIeta);
    digiHcal3DHist.push_back(digiHcal2DHist_depth);
  }

  return digiHcal3DHist;
}

// Perform inference for a given dqm map
std::vector<std::vector<float>> OnlineDQMDigiAD::Inference(std::vector<float> &digiHcalMapTW,
                                                           std::vector<float> &numEvents,
                                                           std::vector<float> &adThr,
                                                           std::vector<float> &input_model_state_memory_e_0_0,
                                                           std::vector<float> &input_model_state_memory_e_0_1,
                                                           std::vector<float> &input_model_state_memory_e_1_0,
                                                           std::vector<float> &input_model_state_memory_e_1_1,
                                                           std::vector<float> &input_model_state_memory_d_0_0,
                                                           std::vector<float> &input_model_state_memory_d_0_1,
                                                           std::vector<float> &input_model_state_memory_d_1_0,
                                                           std::vector<float> &input_model_state_memory_d_1_1) {
  /**************** Preprocessing ******************/
  // Create input tensor (including size and value) from the loaded inputs
  // Compute the product of all input dimension
  // Assign memory for input tensor
  // inputTensors will be used by the Session Run for inference

  input_values.clear();
  input_values.emplace_back(digiHcalMapTW);
  input_values.emplace_back(numEvents);
  input_values.emplace_back(adThr);
  input_values.emplace_back(input_model_state_memory_e_0_0);
  input_values.emplace_back(input_model_state_memory_e_0_1);
  input_values.emplace_back(input_model_state_memory_e_1_0);
  input_values.emplace_back(input_model_state_memory_e_1_1);
  input_values.emplace_back(input_model_state_memory_d_0_0);
  input_values.emplace_back(input_model_state_memory_d_0_1);
  input_values.emplace_back(input_model_state_memory_d_1_0);
  input_values.emplace_back(input_model_state_memory_d_1_1);

  /**************** Inference ******************/

  output_values = ort_mSession->run(input_names, input_values, input_shapes, output_names, batch_size);

  return output_values;
}

// AD method to be called by the CMS system
std::vector<std::vector<float>> OnlineDQMDigiAD::Inference_CMSSW(
    const std::vector<std::vector<float>> &digiHcal2DHist_depth_1,
    const std::vector<std::vector<float>> &digiHcal2DHist_depth_2,
    const std::vector<std::vector<float>> &digiHcal2DHist_depth_3,
    const std::vector<std::vector<float>> &digiHcal2DHist_depth_4,
    const std::vector<std::vector<float>> &digiHcal2DHist_depth_5,
    const std::vector<std::vector<float>> &digiHcal2DHist_depth_6,
    const std::vector<std::vector<float>> &digiHcal2DHist_depth_7,
    const float LS_numEvents,
    const float flagDecisionThr)

{
  /**************** Prepare data ******************/
  // merging all 2d hist into one 3d depth[ieta[iphi]]
  std::vector<std::vector<std::vector<float>>> digiHcal2DHist_depth_all;

  if (hcal_subsystem_name == "he") {
    digiHcal2DHist_depth_all.push_back(digiHcal2DHist_depth_1);
    digiHcal2DHist_depth_all.push_back(digiHcal2DHist_depth_2);
    digiHcal2DHist_depth_all.push_back(digiHcal2DHist_depth_3);
    digiHcal2DHist_depth_all.push_back(digiHcal2DHist_depth_4);
    digiHcal2DHist_depth_all.push_back(digiHcal2DHist_depth_5);
    digiHcal2DHist_depth_all.push_back(digiHcal2DHist_depth_6);
    digiHcal2DHist_depth_all.push_back(digiHcal2DHist_depth_7);
  }

  else if (hcal_subsystem_name == "hb") {
    digiHcal2DHist_depth_all.push_back(digiHcal2DHist_depth_1);
    digiHcal2DHist_depth_all.push_back(digiHcal2DHist_depth_2);
    digiHcal2DHist_depth_all.push_back(digiHcal2DHist_depth_3);
    digiHcal2DHist_depth_all.push_back(digiHcal2DHist_depth_4);
  }

  // convert the 3d depth[ieta[iphi]] vector into 1d and commbined
  std::vector<float> digiHcalMapTW = PrepareONNXDQMMapVectors(digiHcal2DHist_depth_all);

  std::vector<float> adThr{flagDecisionThr};  // AD decision threshold, increase to reduce sensitivity
  std::vector<float> numEvents{LS_numEvents};

  // call model inference
  /**************** Inference ******************/
  std::vector<std::vector<float>> output_tensors = Inference(digiHcalMapTW,
                                                             numEvents,
                                                             adThr,
                                                             input_model_state_memory_e_0_0,
                                                             input_model_state_memory_e_0_1,
                                                             input_model_state_memory_e_1_0,
                                                             input_model_state_memory_e_1_1,
                                                             input_model_state_memory_d_0_0,
                                                             input_model_state_memory_d_0_1,
                                                             input_model_state_memory_d_1_0,
                                                             input_model_state_memory_d_1_1);

  // auto output_tensors = Inference(digiHcalMapTW, numEvents, adThr);
  //std::cout << "******* model inference is success *******" << std::endl;

  /**************** Output post processing ******************/
  //  split outputs into ad output vectors and state_memory vectors
  std::string state_output_name_tag = "rnn_hidden";
  std::vector<std::vector<float>> ad_model_output_vectors, ad_model_state_vectors;
  for (size_t i = 0; i < output_tensors.size(); i++) {
    std::string output_names_startstr = output_names[i].substr(
        2, state_output_name_tag.length());  // Extract the same number of characters as str2 from mOutputNames
    if (output_names_startstr == state_output_name_tag) {
      ad_model_state_vectors.emplace_back(output_tensors[i]);
    } else {
      ad_model_output_vectors.emplace_back(output_tensors[i]);
    }
  }

  if (ad_model_output_vectors.size() == num_state_vectors) {
    input_model_state_memory_e_0_0 = ad_model_state_vectors[0];
    input_model_state_memory_e_0_1 = ad_model_state_vectors[1];
    input_model_state_memory_e_1_0 = ad_model_state_vectors[2];
    input_model_state_memory_e_1_1 = ad_model_state_vectors[3];
    input_model_state_memory_d_0_0 = ad_model_state_vectors[4];
    input_model_state_memory_d_0_1 = ad_model_state_vectors[5];
    input_model_state_memory_d_1_0 = ad_model_state_vectors[6];
    input_model_state_memory_d_1_1 = ad_model_state_vectors[7];
  } else {
    std::cout << "Warning: the number of output state vectors does NOT equals to expected!. The states are set to  "
                 "default values."
              << std::endl;
    InitializeState();
  }

  // # if onnx is returning serialized 1d vectors instead of vector of 3d vectors
  // aml score and flag are at index 5 and 7 of the vector ad_model_output_vectors: anomaly score: ad_model_output_vectors[5], anomaly flags: ad_model_output_vectors[7]
  /*
      selOutputIdx: index to select of the onnx output. e.g. 5 is the anomaly score and 7 is the anomaly flag (1 is with anomaly, 0 is healthy)
      std::vector<std::vector<std::vector<float>>> digiHcal3DHist_ANOMALY_FLAG = ONNXOutputToDQMHistMap(ad_model_output_vectors, 7)
      std::vector<std::vector<std::vector<float>>> digiHcal3DHist_ANOMALY_SCORE = ONNXOutputToDQMHistMap(ad_model_output_vectors, 5)
      */

  // reduce counter for each ls call. due to onnx double datatype handling limitation that might cause precision error to propagate.
  if (--model_state_refresh_counter == 0)
    InitializeState();

  return ad_model_output_vectors;
}