.. _program_listing_file_larflow_PrepFlowMatchData_PrepMatchTriplets.cxx:

Program Listing for File PrepMatchTriplets.cxx
==============================================

|exhale_lsh| :ref:`Return to documentation for file <file_larflow_PrepFlowMatchData_PrepMatchTriplets.cxx>` (``larflow/PrepFlowMatchData/PrepMatchTriplets.cxx``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #include "PrepMatchTriplets.h"
   #include "FlowTriples.h"
   
   #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
   #include <numpy/ndarrayobject.h>
   
   #include <sstream>
   #include <ctime>
   #include <algorithm>    // std::shuffle
   #include <random>       // std::default_random_engine
   #include <chrono>       // std::chrono::system_clock
   
   #include "larlite/core/LArUtil/LArProperties.h"
   #include "larlite/core/LArUtil/Geometry.h"
   #include "larcv/core/PyUtil/PyUtils.h"
   #include "larflow/LArFlowConstants/LArFlowConstants.h"
   #include "ublarcvapp/UBWireTool/UBWireTool.h"
   #include "ublarcvapp/UBImageMod/EmptyChannelAlgo.h"
   
   #include "TRandom3.h"
   
   namespace larflow {
   namespace prep {
   
     bool PrepMatchTriplets::_setup_numpy = false;
     
     void PrepMatchTriplets::clear()
     {
   
       _imgmeta_v.clear();
       _sparseimg_vv.clear();
       _triplet_v.clear();
       _truth_v.clear();
       _truth_2plane_v.clear();
       _weight_v.clear();
       _flowdir_v.clear();
       _triarea_v.clear();
       _pos_v.clear();
       
     }
     
     void PrepMatchTriplets::process( larcv::IOManager& iolcv,
                                      std::string wire_producer,
                                      std::string chstatus_producer,
                                      const float adc_threshold,
                                      const bool check_wire_intersection )
     {
       
       // create bad channel maker algo
       ublarcvapp::EmptyChannelAlgo badchmaker;
   
       // get images
       larcv::EventImage2D* ev_adc
         = (larcv::EventImage2D*)iolcv.get_data(larcv::kProductImage2D,wire_producer);
       
       // get chstatus
       larcv::EventChStatus* ev_badch
         = (larcv::EventChStatus*)iolcv.get_data(larcv::kProductChStatus,chstatus_producer);
       
       std::vector<larcv::Image2D> badch_v = badchmaker.makeBadChImage( 4, 3, 2400, 6*1008, 3456, 6, 1, *ev_badch );
       
       process( ev_adc->Image2DArray(), badch_v, adc_threshold, check_wire_intersection );
       
     } 
     
     void PrepMatchTriplets::process( const std::vector<larcv::Image2D>& adc_v,
                                      const std::vector<larcv::Image2D>& badch_v,
                                      const float adc_threshold,
                                      const bool check_wire_intersection )
     {
   
       std::clock_t start = std::clock();
       
       // first we make a common sparse image
       _sparseimg_vv = larflow::prep::FlowTriples::make_initial_sparse_image( adc_v, adc_threshold );
       _imgmeta_v.clear();
       for ( auto const& img : adc_v )
         _imgmeta_v.push_back( img.meta() );
         
   
       // then we make the flow triples for each of the six flows
       // we order the flows based on the quality of the planes, with the Y plane being better
       
       FlowDir_t flow_order[] = { kY2V, kY2U, kV2Y, kU2Y, kU2V, kV2U };
       
       std::vector< larflow::prep::FlowTriples > triplet_v( larflow::kNumFlows );
       for (int flowindex=0; flowindex<(int)larflow::kNumFlows; flowindex++) {
   
         // if ( flowindex!=kV2Y )
         //   continue;
         
         int sourceplane, targetplane;
         larflow::LArFlowConstants::getFlowPlanes( (larflow::FlowDir_t)flowindex, sourceplane, targetplane );
   
         // save pixel position, rather than index: we will reindex after adding dead channels
         triplet_v[flowindex]  = FlowTriples( sourceplane, targetplane,
                                              adc_v, badch_v,
                                              _sparseimg_vv, 10.0, false );
       }
   
       // collect the unique dead channel additions to each plane
       std::set< std::pair<int,int> > deadpixels_to_add[ adc_v.size() ];
   
       // add the unique dead channel additions
       std::cout << "sparse pixel totals before deadch additions: "
                 << "(" << _sparseimg_vv[0].size() << "," << _sparseimg_vv[1].size() << "," << _sparseimg_vv[2].size() << ")"
                 << std::endl;
   
       for (int flowindex=0; flowindex<(int)larflow::kNumFlows; flowindex++) {
         
         // if ( flowindex!=kV2Y )
         //   continue;
   
         auto& triplet = triplet_v[flowindex];
         int otherplane = triplet.get_other_plane_index();
         std::vector<FlowTriples::PixData_t>& pix_v = triplet.getDeadChToAdd()[ otherplane ];
         
         for ( auto& pix : pix_v ) {
           auto it = deadpixels_to_add[ otherplane ].find( std::pair<int,int>( pix.row, pix.col ) );
           if ( it==deadpixels_to_add[ otherplane ].end() ) {
             // unique, add to sparse image
             pix.val = 1.0;
             _sparseimg_vv[otherplane].push_back( pix );
           }
         }
       }
       
       std::cout << "sparse pixel totals after deadch additions: "
                 << "(" << _sparseimg_vv[0].size() << "," << _sparseimg_vv[1].size() << "," << _sparseimg_vv[2].size() << ")"
                 << std::endl;
       
       // sort all pixels
       for ( auto& pix_v : _sparseimg_vv ) {
         std::sort( pix_v.begin(), pix_v.end() );
       }
       
       // condense and reindex matches
       std::set< std::vector<int> > triplet_set;
       _triplet_v.clear();
       _flowdir_v.clear();
       _triarea_v.clear();
       _pos_v.clear();
       _triplet_v.reserve( 200000 );
       _flowdir_v.reserve( 200000 );
       _triarea_v.reserve( 200000 );
       _pos_v.reserve( 200000 );
   
       int   n_not_crosses = 0;
       float n_bad_triarea = 0;
   
       for (int iflow=0; iflow<(int)larflow::kNumFlows; iflow++) {
         
         larflow::FlowDir_t flowdir=flow_order[iflow];
   
         // for debug
         // if ( flowdir!=kY2V )
         //   continue;
         
         auto& triplet_data = triplet_v[ (int)flowdir ];
         
         int srcplane = triplet_data.get_source_plane_index();
         int tarplane = triplet_data.get_target_plane_index();
         int othplane = triplet_data.get_other_plane_index();
   
         std::cout << "[PrepMatchTriplets] combine matches from flow "
                   << larflow::LArFlowConstants::getFlowName(flowdir)
                   << " planes={" << srcplane << " -> " << tarplane << ", " << othplane << "}"
                   << std::endl;            
         
         for ( auto& trip : triplet_data.getTriples() ) {
   
           std::vector<FlowTriples::PixData_t> pix_v(3);
           pix_v[ srcplane ] = FlowTriples::PixData_t( trip[3], trip[0], 0.0 );
           pix_v[ tarplane ] = FlowTriples::PixData_t( trip[3], trip[1], 0.0 );
           pix_v[ othplane ] = FlowTriples::PixData_t( trip[3], trip[2], 0.0 );
   
           auto it_src = std::lower_bound( _sparseimg_vv[srcplane].begin(), _sparseimg_vv[srcplane].end(), pix_v[ srcplane ] );
           auto it_tar = std::lower_bound( _sparseimg_vv[tarplane].begin(), _sparseimg_vv[tarplane].end(), pix_v[ tarplane ] );
           auto it_oth = std::lower_bound( _sparseimg_vv[othplane].begin(), _sparseimg_vv[othplane].end(), pix_v[ othplane ] );
   
           if ( it_src==_sparseimg_vv[srcplane].end()
                || it_tar==_sparseimg_vv[tarplane].end()
                || it_oth==_sparseimg_vv[othplane].end() ) {
             std::stringstream ss;
             ss << "Did not find one of sparse image pixels for col triplet=(" << trip[0] << "," << trip[1] << "," << trip[2] << ")";
             ss << " found-index=("
                << it_src-_sparseimg_vv[srcplane].begin() << ","
                << it_tar - _sparseimg_vv[tarplane].begin() << ","
                << it_oth - _sparseimg_vv[othplane].begin() << ")"
                << std::endl;
             throw std::runtime_error( ss.str() );
           }
           if ( (*it_src).row!=trip[3]
                || (*it_tar).row!=trip[3]
                || (*it_oth).row!=trip[3]
                || (*it_src).col!=trip[0]
                || (*it_tar).col!=trip[1]
                || (*it_oth).col!=trip[2] ) {             
             
             std::stringstream ss;
             ss << "found the wrong pixel. searching for triplet=(" << trip[0] << "," << trip[1] << "," << trip[2] << "," << trip[3] << ") and got"
                << " src(c,r)=(" << (*it_src).col << "," << (*it_src).row << ")"
                << " tar(c,r)=(" << (*it_tar).col << "," << (*it_tar).row << ")"
                << " oth(c,r)=(" << (*it_oth).col << "," << (*it_oth).row << ")"
                << std::endl;
             throw std::runtime_error( ss.str() );
           }
                
   
           std::vector<int> imgcoord_v(4);
           imgcoord_v[ srcplane ] = trip[0];
           imgcoord_v[ tarplane ] = trip[1];
           imgcoord_v[ othplane ] = trip[2];
           imgcoord_v[ 3 ]        = trip[3];
   
           // check triplet 3d consistency, using ublarcvapp tool
           int crosses = 1;
           double tri_area = 0.;
           std::vector<float> intersection = {0,0,0};
           
           if ( check_wire_intersection ) {
   
             UInt_t src_ch = larutil::Geometry::GetME()->PlaneWireToChannel( (UInt_t)srcplane, (UInt_t)trip[0]);
             UInt_t tar_ch = larutil::Geometry::GetME()->PlaneWireToChannel( (UInt_t)tarplane, (UInt_t)trip[1] );
             Double_t y,z;
             bool crosses  = larutil::Geometry::GetME()->ChannelsIntersect( src_ch, tar_ch, y, z );
   
             if ( !crosses ) {
               n_not_crosses++;
               continue;
             }          
             intersection[0] = (adc_v[0].meta().pos_y( imgcoord_v[3] )-3200.0)*0.5*larutil::LArProperties::GetME()->DriftVelocity();
             intersection[1] = y;
             intersection[2] = z;
             
           }
           
           
           auto it_trip = triplet_set.find( imgcoord_v );
           if ( it_trip==triplet_set.end() ) {
             triplet_set.insert( imgcoord_v );
   
             std::vector<int> imgindex_v(4);
             imgindex_v[ srcplane ] = it_src - _sparseimg_vv[srcplane].begin();
             imgindex_v[ tarplane ] = it_tar - _sparseimg_vv[tarplane].begin();
             imgindex_v[ othplane ] = it_oth - _sparseimg_vv[othplane].begin();
             imgindex_v[ 3 ]        = trip[3];
   
             _triplet_v.push_back( imgindex_v );
             _flowdir_v.push_back( flowdir );
             if ( check_wire_intersection) {
               _triarea_v.push_back( tri_area );
               _pos_v.push_back( intersection );
             }
           }
           
         }
   
       }
   
       std::clock_t end = std::clock();
       std::cout << "[PrepMatchTriplets] made total of " << _triplet_v.size()
                 << " unique index triplets. time elapsed=" << float(end-start)/float(CLOCKS_PER_SEC)
                 << std::endl;
       std::cout << "[PrepMatchTriplets] number removed for not intersecting: " << n_not_crosses << std::endl;
   
     }//end of process method
   
     std::vector<TH2D> PrepMatchTriplets::plot_sparse_images( std::string hist_stem_name )
     {
       std::vector<TH2D> out_v;
       for ( int p=0; p<(int)_imgmeta_v.size(); p++ ) {
         std::stringstream ss;
         ss << "htriples_plane" << p << "_" << hist_stem_name;
         auto const& meta = _imgmeta_v[p];
         TH2D hist( ss.str().c_str(), "",
                    meta.cols(), meta.min_x(), meta.max_x(),
                    meta.rows(), meta.min_y(), meta.max_y() );
   
         for ( auto const& pix : _sparseimg_vv[p] ) {
           hist.SetBinContent( pix.col+1, pix.row+1, pix.val );
         }
         
         out_v.emplace_back(std::move(hist));
       }
       return out_v;
     }
   
     void PrepMatchTriplets::make_truth_vector( const std::vector<larcv::Image2D>& larflow_v )
     {
   
       _truth_v.clear();
       _truth_2plane_v.clear();
       _truth_v.resize( _triplet_v.size(), 0 );
       _truth_2plane_v.resize( _triplet_v.size() );
   
       const int true_match_span = 3;
       const int min_required_connections = 1;
   
       int ntriplet_truth = 0;
       std::vector< int > ndoublet_truth( (int)larflow::kNumFlows, 0 );
       
       for ( size_t itrip=0; itrip<_triplet_v.size(); itrip++ ) {
         // for each triplet, we look for truth flows that connect the planes
         auto const& triplet = _triplet_v[itrip];
         larflow::FlowDir_t flow_dir_origin = _flowdir_v[itrip];
         _truth_2plane_v[itrip].resize( (int)larflow::kNumFlows, 0 );
         
         // for debug
         //if ( flow_dir_origin!=kY2V ) continue;
   
         int srcplane, tarplane;
         larflow::LArFlowConstants::getFlowPlanes( flow_dir_origin, srcplane, tarplane );
         int othplane = larflow::LArFlowConstants::getOtherPlane( srcplane, tarplane );
         
         std::vector< const FlowTriples::PixData_t* > pix_v( _sparseimg_vv.size() );
         pix_v[srcplane] = &_sparseimg_vv[srcplane][ triplet[srcplane] ];
         pix_v[tarplane] = &_sparseimg_vv[tarplane][ triplet[tarplane] ];
         pix_v[othplane] = &_sparseimg_vv[othplane][ triplet[othplane] ]; 
   
         int ngood_connections = 0;
         float pixflow  = larflow_v[(int)flow_dir_origin].pixel( pix_v[srcplane]->row, pix_v[srcplane]->col );
         int target_col = pix_v[srcplane]->col + (int)pixflow;
         if ( fabs(pixflow)!=0 && pixflow>-3999 && abs(target_col-pix_v[tarplane]->col)<true_match_span ) {
           ngood_connections++;
           _truth_2plane_v[itrip][(int)flow_dir_origin] = 1;
           ndoublet_truth[(int)flow_dir_origin]++;
         }
         
         if ( ngood_connections>=min_required_connections ) {
           _truth_v[itrip] = 1;
           ntriplet_truth++;
         }
         else {
           _truth_v[itrip] = 0;
         }
               
       }//end of trips loop
       
       std::cout << "[PrepMatchTriplets::make_truth_vector] " << std::endl;
       std::cout << "  number of triplets: " << _triplet_v.size() << std::endl;
       std::cout << "  number of sparse pixels: [ "
                 << _sparseimg_vv[0].size() << ", "
                 << _sparseimg_vv[1].size() << ", "
                 << _sparseimg_vv[2].size() << " ]"
                 << std::endl;
       std::cout << "  number of true-match triplets: " << ntriplet_truth << std::endl;
       std::cout << "  doublet truth: [";
       for (auto& n : ndoublet_truth ) std::cout << " " << n << ",";
       std::cout << " ]" << std::endl;
       
     }
   
     std::vector<TH2D> PrepMatchTriplets::plot_truth_images( std::string hist_stem_name )
     {
       std::vector<TH2D> out_v;
   
       for ( int p=0; p<(int)_imgmeta_v.size(); p++ ) {
         std::stringstream ss;
         ss << "htriples_truth_plane" << p << "_" << hist_stem_name;
         auto const& meta = _imgmeta_v[p];
         TH2D hist( ss.str().c_str(), "",
                    meta.cols(), meta.min_x(), meta.max_x(),
                    meta.rows(), meta.min_y(), meta.max_y() );
   
         out_v.emplace_back( std::move(hist) );
       }
       
       for (int i=0; i<(int)_triplet_v.size(); i++ ) {
         auto& trip  = _triplet_v[i];
         auto& truth = _truth_v[i];
         std::vector< const FlowTriples::PixData_t* > pix_v( trip.size(), 0 );
         for (int p=0; p<(int)_imgmeta_v.size(); p++ ) {
           pix_v[p] = &_sparseimg_vv[p][ trip[p] ];
           int col = pix_v[p]->col+1;
           int row = pix_v[p]->row+1;
           if ( out_v[p].GetBinContent( col+1, row+1 )<10 )
             out_v[p].SetBinContent( col+1, row+1, 1 + 10*truth );
         }
       }
       return out_v;
     }
   
     PyObject* PrepMatchTriplets::make_sparse_image( int plane ) {
       
       if ( !_setup_numpy ) {
         import_array1(0);
         _setup_numpy = true;
       }
   
       npy_intp* dims = new npy_intp[2];
       dims[0] = (int)_sparseimg_vv[plane].size();
   
       // if we want truth, we include additional value with 1=correct match, 0=false    
       dims[1] = 3;
   
       // output array
       PyArrayObject* array = (PyArrayObject*)PyArray_SimpleNew( 2, dims, NPY_FLOAT );
   
       for ( size_t idx=0; idx<_sparseimg_vv[plane].size(); idx++ ) {
         *((float*)PyArray_GETPTR2( array, (int)idx, 0)) = (float)_sparseimg_vv[plane][idx].row;
         *((float*)PyArray_GETPTR2( array, (int)idx, 1)) = (float)_sparseimg_vv[plane][idx].col;
         *((float*)PyArray_GETPTR2( array, (int)idx, 2)) = (float)_sparseimg_vv[plane][idx].val;      
       }
   
       
       return (PyObject*)array;
     }
   
     PyObject* PrepMatchTriplets::make_2plane_match_array( larflow::FlowDir_t kdir,
                                                           const int max_num_samples,
                                                           const std::vector<int>& idx_v,
                                                           const int start_idx,
                                                           const bool withtruth,
                                                           int& nsamples )
     {
   
       if ( !_setup_numpy ) {
         import_array1(0);
         _setup_numpy = true;
       }
   
       int nd = 2;
       int ndims2 = 3;
       npy_intp dims[] = { max_num_samples, ndims2 };
   
       // output array
       PyArrayObject* array = (PyArrayObject*)PyArray_SimpleNew( nd, dims, NPY_LONG );
   
       int srcplane,tarplane;
       larflow::LArFlowConstants::getFlowPlanes( kdir, srcplane, tarplane );
   
       // number of pairs we've stored
       nsamples = 0;
       
       int end_idx = start_idx + max_num_samples;
       end_idx = ( end_idx>(int)idx_v.size() )   ?  idx_v.size() : end_idx; // cap to number of indices
   
       std::cout << "[PrepMatchTriplets::make_2plane_match_array] src=" << srcplane << " tar=" << tarplane << " withtruth=" << withtruth << " "
                 << "make numpy array with indices from triplets[" << start_idx << ":" << end_idx << "]" << std::endl;
       
       for ( int idx=start_idx; idx<end_idx; idx++ ) {
         int tripidx = idx_v[idx];
         *((long*)PyArray_GETPTR2( array, nsamples, 0)) = (long)_triplet_v[tripidx][srcplane];
         *((long*)PyArray_GETPTR2( array, nsamples, 1)) = (long)_triplet_v[tripidx][tarplane];
         if ( withtruth ) {
           *((long*)PyArray_GETPTR2( array, nsamples, 2)) = (long)_truth_2plane_v[tripidx][(int)kdir];
         }
         else {
           *((long*)PyArray_GETPTR2( array, nsamples, 2)) = 0;
         }
         nsamples++;
         if (nsamples==max_num_samples)
           break;
   
       }//end of indices loop
   
       std::cout << "[PrepMatchTriplets::make_2plane_match_array] nsamples=" << nsamples << std::endl;
   
       // zero rest of array
       if ( nsamples<max_num_samples ) {
         for ( size_t i=nsamples; i<max_num_samples; i++ ) {
           for (int j=0; j<dims[1]; j++) {
             *((long*)PyArray_GETPTR2( array, i, j)) = 0;
           }
         }
       }
       
       // return the array
       return (PyObject*)array;
       
     }
   
     PyObject* PrepMatchTriplets::sample_2plane_matches( larflow::FlowDir_t kdir,
                                                         const int& nsamples,
                                                         int& nfilled,
                                                         bool withtruth ) {
   
       std::vector<int> idx_v( _triplet_v.size() );
       for ( size_t i=0; i<_triplet_v.size(); i++ ) idx_v[i] = (int)i;
       unsigned seed =  std::chrono::system_clock::now().time_since_epoch().count();
       shuffle (idx_v.begin(), idx_v.end(), std::default_random_engine(seed));
   
       return make_2plane_match_array( kdir, nsamples, idx_v, 0, withtruth, nfilled );
   
     }
   
     PyObject* PrepMatchTriplets::get_chunk_2plane_matches( larflow::FlowDir_t kdir,
                                                            const int& start_index,
                                                            const int& max_num_pairs,
                                                            int& last_index,
                                                            int& num_pairs_filled,
                                                            bool with_truth ) {
       
       std::vector<int> idx_v( max_num_pairs, 0 );
       last_index = start_index + max_num_pairs;
       last_index = ( last_index>(int)_triplet_v.size() ) ? (int)_triplet_v.size() : last_index;
       
       for ( int i=start_index; i<last_index; i++ ) {
         idx_v[i-start_index] = (int)i;
       }
   
       return make_2plane_match_array( kdir, max_num_pairs, idx_v, 0, with_truth, num_pairs_filled );
   
     }
   
     PyObject* PrepMatchTriplets::make_triplet_array( const int max_num_samples,
                                                      const std::vector<int>& idx_v,
                                                      const int start_idx,
                                                      const bool withtruth,
                                                      int& nsamples )
     {
   
       if ( !_setup_numpy ) {
         import_array1(0);
         _setup_numpy = true;
       }
   
       int nd = 2;
       int ndims2 = 5;
       npy_intp dims[] = { max_num_samples, ndims2 };
   
       // output array
       PyArrayObject* array = (PyArrayObject*)PyArray_SimpleNew( nd, dims, NPY_LONG );
   
       // number of pairs we've stored
       nsamples = 0;
       
       int end_idx = start_idx + max_num_samples;
       end_idx = ( end_idx>(int)idx_v.size() )   ?  idx_v.size() : end_idx; // cap to number of indices
   
       std::cout << "[PrepMatchTriplets::make_triplet_array] withtruth=" << withtruth << " "
                 << "make numpy array with indices from triplets[" << start_idx << ":" << end_idx << "]"
                 << std::endl;
       
       for ( int idx=start_idx; idx<end_idx; idx++ ) {
         int tripidx = idx_v[idx];
         for (size_t p=0; p<3; p++ )
           *((long*)PyArray_GETPTR2( array, nsamples, p)) = (long)_triplet_v[tripidx][p];
         if ( withtruth ) {
           *((long*)PyArray_GETPTR2( array, nsamples, 3)) = (long)_truth_v[tripidx];
           *((long*)PyArray_GETPTR2( array, nsamples, 4)) = (long)tripidx;
         }
         else {
           *((long*)PyArray_GETPTR2( array, nsamples, 3)) = 0;        
           *((long*)PyArray_GETPTR2( array, nsamples, 4)) = (long)tripidx;
         }
         nsamples++;
         if (nsamples==max_num_samples)
           break;
       }//end of indices loop
   
       std::cout << "[PrepMatchTriplets::make_triplet_array] nsamples=" << nsamples << std::endl;
   
       // zero rest of array
       if ( nsamples<max_num_samples ) {
         for ( size_t i=nsamples; i<max_num_samples; i++ ) {
           for (int j=0; j<dims[1]; j++) {
             *((long*)PyArray_GETPTR2( array, i, j)) = 0;
           }
         }
       }
   
       
       
       // return the array
       return (PyObject*)array;
       
     }
   
   
     PyObject* PrepMatchTriplets::sample_triplet_matches( const int& nsamples,
                                                          int& nfilled,
                                                          bool withtruth ) {
   
       std::vector<int> idx_v( _triplet_v.size() );
       for ( size_t i=0; i<_triplet_v.size(); i++ ) idx_v[i] = (int)i;
       unsigned seed =  std::chrono::system_clock::now().time_since_epoch().count();
       shuffle (idx_v.begin(), idx_v.end(), std::default_random_engine(seed));
       
       return make_triplet_array( nsamples, idx_v, 0, withtruth, nfilled );
   
     }
     
     PyObject* PrepMatchTriplets::get_chunk_triplet_matches( const int& start_index,
                                                             const int& max_num_pairs,
                                                             int& last_index,
                                                             int& num_pairs_filled,
                                                             bool with_truth ) {
       
       last_index = start_index + max_num_pairs;
       last_index = ( last_index>(int)_triplet_v.size() ) ? (int)_triplet_v.size() : last_index;
       num_pairs_filled = last_index-start_index;
       std::vector<int> idx_v( num_pairs_filled, 0 );    
       
       for ( int i=start_index; i<last_index; i++ ) {
         idx_v[i-start_index] = (int)i;
       }
   
       return make_triplet_array( max_num_pairs, idx_v, 0, with_truth, num_pairs_filled );
   
     }
   
     PyObject* PrepMatchTriplets::sample_hard_example_matches( const int& nsamples,
                                                               const int& nhard_samples,
                                                               PyObject* triplet_scores,                                                            
                                                               int& nfilled,
                                                               bool withtruth ) {
   
       // if number of triplets less than the requested triplet sample, just pass the triplets back
       if (nsamples<=(int)_triplet_v.size()) {
         return sample_triplet_matches( nsamples, nfilled, withtruth );
       }
       
       std::vector<int> idx_v( _triplet_v.size(), 0 );
       for ( size_t i=0; i<_triplet_v.size(); i++ ) idx_v[i] = (int)i;
       unsigned seed =  std::chrono::system_clock::now().time_since_epoch().count();
       shuffle (idx_v.begin(), idx_v.end(), std::default_random_engine(seed));
   
   
       TRandom3 sampler( seed );
       // now we sample for bad events
       int nbad = 0;
       int nsaved =  0;
       std::set<int> saved_idx_set;
       std::vector<int> saved_idx_v( _triplet_v.size(), 0 );
       for (size_t i=0; i<_triplet_v.size(); i++) {
         int idx = idx_v[i];
         float past_score = *((float*)PyArray_GETPTR1( (PyArrayObject*)triplet_scores, (int)idx ));
         float weight = fabs(past_score-(float)_truth_v[idx]);
         if (nbad<nhard_samples) {
           // first try to get a set number of bad examples
           if ( sampler.Uniform()>weight ) {
             saved_idx_v[nsaved] = idx;        
             saved_idx_set.insert(idx);
             nsaved++;
             if ( weight>0.5 )
               nbad++;
           }
         }
         else {
           saved_idx_v[nsaved] = idx;
           saved_idx_set.insert(idx);
           nsaved++;
         }
       }
   
       if ( nsaved<nsamples ) {
         // go back and grab triplets we haven't used
         for (size_t i=0; i<_triplet_v.size(); i++) {
           int idx = idx_v[i];        
           if ( saved_idx_set.find(idx)!=saved_idx_set.end() ) {
             continue;
           }
           saved_idx_v[nsaved] = idx;
           saved_idx_set.insert(idx);
           nsaved++;
           if ( nsaved==nsamples )
             break;
         }
       }
       saved_idx_v.resize(nsaved);
   
       return make_triplet_array( nsamples, saved_idx_v, 0, withtruth, nfilled );
   
     }
     
    
   }  
   }