FillkNN.cxx

Go to the documentation of this file.

// $Id: FillkNN.cxx,v 1.22 2009/09/22 18:32:34 jyuko Exp $

// C++
#include <cassert>
#include <climits>
#include <cmath>
#include <iterator>
#include <iomanip>
#include <sstream>

// ROOT
#include "TFile.h"
#include "TMath.h"

// MINOS
#include "Registry/Registry.h"
#include "Util/UtilString.h"

// Local
#include "PhysicsNtuple/Default.h"
#include "PhysicsNtuple/Factory.h"
#include "LikeModule.h"

// Local
#include "FillkNN.h"

REGISTER_ANP_OBJECT(AlgEvent,FillkNN)
REGISTER_ANP_OBJECT(AlgSnarl,FillkNN)

using namespace std;

//------------------------------------------------------------------------------------------
Anp::FillkNN::FillkNN()
   :fNAbortFind(0),
    fNTrackMiss(0),    
    fNTrackFill(0),
    fNEventMiss(0),
    fNEventFill(0),
    fDirName("fill_knn"),
    fFilePath(),
    fMaxNEvent(0),
    fKeyBase(5200),
    fKeySignal(1),
    fKeyTruth(7000),
    fTrimDelta(0),
    fKNeighbor(81),
    fKNeighborMod(20),
    fAddAll(false),
    fAddEvent(true),
    fCreateNew(false),
    fPrint(false),
    fRegression(false),
    fTrim(true),
    fUseKernel(false),
    fUseTrack(true),
    fModule(0)
{ 
}

//------------------------------------------------------------------------------------------
Anp::FillkNN::~FillkNN()
{
   if(fModule)
   {
      delete fModule;
      fModule = 0;
   }
}

//-----------------------------------------------------------------------------
bool Anp::FillkNN::Run(Record &record)
{
   //
   // Compute kNN variable for events and/or tracks
   //
   if(!fModule) return true;

   map<short, const DataVec> dmap;

   if(fUseTrack)
   {
      for(TrackIterator itrack = record.TrackBegIterator(); itrack != record.TrackEndIterator(); ++itrack)
      {

         if(fCreateNew)
         {
            if(FillkNN::Fill(itrack -> DataBeg(), itrack -> DataEnd(), itrack -> Weight()))
            {
               ++fNTrackFill;
            }
            else
            {
               ++fNTrackMiss;
            }
         }
         else
         {
            const DataVec dvec = FillkNN::Find(itrack -> DataBeg(), itrack -> DataEnd(), fKeyBase);
         
            if(dvec.empty())
            {
               ++fNTrackMiss;
            }
            else
            {
               ++fNTrackFill;
            }
            
            for(DataVec::const_iterator dit = dvec.begin(); dit != dvec.end(); ++dit)
            {
               itrack -> Add(dit -> Key(), dit -> Data());
            }
            
            if(!dmap.insert(map<short, const DataVec>::value_type(itrack -> TrackIndex(), dvec)).second)
            {
               cerr << "FillkNN::Run - duplicate track index " << itrack -> TrackIndex() << endl;
            }
         }
      }
   }

   for(EventIterator event = record.EventBegIterator(); event != record.EventEndIterator(); ++event)
   {
      if(fUseTrack)
      {
         if(!fCreateNew && fAddEvent)
         {
            const TrackIter track = Anp::LongestTrack(*event, record);
            if(track == record.TrackEnd())
            {
               continue;
            }
            
            const map<short, const DataVec>::const_iterator dit = dmap.find(track -> TrackIndex());
            if(dit == dmap.end())
            {
               cerr << "FillkNN::Run - failed to find track index " << track -> TrackIndex() << endl;
               continue;
            }
            
            const DataVec &dvec = dit -> second;
            
            for(DataVec::const_iterator dit = dvec.begin(); dit != dvec.end(); ++dit)
            {
               event -> Add(dit -> Key(), dit -> Data());
            }
         }
      }
      else if(fCreateNew)
      {
         if(FillkNN::Fill(event -> DataBeg(), event -> DataEnd(), event -> Weight()))
         {
            ++fNEventFill;
         }
         else
         {
            ++fNEventMiss;
         }
      }
      else
      {
         const DataVec dvec = FillkNN::Find(event -> DataBeg(), event -> DataEnd(), fKeyBase);
         
         if(dvec.empty())
         {
            ++fNEventMiss;
         }
         else
         {
            ++fNEventFill;
         }
         
         for(DataVec::const_iterator dit = dvec.begin(); dit != dvec.end(); ++dit)
         {
            event -> Add(dit -> Key(), dit -> Data());
         }
      }
   }   

   return true;
}

//-----------------------------------------------------------------------------
bool Anp::FillkNN::Run(Event &event, const Record &record, const bool pass)
{
   //
   // Compute kNN variable for events
   //
   if(!pass || !fModule) return true;

   if(fUseTrack)
   {
      const TrackIter track = Anp::LongestTrack(event, record);
      if(track == record.TrackEnd())
      {
         return true;
      }
      
      if(fCreateNew)
      {  
         if(FillkNN::Fill(track -> DataBeg(), track -> DataEnd(), event.Weight()))
         {
            ++fNTrackFill;
         }
         else
         {
            ++fNTrackMiss;  
         }
      }
      else
      {
         const DataVec dvec = FillkNN::Find(track -> DataBeg(), track -> DataEnd(), fKeyBase);
         
         if(dvec.empty())
         {
            ++fNTrackMiss;
         }
         else
         {
            ++fNTrackFill;
         }

         for(DataVec::const_iterator dit = dvec.begin(); dit != dvec.end(); ++dit)
         {
            event.Add(dit -> Key(), dit -> Data());
         }
      }
   }
   else
   {
      if(fCreateNew)
      {
         if(FillkNN::Fill(event.DataBeg(), event.DataEnd(), event.Weight()))
         {
            ++fNEventFill;
         }
         else
         {
            ++fNEventMiss;
         }
      }
      else
      {
         const DataVec dvec = FillkNN::Find(event.DataBeg(), event.DataEnd(), fKeyBase);
         
         if(dvec.empty())
         {
            ++fNEventMiss;
         }
         else
         {
            ++fNEventFill;
         }

         for(DataVec::const_iterator dit = dvec.begin(); dit != dvec.end(); ++dit)
         {
            event.Add(dit -> Key(), dit -> Data());
         }
      }
   }
   
   return true;
}

//-----------------------------------------------------------------------------
void Anp::FillkNN::Config(const Registry &reg)
{
   //
   // Configure self
   //

   fKeys = Anp::ReadList<short>(reg, "FillkNNKeyList", ", ");
   fTgts = Anp::ReadList<short>(reg, "FillkNNTgtList", ", ");

   //
   // Read TDirectory name from Registry OR create name from list of keys
   //
   Anp::Read(reg, "FillkNNFilePath", fFilePath);
   Anp::Read(reg, "FillkNNDirName",  fDirName);

   if((fDirName == "fill_knn" || fDirName.empty()) && !fKeys.empty())
   {
      //
      // Make directory name based on input list of keys
      //
      stringstream dname;
      dname << "knn";
      for(ShortVec::const_iterator it = fKeys.begin(); it != fKeys.end(); ++it)
      {
         dname << "_" << *it;
      }
      
      fDirName = dname.str();
   }

   //
   // Read "int" Registry keys, if exist, directly to internal variables
   //
   reg.Get("FillkNNMaxNEvent", fMaxNEvent);
   reg.Get("FillkNNKeyBase",   fKeyBase);
   reg.Get("FillkNNKeySignal", fKeySignal);
   reg.Get("FillkNNKeyTruth",  fKeyTruth);
   reg.Get("FillkNNTrimDelta", fTrimDelta);

   int value_int = 0;
   if(reg.Get("FillkNNKNeighbor", value_int) && value_int >= 0)
   {
      fKNeighbor = static_cast<unsigned int>(value_int);
   }

   value_int = 0;
   if(reg.Get("FillkNNKNeighborMod", value_int) && value_int > 0)
   {
      fKNeighborMod = static_cast<unsigned int>(value_int);
   }

   //
   // Read "bool" Registry keys, if exist, directly to internal variables
   //
   Anp::Read(reg, "FillkNNAddAll",     fAddAll);
   Anp::Read(reg, "FillkNNAddEvent",   fAddEvent);
   Anp::Read(reg, "FillkNN",           fCreateNew);
   Anp::Read(reg, "FillkNNPrint",      fPrint);
   Anp::Read(reg, "FillkNNRegression", fRegression);
   Anp::Read(reg, "FillkNNTrim",       fTrim);
   Anp::Read(reg, "FillkNNUseKernel",  fUseKernel);
   Anp::Read(reg, "FillkNNUseTrack",   fUseTrack);

   if(reg.KeyExists("PrintConfig"))
   {
      cout << "FillkNN::Config" << endl;

      Anp::PrintList<short>("   KeyList", fKeys);
      Anp::PrintList<short>("   TgtList", fTgts);

      cout << "   DirName = " << fDirName << endl
           << "   FilePath = " << fFilePath << endl
           << "   MaxNEvent = " << fMaxNEvent << endl
           << "   KeyBase = " << fKeyBase << endl
           << "   KeySignal = " << fKeySignal << endl
           << "   KeyTruth = " << fKeyTruth << endl
           << "   KNeighbor = " << fKNeighbor << endl
           << "   KNeighborMod = " << fKNeighborMod << endl
           << "   AddAll = " << fAddAll << endl
           << "   AddEvent = " << fAddEvent << endl
           << "   CreateNew = " << fCreateNew << endl
           << "   Print = " << fPrint << endl
           << "   Regression = " << fRegression << endl
           << "   Trim = " << fTrim << endl
           << "   TrimDelta = " << fTrimDelta << endl
           << "   UseKernel = " << fUseKernel << endl
           << "   UseTrack = " << fUseTrack << endl;
   }

   //
   // Delete previous kNN module
   //
   if(fModule)
   {
      delete fModule;
   }  
}

//------------------------------------------------------------------------------------------
bool Anp::FillkNN::Init(const Header &)
{
   //
   // Configure self and (maybe) create new kNN classification module
   //

   if(fModule)
   {
      cerr << "FillkNN::Init - knn module already exists" << endl;
      return true;
   }
   else
   {
      fModule = new Lit::LikeModule();
   }

   //
   // Do not read file when creating new knn prototype
   //
   if(fCreateNew) return true;

   cout << "FillkNN::Init - attempting to read prototype from:\n   "
        << fFilePath << ":" << fDirName << endl;

   TFile *file = TFile::Open(fFilePath.c_str(), "READ");   
   if(!file || !(file -> IsOpen()))
   {
      cerr << "FillkNN::Init - failed to open ROOT file" << endl;
      return false;
   }

   //
   // Get TDirectory from file
   //
   TDirectory *dir = dynamic_cast<TDirectory *>(file -> Get(fDirName.c_str()));
   
   if(dir && fModule -> Read(dir, fMaxNEvent))
   {
      // with metric option variables are rescaled to same width
      // 7 is number of kd-tree levels that are used to split
      // variables      

      //
      // Get keys for stored events (eg background=0 and signal=1, etc)
      //      
      const Lit::EventVec &evec = fModule -> GetEventVec();
      
      for(Lit::EventVec::const_iterator it = evec.begin(); it != evec.end(); ++it)
      {
         if(std::find(fTypes.begin(), fTypes.end(), it -> GetType()) == fTypes.end())
         {
            fTypes.push_back(it -> GetType());
         }
      }
      
      std::sort(fTypes.begin(), fTypes.end());
      
      for(unsigned int i = 0; i < fTypes.size(); ++i)
      {
         if(i == 0)
         {
            cout << "FillkNN::Init - prototypes: ";
         }
         if(i + 1 != fTypes.size())
         {
            cout << fTypes[i] << ", ";
         }
         else
         {
            cout << fTypes[i] << endl;
         }
      }

      if(fTrim)
      {
         fModule -> Fill(7, "metric trim erase");
      }
      else
      {
         fModule -> Fill(7, "metric erase");
      }
   }
   else
   {
      cerr << "FillkNN::Init - failed to read LikeModule" << endl;
      delete fModule; fModule = 0;      
   }

   //
   // Close ROOT file which was opened for reading
   //
   file -> Close();

   return true;
}

//------------------------------------------------------------------------------------------
void Anp::FillkNN::End(const DataBlock &)
{
   if(fModule && fCreateNew)
   {
      TFile *file = TFile::Open(fFilePath.c_str(), "UPDATE");
      
      if(file && file -> IsOpen())
      {
         // trim option counts number of events of each type
         // and trims each event class to have the same 
         // number of events
         // 7 is number of kd-tree levels that are filled with empty events
         // in order to optimize kd-tree search algorithm
         
         if(fTrim)
         {
            fModule -> Fill(7, "trim");
         }
         else
         {
            fModule -> Fill(7, "");
         }

         cout << "FillkNN::End - writing prototype to\n   " 
              << fFilePath << ":" << fDirName << endl;

         if(file -> Get(fDirName.c_str()))
         {
            const string delname = fDirName + ";1";
            file -> Delete(delname.c_str());        
            cout << "FillkNN::End - deleted TDirectory " << delname << endl;
         }

         fModule -> Write(file -> mkdir(fDirName.c_str()));

         file -> Write();
         file -> Close();
      }
   }

   if(fModule)
   {
      delete fModule;
      fModule = 0;
   }

   cout << "   NAbortFind = " << fNAbortFind << endl
        << "   NTrackMiss = " << fNTrackMiss << endl
        << "   NTrackFill = " << fNTrackFill << endl
        << "   NEventMiss = " << fNEventMiss << endl
        << "   NEventFill = " << fNEventFill << endl;
}

//-----------------------------------------------------------------------------
bool Anp::FillkNN::Fill(DataIter beg, DataIter end, const double weight)
{
   //
   // Add variables for new protype
   //

   if(!fModule)
   {
      cerr << "FillkNN::Fill - knn module does not exist" << endl;
      return false;
   }

   short etype = 0;
   if(!fRegression)
   {
      assert(fKeyTruth < SHRT_MAX && "int to short overflow");
      
      const DataIter fit = std::find(beg, end, static_cast<short>(fKeyTruth));
      if(fit == end)
      {
         return false;
      }
      
      assert(fit -> Data() < SHRT_MAX && "int to short overflow");
      etype = static_cast<short>(fit -> Data());

      if(fTrim && fTrimDelta > 0)
      {
         //
         // Get iterator to current event type
         //
         ShortMap::iterator cur_it = fCount.insert(ShortMap::value_type(etype, 0)).first;
      
         //
         // Find event type with minimum number of entries
         //
         ShortMap::const_iterator min_it = fCount.end();

         for(ShortMap::const_iterator cit = fCount.begin(); cit != fCount.end(); ++cit)
         {
            if(min_it == fCount.end())
            {
               min_it = cit;
            }
            else if(cit -> second < min_it -> second)
            {
               min_it = cit;
            }
         }
         
         if(min_it != fCount.end() && cur_it -> second > fTrimDelta + min_it -> second)
         {
            return false;
         }
         
         ++(cur_it -> second);
      }
   }
  
   const vector<float> dvec = FillkNN::Get(beg, end, fKeys);
   const vector<float> tvec = FillkNN::Get(beg, end, fTgts);

   if(dvec.size() != fKeys.size())
   {
      return false;
   }      

   fModule -> Add(Lit::Event(dvec, weight, etype, tvec));

   return true;
}

//-----------------------------------------------------------------------------
const Anp::DataVec Anp::FillkNN::Find(DataIter beg, DataIter end, const short keybase) const
{
   //
   // Find knn and fill classification variables
   //
   DataVec dvec;

   if(!fModule)
   {
      cerr << "FillkNN::Fill - knn module does not exist" << endl;
      return dvec;
   }

   //
   // !!!KLUDGE!!! Stop now if knn radius already exists
   //
   if(std::find(dvec.begin(), dvec.end(), keybase) != dvec.end())
   {
      ++fNAbortFind;
      return dvec;
   }

   //
   // Read all required input variables 
   //
   const vector<float> kvec = FillkNN::Get(beg, end, fKeys);
   if(kvec.size() != fKeys.size())
   {
      return dvec;
   }

   fModule -> Find(Lit::Event(kvec, 1.0, 0), fKNeighbor);
   
   if(fPrint)
   {
      cout << "Printing kNN result for key " << keybase << endl;
      fModule -> GetResult().Print();
   }
   
   const Lit::List &rlist = fModule -> GetResult().GetList();
   
   if(rlist.empty())
   {
      cerr << "FillkNN::Find - kNN result list is empty" << endl;
      return dvec;
   }

   if(fRegression)
   {
      double totalw = 0.0;
      unsigned int icount = 0;
      vector<double> tvec;

      for(Lit::List::const_iterator lit = rlist.begin(); lit != rlist.end(); ++lit)
      {
         const Lit::Event &event_ = lit->first->GetEvent();

         ++icount;
         totalw += event_.GetWeight();

         if(tvec.empty())
         {
            tvec.insert(tvec.end(), event_.GetNTgt(), 0.0);
         }

         if(event_.GetNTgt() < 1)
         {
            cerr << "FillkNN::Find - target list is empty" << endl;
         }
         
         for(unsigned int ivar = 0; ivar < event_.GetNTgt(); ++ivar)
         {
            tvec[ivar] += event_.GetTgt(ivar)*event_.GetWeight();
         }
      
         //
         // Continue to next event when modulus is not zero...
         //
         if(icount % fKNeighborMod != 0) continue;

         //
         // Check total weight sum
         //
         if(!(totalw > 0.0))
         {
            cerr << "FillkNN::Find - non positive total weight" << endl;
            continue;
         }
         
         if(fPrint)
         {
            cout << "FillkNN::Find - saving regression result for icount = " << icount << endl;
         }

         for(unsigned int ivar = 0; ivar < tvec.size(); ++ivar)
         {
            dvec.push_back(Anp::Data(keybase + icount + ivar, tvec[ivar]/totalw));
         }
      }
   }
   else
   {
      double totalw = 0.0;
      unsigned int count = 0;
      map<short, double> pmap;
      
      for(Lit::List::const_iterator lit = rlist.begin(); lit != rlist.end(); ++lit)
      {
         const Lit::Node &node = *(lit -> first);
         
         //
         // ignore Monte-Carlo event with zero distance (probably same event)
         //
         if(!(lit -> second > 0.0)) continue;
         
         //
         // Count weight for each truth category
         //
         map<short, double>::iterator pit = 
            pmap.insert(map<short, double>::value_type(node.GetType(), 0.0)).first;
         
         pit -> second += node.GetWeight();
         totalw += node.GetWeight();
         
         ++count;

         //
         // Continue to next event when modulus is not zero...
         //
         if(count % fKNeighborMod != 0) continue;
         
         if(!(totalw > 0.0))
         {
            cerr << "FillkNN::Find - non positive total weight" << endl;
            continue;
         }
         
         //
         // Add knn probability value 
         //
         for(ShortVec::const_iterator itype = fTypes.begin(); itype != fTypes.end(); ++itype)
         {
            const short type = *itype;         
            double prob = 0.0;
            
            map<short, double>::const_iterator pit = pmap.find(type);
            if(pit != pmap.end())
            {
               prob = (pit -> second)/totalw;
            }
            
            if(fAddAll)
            {
               dvec.push_back(Anp::Data(keybase + count + type, prob));
            }
            else if(type == fKeySignal)
            {
               dvec.push_back(Anp::Data(keybase + count, prob));
            }
         }
      }

      //
      // Add radius of knn neighborhood
      //
      dvec.push_back(Anp::Data(keybase, std::sqrt(rlist.back().second)));
   }

   return dvec;
}

//------------------------------------------------------------------------------------------
double Anp::FillkNN:: PidKer(Lit::List::const_iterator beg, Lit::List::const_iterator end) const
{
   //
   // 
   //
   double maxradius = -1.0;

   for(Lit::List::const_iterator lit = beg; lit != end; ++lit)
   {
      if(maxradius < lit -> second)
      {
         maxradius = lit -> second;
      }
   }

   if(!(maxradius > 0.0))
   {
      cerr << "FillkNN::PidKer - maximum radius is not positive" << endl;
      return -100.0;
   }

   maxradius = 1.0/std::sqrt(maxradius);

   double count_signal = 0.0, count_all = 0.0;

   for(Lit::List::const_iterator lit = beg; lit != end; ++lit)
   {
      const Lit::Node &node = *(lit -> first);

      const double weight = Kernel(std::sqrt(lit -> second)*maxradius) * node.GetWeight();

      if(weight < 0.0)
      {
         cerr << "FillkNN::PidKer - negative event weight " << weight << endl;
         continue;
      }
      
      count_all += weight;

      if(node.GetType() == fKeySignal)
      {
         count_signal += weight;
      }
   }

   if(count_all > 0.0)
   {
      return count_signal/count_all;
   }

   cerr << "FillkNN:: PidKer - total event count is not positive" << endl;

   return -100.0;
}

//------------------------------------------------------------------------------------------
const vector<float> Anp::FillkNN::Get(DataIter beg, DataIter end, const ShortVec &kvec) const
{
   // extract data with keys matching keys stored in ShortVec kvec
   // these data elements are first stored in a map to make sure
   // that data is sorted using keys - thus removing the dependency
   // on the order in which keys are stored in DataVec, with range
   // delimited by beg and end iterators

   map<short, float> fmap;
   
   for(DataIter dit = beg; dit != end; ++dit)
   {
      if(find(kvec.begin(), kvec.end(), dit -> Key()) != kvec.end())
      {
         if(!fmap.insert(map<short, float>::value_type(dit -> Key(), dit -> Data())).second)
         {
            cerr << "FillkNN::Get - failed to add data due to duplicate key" << endl;
         }
      }
   }
   
   if(fPrint && !kvec.empty() && fmap.size() != kvec.size())
   {
      cout << "FillkNN::Get - failed to read all keys" << endl;
      
      for(unsigned int i = 0; i < kvec.size(); ++i)
      {
         map<short, float>::const_iterator fit = fmap.find(kvec[i]);
         if(fit != fmap.end())
         {
            cout << "   (key, data) = (" << kvec[i] << ", " << fit -> second << ")" << endl;
         }
         else
         {
            cout << "   (key, data) = (" << kvec[i] << ", unknown)" << endl;
         }
      }
   }

   vector<float> fvec;

   for(map<short, float>::const_iterator it = fmap.begin(); it != fmap.end(); ++it)
   {
      fvec.push_back(it -> second);
   }

   return fvec;
}

//-----------------------------------------------------------------------------------------
double Anp::FillkNN::Kernel(const double value) const
{
   //
   // return distance using unit normalize kernel: 15*(1-x*x)^2
   //

   const double avalue = std::fabs(value);

   if(!(avalue < 1.0))
   {
      return 0.0;
   }  

   const double prod = 1.0 - avalue*avalue;

   return 15.0*prod*prod/16.0;
}

---