standalone/test/binnor_demo.cxx

0001 #include <random>
0002 #include <chrono>
0003
0004 #include "../../interface/radix_sort.h"
0005 #include "../../src/radix_sort.cc"
0006 #include "../../interface/binnor.h"
0007
0008 // build as:
0009 // c++ -o binnor_demo -O3 -mavx2 -std=c++17 -I../../../.. binnor_demo.cxx
0010
0011 bool print_axis_binnor = true;
0012 bool print_bin_contents = true;
0013
0014 void run_scan();
0015 void run_test(const int NLoop, int NN, const bool use_radix = true, const bool keep_cons = false);
0016
0017 using namespace mkfit;
0018
0019 int main()
0020 {
0021   bool use_radix = true;
0022   bool keep_cons = false;
0023   run_test(100, 10000, use_radix, keep_cons);
0024
0025   // Run scan over low N-elements to see where radix starts to make more sense.
0026   // NOTE: In binnor.h, finalize_registration(), comment-out the existing heuristics
0027   // that uses std::sort over radix when NN < 256 (128 with keep_cons).
0028   // run_scan();
0029
0030 }
0031
0032 // -------------------------------------------------------------
0033
0034 void run_scan()
0035 {
0036   print_axis_binnor = false;
0037   print_bin_contents = false;
0038
0039   const bool radix[] = { true, true, false, false };
0040   const bool keepc[] = { false, true, false, true };
0041
0042   for (int t = 0; t < 4; ++t) {
0043     printf("RADIX=%s  KEEP-CONS=%s\n", radix[t] ? "true" : "false", keepc[t] ? "true" : "false");
0044     for (int i = 100; i <= 400; i += 10) {
0045       run_test(100000, i, radix[t], keepc[t]);
0046     }
0047     printf("\n");
0048   }
0049 }
0050
0051 // -------------------------------------------------------------
0052
0053 void run_test(const int NLoop, int NN, const bool use_radix, const bool keep_cons) {
0054   constexpr float PI = 3.14159265358979323846;
0055   constexpr float TwoPI = 6.28318530717958647692;
0056   constexpr float PIOver2 = PI / 2.0f;
0057   constexpr float PIOver4 = PI / 4.0f;
0058
0059   axis_pow2_u1<float, unsigned short, 16, 8> phi(-PI, PI);
0060   axis<float, unsigned short, 16, 8> eta(-2.6, 2.6, 20u);
0061   binnor<unsigned int, decltype(phi), decltype(eta), 18, 14> b(phi, eta, use_radix, keep_cons);
0062
0063   if (print_axis_binnor) {
0064     printf("Axis phi: M-bits=%d, N-bits=%d  Masks M:0x%x N:0x%x\n", phi.c_M, phi.c_N, phi.c_M_mask, phi.c_N_mask);
0065     printf("Axis eta: M-bits=%d, N-bits=%d  m_bins=%d n_bins=%d\n", eta.c_M, eta.c_N, eta.size_of_M(), eta.size_of_N());
0066     printf("Have binnor, size of vec = %zu, sizeof(C_pair) = %d\n", b.m_bins.size(), sizeof(decltype(b)::C_pair));
0067   }
0068   /*
0069     for (float p = -TwoPI; p < TwoPI; p += TwoPI / 15.4f) {
0070         printf("  phi=%-9f m=%5d n=%3d m2n=%3d n_safe=%3d\n", p,
0071                phi.from_R_to_M_bin(p), phi.from_R_to_N_bin(p),
0072                phi.from_M_bin_to_N_bin( phi.from_R_to_M_bin(p) ),
0073                phi.from_R_to_N_bin_safe(p) );
0074     }
0075     */
0076
0077   std::mt19937 rnd(std::random_device{}());
0078   std::uniform_real_distribution<float> d_phi(-PI, PI);
0079   std::uniform_real_distribution<float> d_eta(-2.55, 2.55);
0080
0081   struct track {
0082     float phi, eta;
0083   };
0084   std::vector<track> tracks;
0085   tracks.reserve(NN);
0086   for (int i = 0; i < NN; ++i) {
0087     tracks.push_back({d_phi(rnd), d_eta(rnd)});
0088     // printf("made track %3d:  phi=%f  eta=%f\n", i, tracks.back().phi, tracks.back().eta);
0089   }
0090
0091   std::chrono::duration<long long, std::nano> duration, fill_duration, sort_duration;
0092   duration = fill_duration = sort_duration = std::chrono::nanoseconds::zero();
0093
0094   for (int n_loop = 0; n_loop < NLoop; ++n_loop) {
0095     b.reset_contents(false); // do not shrink the vectors
0096
0097     auto start = std::chrono::high_resolution_clock::now();
0098
0099     b.begin_registration(NN);  // NN optional, reserves construction vector
0100
0101     for (int i = 0; i < NN; ++i) {
0102       b.register_entry(tracks[i].phi, tracks[i].eta);
0103     }
0104
0105     auto reg_start = std::chrono::high_resolution_clock::now();
0106
0107     b.finalize_registration();
0108
0109     auto stop = std::chrono::high_resolution_clock::now();
0110
0111     duration += std::chrono::duration_cast<std::chrono::nanoseconds>(stop - start);
0112     fill_duration += std::chrono::duration_cast<std::chrono::nanoseconds>(reg_start - start);
0113     sort_duration += std::chrono::duration_cast<std::chrono::nanoseconds>(stop - reg_start);
0114   }
0115
0116   // for (int i = 0; i < NN; ++i)
0117   // {
0118   //     const track &t = tracks[ b.m_ranks[i] ];
0119   //     printf("%3d  %3d  phi=%f  eta=%f\n", i, b.m_ranks[i], t.phi, t.eta);
0120   // }
0121
0122   if (print_bin_contents) {
0123     printf("\n\n--- Single bin access for (phi, eta) = (0,0):\n\n");
0124     auto nbin = b.get_n_bin(0.f, 0.f);
0125     auto cbin = b.get_content(0.f, 0.f);
0126     printf("For (phi 0, eta 0; %u, %u) got first %d, count %d\n", nbin.bin1(), nbin.bin2(), cbin.first, cbin.count);
0127     for (auto i = cbin.first; i < cbin.first + cbin.count; ++i) {
0128       const track &t = tracks[b.m_ranks[i]];
0129       printf("%3d  %3d  phi=%f  eta=%f\n", i, b.m_ranks[i], t.phi, t.eta);
0130     }
0131
0132     printf("\n\n--- Range access for phi=[(-PI+0.02 +- 0.1], eta=[1.3 +- .2]:\n\n");
0133     auto phi_rng = phi.from_R_rdr_to_N_bins(-PI + 0.02, 0.1);
0134     auto eta_rng = eta.from_R_rdr_to_N_bins(1.3, .2);
0135     printf("phi bin range: %u, %u; eta %u, %u\n", phi_rng.begin, phi_rng.end, eta_rng.begin, eta_rng.end);
0136     for (auto i_phi = phi_rng.begin; i_phi != phi_rng.end; i_phi = phi.next_N_bin(i_phi)) {
0137       for (auto i_eta = eta_rng.begin; i_eta != eta_rng.end; i_eta = eta.next_N_bin(i_eta)) {
0138         printf(" at i_phi=%u, i_eta=%u\n", i_phi, i_eta);
0139         const auto cbin = b.get_content(i_phi, i_eta);
0140         for (auto i = cbin.begin(); i < cbin.end(); ++i) {
0141           const track &t = tracks[b.m_ranks[i]];
0142           printf("   %3d  %3d  phi=%f  eta=%f\n", i, b.m_ranks[i], t.phi, t.eta);
0143         }
0144       }
0145     }
0146     printf("\n");
0147   }
0148
0149   printf("Processing time for %d times %d points: %f sec (filling %f sec, sort + binning %f sec)\n",
0150          NLoop, NN,
0151          1.e-9 * duration.count(),
0152          1.e-9 * fill_duration.count(),
0153          1.e-9 * sort_duration.count());
0154 }