GBTS hit bidding and dropout (#1296)

* seed hit dropout and bidding

* seed ambiwq

* move configs to struct

* formating

* adding seed sorting

* better dropout logic

* better naming

* update names

* add dropout for low-eta if all seeds are incompatable

* formatting

* removing seed sorting test adding use_dropout flag

* formatting

* reverting testing change in seed fitting

* fixing fp63->fp32 warnings

---------

Co-authored-by: Mark <markbrettel@cern.ch>
Co-authored-by: Stephen Nicholas Swatman <stephen.nicholas.swatman@cern.ch>

Author: 3 people

core/include/traccc/gbts_seeding/gbts_seeding_config.hpp

@@ -125,6 +125,16 @@ struct gbts_seed_extraction_params {
     float max_z0 = 160.0f;
 };
 
+struct gbts_seed_ambi_params {
+    // sample multiple triplets when forming seeds to hedge against outliers
+    bool use_dropout = true;
+    // these curvatures are in 1/m
+    float dropout_dcurv_m = 0.007f;
+    float force_dropout_max_curv_m = 0.03f;
+    float best_hit_frac = 0.49f;
+    float tight_bid_cot_threshold = 1.0f;
+};
+
 struct gbts_seedfinder_config {
     bool setLinkingScheme(
         const std::vector<std::pair<int, std::vector<int>>>& binTables,
@@ -142,8 +152,8 @@ struct gbts_seedfinder_config {
 
     // tuned for 900 MeV pT cut and scaled by input minPt
     gbts_graph_building_params graph_building_params{};
-
     gbts_seed_extraction_params seed_extraction_params{};
+    gbts_seed_ambi_params seed_ambi_params{};
 
     // node making bin counts
     unsigned int n_eta_bins = 0;  // calculated from input layerInfo

device/cuda/src/gbts_seeding/gbts_seeding_algorithm.cu

@@ -105,6 +105,7 @@ struct gbts_ctx {
     int2* d_seed_proposals{};  // int quality and final mini_state_idx
     // first 32 bits are seed quality second 32 bits are seed_proposals index
     unsigned long long int* d_edge_bids{};
+    unsigned long long int* d_hit_bids{};
     // 0 as default/is real seed, 1 as maybe seed,
     //-1 as maybe fake seed, -2 as fake
     char* d_seed_ambiguity{};
@@ -823,9 +824,15 @@ gbts_seeding_algorithm::output_type gbts_seeding_algorithm::operator()(
 
     cudaStreamSynchronize(stream);
 
+    error = cudaGetLastError();
+    if (error != cudaSuccess) {
+        TRACCC_ERROR(
+            "seed extraction: CUDA error: " << cudaGetErrorString(error));
+        return {0, m_mr.main};
+    }
+
     cudaFree(ctx.d_levels);
     cudaFree(ctx.d_outgoing_paths);
-    cudaFree(ctx.d_reducedSP);
 
     if (nProps == 0) {
         return {0, m_mr.main};
@@ -865,26 +872,47 @@ gbts_seeding_algorithm::output_type gbts_seeding_algorithm::operator()(
     cudaFree(ctx.d_graph_building_params);
 
     // 8. convert to 3sp seeds and make output buffer
-
+    // allocate extra seed space for hit permutation
     edm::seed_collection::buffer output_seeds(
-        ctx.nSeeds, m_mr.main, vecmem::data::buffer_type::resizable);
-    m_copy.get().setup(output_seeds)->ignore();
+        2 * ctx.nSeeds, m_mr.main, vecmem::data::buffer_type::resizable);
+    m_copy.get().setup(output_seeds)->wait();
+
+    cudaMalloc(&ctx.d_hit_bids, ctx.nSp * sizeof(unsigned long long int));
+    cudaMemsetAsync(ctx.d_hit_bids, 0, ctx.nSp * sizeof(unsigned long long int),
+                    stream);
+
+    nThreads = 128;
+    nBlocks = 1 + (ctx.nSeeds - 1) / nThreads;
+
+    kernels::seeds_bid_for_hits<<<nBlocks, nThreads, 0, stream>>>(
+        ctx.d_output_graph, ctx.d_seed_proposals, ctx.d_path_store,
+        ctx.d_seed_ambiguity, ctx.d_hit_bids, nProps,
+        1 + 2 + m_config.max_num_neighbours);
 
     kernels::gbts_seed_conversion_kernel<<<nBlocks, nThreads, 0, stream>>>(
         ctx.d_seed_proposals, ctx.d_seed_ambiguity, ctx.d_path_store,
-        ctx.d_output_graph, output_seeds, nProps, m_config.max_num_neighbours);
+        ctx.d_output_graph, ctx.d_reducedSP, output_seeds, ctx.d_hit_bids,
+        nProps, m_config.max_num_neighbours,
+        m_config.seed_ambi_params.dropout_dcurv_m,
+        m_config.seed_ambi_params.force_dropout_max_curv_m,
+        m_config.seed_ambi_params.best_hit_frac,
+        m_config.seed_ambi_params.tight_bid_cot_threshold,
+        m_config.seed_ambi_params.use_dropout);
 
     cudaStreamSynchronize(stream);
 
+    ctx.nSeeds = m_copy.get().get_size(output_seeds);
+
+    cudaFree(ctx.d_reducedSP);
     cudaFree(ctx.d_output_graph);
     cudaFree(ctx.d_path_store);
     cudaFree(ctx.d_seed_proposals);
     cudaFree(ctx.d_seed_ambiguity);
+    cudaFree(ctx.d_hit_bids);
 
     error = cudaGetLastError();
-
     if (error != cudaSuccess) {
-        TRACCC_ERROR("seed-extracting kalman filter: CUDA error: "
+        TRACCC_ERROR("Seed ambiguity solving: CUDA error: "
                      << cudaGetErrorString(error));
         return {0, m_mr.main};
     }

device/cuda/src/gbts_seeding/kernels/GbtsGraphProcessingKernels.cuh

@@ -693,10 +693,91 @@ void __global__ seeds_rebid_for_edges(int2* d_path_store,
     }
 }
 
+void __global__ seeds_bid_for_hits(int* d_output_graph, int2* d_seed_proposals,
+                                   int2* d_path_store, char* d_seed_ambiguity,
+                                   unsigned long long int* d_hit_bids,
+                                   const unsigned int nProps, int edge_size) {
+
+    for (unsigned int prop_idx = threadIdx.x + blockDim.x * blockIdx.x;
+         prop_idx < nProps; prop_idx += gridDim.x * blockDim.x) {
+        if (d_seed_ambiguity[prop_idx] == -2) {
+            continue;
+        }
+        int2 prop = d_seed_proposals[prop_idx];
+        unsigned long long int seed_bid =
+            (static_cast<unsigned long long int>(prop.x) << 32) |
+            (static_cast<unsigned long long int>(prop_idx));
+
+        int2 path = make_int2(0, prop.y);
+        while (path.y >= 0) {
+            path = d_path_store[path.y];
+            int sp_idx = d_output_graph[traccc::device::gbts_consts::node1 +
+                                        edge_size * path.x];
+            atomicMax(&d_hit_bids[sp_idx], seed_bid);
+        }
+        int sp_idx = d_output_graph[traccc::device::gbts_consts::node2 +
+                                    edge_size * path.x];
+        atomicMax(&d_hit_bids[sp_idx], seed_bid);
+    }
+}
+
+inline __device__ float2 estimate_params(float4 sps[3]) {
+
+    // conformal mapping with the center at the middle spacepoint
+
+    float u[2], v[2];
+
+    const float x0 = sps[1].x;
+    const float y0 = sps[1].y;
+
+    const float r0 = sqrtf(x0 * x0 + y0 * y0);
+
+    const float cosA = x0 / r0;
+
+    const float sinA = y0 / r0;
+
+    for (unsigned int k = 0; k < 2; k++) {
+
+        int sp_idx = (k == 1) ? 2 : k;
+
+        const float dx = sps[sp_idx].x - x0;
+
+        const float dy = sps[sp_idx].y - y0;
+
+        const float r2_inv = 1.0f / (dx * dx + dy * dy);
+
+        const float xn = dx * cosA + dy * sinA;
+
+        const float yn = -dx * sinA + dy * cosA;
+
+        u[k] = xn * r2_inv;
+        v[k] = yn * r2_inv;
+    }
+
+    const float du = u[0] - u[1];
+    if (du == 0.0f) {
+        return make_float2(0.0f, 0.0f);
+    }
+
+    const float A = (v[0] - v[1]) / du;
+
+    const float B = v[1] - A * u[1];
+
+    // signed curvature in 1/m
+    const float curv = 1000.0f * B / sqrtf(1 + A * A);
+    const float cot_t = (sps[2].z - sps[1].z) /
+                        (sqrtf(sps[2].x * sps[2].x + sps[2].y * sps[2].y) - r0);
+    return make_float2(curv, cot_t);
+}
+
 void __global__ gbts_seed_conversion_kernel(
     int2* d_seed_proposals, char* d_seed_ambiguity, int2* d_path_store,
-    int* d_output_graph, edm::seed_collection::view output_seeds,
-    const unsigned int nProps, const unsigned int max_num_neighbours) {
+    int* d_output_graph, float4* d_sp_params,
+    edm::seed_collection::view output_seeds, unsigned long long int* d_hit_bids,
+    const unsigned int nProps, const unsigned int max_num_neighbours,
+    const float dcurv_cut_m, const float force_dropout_max_curv_m,
+    const float best_hit_frac, const float tight_bid_cot_threshold,
+    const bool use_dropout) {
 
     int edge_size = 2 + 1 + max_num_neighbours;
     edm::seed_collection::device seeds_device(output_seeds);
@@ -707,23 +788,88 @@ void __global__ gbts_seed_conversion_kernel(
             // drop seeds that lost the bidding
             continue;
         }
+        // collect seed hits and reject those that lost the hit bidding
+        char best_for_hit = 0;
         Tracklet seed;
         seed.size = 0;
         // dummy path to start the loop
-        int2 path = make_int2(0, d_seed_proposals[prop_idx].y);
+        int2 prop = d_seed_proposals[prop_idx];
+        int2 path = make_int2(0, prop.y);
         while (path.y >= 0) {
             path = d_path_store[path.y];
             seed.nodes[seed.size++] =
                 d_output_graph[traccc::device::gbts_consts::node1 +
                                edge_size * path.x];
+            best_for_hit +=
+                (prop_idx ==
+                 (d_hit_bids[seed.nodes[seed.size - 1]] & 0xFFFFFFFFLL));
         }
         seed.nodes[seed.size++] =
             d_output_graph[traccc::device::gbts_consts::node2 +
                            edge_size * path.x];
-        // sample begining, middle, end sp from tracklet for now
-        seeds_device.push_back({seed.nodes[seed.size - 1],
-                                seed.nodes[(1 + seed.size) / 2 - 1],
-                                seed.nodes[0]});
+        best_for_hit += (prop_idx == (d_hit_bids[seed.nodes[seed.size - 1]] &
+                                      0xFFFFFFFFLL));
+
+        if ((best_for_hit < best_hit_frac * seed.size)) {
+            continue;
+        }
+        char diff_code = 0;
+        bool force_dropout = false;
+        if (use_dropout) {
+            float4 sps[3];
+            // seed 1
+            sps[0] = d_sp_params[seed.nodes[seed.size - 1]];
+            sps[1] = d_sp_params[seed.nodes[(seed.size - 1) / 2 + 1]];
+            sps[2] = d_sp_params[seed.nodes[0]];
+            float2 curv_cot_1 = estimate_params(sps);
+            // seed 2
+            sps[1] = d_sp_params[seed.nodes[(seed.size - 1) / 2]];
+            float2 curv_cot_2 = estimate_params(sps);
+            sps[0] = d_sp_params[seed.nodes[seed.size - 2]];
+            // seed 3
+            float2 curv_cot_3 = estimate_params(sps);
+            // for low eta (higher fake rate) seeds perform a stronger cut
+            if ((best_for_hit < seed.size - 1) &
+                (abs(curv_cot_1.y + curv_cot_2.y + curv_cot_3.y) <
+                 3.0f * tight_bid_cot_threshold) &
+                (seed.size < 5)) {
+                continue;
+            }
+            float diff[3] = {abs(curv_cot_1.x - curv_cot_2.x),
+                             abs(curv_cot_2.x - curv_cot_3.x),
+                             abs(curv_cot_1.x - curv_cot_3.x)};
+            diff_code = 4 * (diff[0] < dcurv_cut_m) +
+                        2 * (diff[1] < dcurv_cut_m) + (diff[2] < dcurv_cut_m);
+            // for high pt the diff may pass dispite bad estimates
+            force_dropout = abs(curv_cot_1.x + curv_cot_2.x + curv_cot_3.x) <
+                            3.0f * force_dropout_max_curv_m;
+            force_dropout |= (abs(curv_cot_1.y + curv_cot_2.y + curv_cot_3.y) <
+                              3.0f * tight_bid_cot_threshold) &
+                             diff_code == 0;
+        }
+        // use one seed from a consistant pair/set + the inconsistant one
+        // sample spacepoints from tracklet to create seeds
+        // include 1st order unless either 2 or 3 are consitant with the other
+        // and 1
+        if (diff_code != 3 & diff_code != 6 | force_dropout) {
+            seeds_device.push_back({seed.nodes[seed.size - 1],
+                                    seed.nodes[(seed.size - 1) / 2 + 1],
+                                    seed.nodes[0]});
+        }
+        // include 2nd order if it consistant with 1 and 3 or only 1 and 3 are
+        // consistant
+        if (diff_code == 1 | diff_code == 6) {
+            seeds_device.push_back({seed.nodes[seed.size - 1],
+                                    seed.nodes[(seed.size - 1) / 2],
+                                    seed.nodes[0]});
+        }
+        // include 3rd order if it is consistant with 1 and 2 or only 1 and 2
+        // are consistant or if only 2 and 3 are consistant
+        if (diff_code == 2 | diff_code == 3 | diff_code == 4 | force_dropout) {
+            seeds_device.push_back({seed.nodes[seed.size - 2],
+                                    seed.nodes[(seed.size - 1) / 2],
+                                    seed.nodes[0]});
+        }
     }
 }