TransformationHandler.cxx

Go to the documentation of this file.

// @(#)root/tmva $Id$
// Author: Andreas Hoecker, Peter Speckmayer, Joerg Stelzer, Helge Voss, Eckhard von Toerne, Jan Therhaag

/**********************************************************************************
 * Project: TMVA - a Root-integrated toolkit for multivariate data analysis       *
 * Package: TMVA                                                                  *
 * Class  : TransformationHandler                                                 *
 * Web    : http://tmva.sourceforge.net                                           *
 *                                                                                *
 * Description:                                                                   *
 *      Implementation (see header for description)                               *
 *                                                                                *
 * Authors (alphabetical):                                                        *
 *      Andreas Hoecker <Andreas.Hocker@cern.ch> - CERN, Switzerland              *
 *      Peter Speckmayer <speckmay@mail.cern.ch>  - CERN, Switzerland             *
 *      Joerg Stelzer   <Joerg.Stelzer@cern.ch>  - CERN, Switzerland              *
 *      Jan Therhaag       <Jan.Therhaag@cern.ch>     - U of Bonn, Germany        *
 *      Eckhard v. Toerne  <evt@uni-bonn.de>     - U of Bonn, Germany             *
 *      Helge Voss      <Helge.Voss@cern.ch>     - MPI-K Heidelberg, Germany      *
 *                                                                                *
 * Copyright (c) 2005-2011:                                                       *
 *      CERN, Switzerland                                                         *
 *      MPI-K Heidelberg, Germany                                                 *
 *      U. of Bonn, Germany                                                       *
 *                                                                                *
 * Redistribution and use in source and binary forms, with or without             *
 * modification, are permitted according to the terms listed in LICENSE           *
 * (http://tmva.sourceforge.net/LICENSE)                                          *
 **********************************************************************************/

/*! \class TMVA::TransformationHandler
\ingroup TMVA
Class that contains all the data information.
*/

#include "TMVA/TransformationHandler.h"

#include "TMVA/Config.h"
#include "TMVA/DataSet.h"
#include "TMVA/DataSetInfo.h"
#include "TMVA/Event.h"
#include "TMVA/MsgLogger.h"
#include "TMVA/Ranking.h"
#include "TMVA/Tools.h"
#include "TMVA/Types.h"
#include "TMVA/VariableDecorrTransform.h"
#include "TMVA/VariableGaussTransform.h"
#include "TMVA/VariableInfo.h"
#include "TMVA/VariableIdentityTransform.h"
#include "TMVA/VariableNormalizeTransform.h"
#include "TMVA/VariablePCATransform.h"
#include "TMVA/VariableRearrangeTransform.h"
#include "TMVA/VariableTransformBase.h"

#include "TAxis.h"
#include "TDirectory.h"
#include "TH1.h"
#include "TH2.h"
#include "TList.h"
#include "TMath.h"
#include "TProfile.h"

#include <vector>
#include <iomanip>

////////////////////////////////////////////////////////////////////////////////
/// constructor

TMVA::TransformationHandler::TransformationHandler( DataSetInfo& dsi, const TString& callerName )
   : fDataSetInfo(dsi),
     fRootBaseDir(0),
     fCallerName (callerName),
     fLogger     ( new MsgLogger(TString("TFHandler_" + callerName).Data(), kINFO) )
{
   // produce one entry for each class and one entry for all classes. If there is only one class,
   // produce only one entry
   fNumC = (dsi.GetNClasses()<= 1) ? 1 : dsi.GetNClasses()+1;

   fVariableStats.resize( fNumC );
   for (Int_t i=0; i<fNumC; i++ ) fVariableStats.at(i).resize(dsi.GetNVariables() + dsi.GetNTargets());
}

////////////////////////////////////////////////////////////////////////////////
/// destructor

TMVA::TransformationHandler::~TransformationHandler()
{
   std::vector<Ranking*>::const_iterator it = fRanking.begin();
   for (; it != fRanking.end(); ++it) delete *it;

   fTransformations.SetOwner();
   delete fLogger;
}

////////////////////////////////////////////////////////////////////////////////

void TMVA::TransformationHandler::SetCallerName( const TString& name )
{
   fCallerName = name;
   fLogger->SetSource( TString("TFHandler_" + fCallerName).Data() );
}

////////////////////////////////////////////////////////////////////////////////

TMVA::VariableTransformBase* TMVA::TransformationHandler::AddTransformation( VariableTransformBase *trf, Int_t cls )
{
   TString tfname = trf->Log().GetName();
   trf->Log().SetSource(TString(fCallerName+"_"+tfname+"_TF").Data());
   fTransformations.Add(trf);
   fTransformationsReferenceClasses.push_back( cls );
   return trf;
}

////////////////////////////////////////////////////////////////////////////////
/// Caches calculated summary statistics of transformed variables.
///
/// \param[in] k index of class
/// \param[in] ivar index of variable
/// \param[in] mean the mean value of the variable
/// \param[in] rms the root-mean-square value of the variable
/// \param[in] min the minimum value of the variable
/// \param[in] max the maximum value of the variable

void TMVA::TransformationHandler::AddStats( Int_t k, UInt_t ivar, Double_t mean, Double_t rms, Double_t min, Double_t max )
{
   if (rms <= 0 || TMath::IsNaN(rms)) {
      Log() << kWARNING << "Variable \"" << Variable(ivar).GetExpression()
            << "\" has zero, negative, or NaN RMS^2: "
            << rms
            << " ==> set to zero. Please check the variable content" << Endl;
      rms = 0;
   }

   VariableStat stat; stat.fMean = mean; stat.fRMS = rms; stat.fMin = min; stat.fMax = max;
   fVariableStats.at(k).at(ivar) = stat;
}

////////////////////////////////////////////////////////////////////////////////
/// overrides the setting for all classes! (this is put in basically for the likelihood-method)
/// be careful with the usage this method

void TMVA::TransformationHandler::SetTransformationReferenceClass( Int_t cls )
{
   for (UInt_t i = 0; i < fTransformationsReferenceClasses.size(); i++) {
      fTransformationsReferenceClasses.at( i ) = cls;
   }
}

////////////////////////////////////////////////////////////////////////////////
/// the transformation

const TMVA::Event* TMVA::TransformationHandler::Transform( const Event* ev ) const
{
   TListIter trIt(&fTransformations);
   std::vector<Int_t>::const_iterator rClsIt = fTransformationsReferenceClasses.begin();
   const Event* trEv = ev;
   while (VariableTransformBase *trf = (VariableTransformBase*) trIt()) {
      if (rClsIt == fTransformationsReferenceClasses.end()) Log() << kFATAL<< "invalid read in TransformationHandler::Transform " <<Endl;
      trEv = trf->Transform(trEv, (*rClsIt) );
      ++rClsIt;
   }
   return trEv;
}

////////////////////////////////////////////////////////////////////////////////

const TMVA::Event* TMVA::TransformationHandler::InverseTransform( const Event* ev, Bool_t suppressIfNoTargets ) const
{
   if (fTransformationsReferenceClasses.empty()){
      //Log() << kWARNING << __FILE__ <<":InverseTransform fTransformationsReferenceClasses is empty" << Endl;
      return ev;
   }
   // the inverse transformation
   TListIter trIt(&fTransformations, kIterBackward);
   std::vector< Int_t >::const_iterator rClsIt = fTransformationsReferenceClasses.end();
   --rClsIt;
   const Event* trEv = ev;
   UInt_t nvars = 0, ntgts = 0, nspcts = 0;
   while (VariableTransformBase *trf = (VariableTransformBase*) trIt() ) { // shouldn't be the transformation called in the inverse order for the inversetransformation?????
      if (trf->IsCreated()) {
         trf->CountVariableTypes( nvars, ntgts, nspcts );
         if( !(suppressIfNoTargets && ntgts==0) )
            trEv = trf->InverseTransform(ev, (*rClsIt) );
      }
      else break;
      --rClsIt;
   }
   return trEv;


   //    TListIter trIt(&fTransformations);
   //    std::vector< Int_t >::const_iterator rClsIt = fTransformationsReferenceClasses.begin();
   //    const Event* trEv = ev;
   //    UInt_t nvars = 0, ntgts = 0, nspcts = 0;
   //    while (VariableTransformBase *trf = (VariableTransformBase*) trIt() ) { // shouldn't be the transformation called in the inverse order for the inversetransformation?????
   //       if (trf->IsCreated()) {
   //          trf->CountVariableTypes( nvars, ntgts, nspcts );
   //    if( !(suppressIfNoTargets && ntgts==0) )
   //       trEv = trf->InverseTransform(ev, (*rClsIt) );
   //       }
   //       else break;
   //       rClsIt++;
   //    }
   //    return trEv;

}

////////////////////////////////////////////////////////////////////////////////
/// computation of transformation

const std::vector<TMVA::Event*>* TMVA::TransformationHandler::CalcTransformations( const std::vector<Event*>& events,
                                                                                   Bool_t createNewVector )
{
   if (fTransformations.GetEntries() <= 0)
      return &events;

   // the transformedEvents are initialised with the initial events and then
   // subsequently replaced with transformed ones. The n-th transformation will
   // and on the events as they look like after the (n-1)-the transformation
   // as intended for the chained transformations
   std::vector<Event*> *transformedEvents = new std::vector<TMVA::Event*>(events.size());
   for ( UInt_t ievt = 0; ievt<events.size(); ievt++)
      transformedEvents->at(ievt) = new Event(*events.at(ievt));

   TListIter trIt(&fTransformations);
   std::vector< Int_t >::iterator rClsIt = fTransformationsReferenceClasses.begin();
   while (VariableTransformBase *trf = (VariableTransformBase*) trIt()) {
      if (trf->PrepareTransformation(*transformedEvents)) {
         for (UInt_t ievt = 0; ievt<transformedEvents->size(); ievt++) {  // loop through all events
            *(*transformedEvents)[ievt] = *trf->Transform((*transformedEvents)[ievt],(*rClsIt));
         }
         ++rClsIt;
      }
   }

   CalcStats(*transformedEvents);

   // plot the variables once in this transformation
   PlotVariables(*transformedEvents);

   //sometimes, the actual transformed event vector is not used for anything but the previous
   //CalcStat and PlotVariables calles, in that case, we delete it again (and return NULL)
   if (!createNewVector) {  // if we don't want that newly created event vector to persist, then delete it
      for ( UInt_t ievt = 0; ievt<transformedEvents->size(); ievt++)
         delete (*transformedEvents)[ievt];
      delete transformedEvents;
      transformedEvents=NULL;
   }

   return transformedEvents; // give back the newly created event collection (containing the transformed events)
}

////////////////////////////////////////////////////////////////////////////////
/// method to calculate minimum, maximum, mean, and RMS for all
/// variables used in the MVA

void TMVA::TransformationHandler::CalcStats (const std::vector<Event*>& events )
{
   UInt_t nevts = events.size();

   if (nevts==0)
      Log() << kFATAL << "No events available to find min, max, mean and rms" << Endl;

   // if transformation has not been succeeded, the tree may be empty
   const UInt_t nvar = events[0]->GetNVariables();
   const UInt_t ntgt = events[0]->GetNTargets();

   Double_t  *sumOfWeights = new Double_t[fNumC];
   Double_t* *x2           = new Double_t*[fNumC];
   Double_t* *x0           = new Double_t*[fNumC];
   Double_t* *varMin       = new Double_t*[fNumC];
   Double_t* *varMax       = new Double_t*[fNumC];

   for (Int_t cls=0; cls<fNumC; cls++) {
      sumOfWeights[cls]=0;
      x2[cls]     = new Double_t[nvar+ntgt];
      x0[cls]     = new Double_t[nvar+ntgt];
      varMin[cls] = new Double_t[nvar+ntgt];
      varMax[cls] = new Double_t[nvar+ntgt];
      for (UInt_t ivar=0; ivar<nvar+ntgt; ivar++) {
         x0[cls][ivar] = x2[cls][ivar] = 0;
         varMin[cls][ivar] = DBL_MAX;
         varMax[cls][ivar] = -DBL_MAX;
      }
   }

   for (UInt_t ievt=0; ievt<nevts; ievt++) {
      const Event* ev  = events[ievt];
      Int_t  cls = ev->GetClass();

      Double_t weight = ev->GetWeight();
      sumOfWeights[cls] += weight;
      if (fNumC > 1 ) sumOfWeights[fNumC-1] += weight; // if more than one class, store values for all classes
      for (UInt_t var_tgt = 0; var_tgt < 2; var_tgt++ ){ // first for variables, then for targets
         UInt_t nloop = ( var_tgt==0?nvar:ntgt );
         for (UInt_t ivar=0; ivar<nloop; ivar++) {
            Double_t x = ( var_tgt==0?ev->GetValue(ivar):ev->GetTarget(ivar) );

            if (x < varMin[cls][(var_tgt*nvar)+ivar]) varMin[cls][(var_tgt*nvar)+ivar]= x;
            if (x > varMax[cls][(var_tgt*nvar)+ivar]) varMax[cls][(var_tgt*nvar)+ivar]= x;

            x0[cls][(var_tgt*nvar)+ivar] += x*weight;
            x2[cls][(var_tgt*nvar)+ivar] += x*x*weight;

            if (fNumC > 1) {
               if (x < varMin[fNumC-1][(var_tgt*nvar)+ivar]) varMin[fNumC-1][(var_tgt*nvar)+ivar]= x;
               if (x > varMax[fNumC-1][(var_tgt*nvar)+ivar]) varMax[fNumC-1][(var_tgt*nvar)+ivar]= x;

               x0[fNumC-1][(var_tgt*nvar)+ivar] += x*weight;
               x2[fNumC-1][(var_tgt*nvar)+ivar] += x*x*weight;
            }
         }
      }
   }


   // set Mean and RMS
   for (UInt_t var_tgt = 0; var_tgt < 2; var_tgt++ ){ // first for variables, then for targets
      UInt_t nloop = ( var_tgt==0?nvar:ntgt );
      for (UInt_t ivar=0; ivar<nloop; ivar++) {
         for (Int_t cls = 0; cls < fNumC; cls++) {
            Double_t mean = x0[cls][(var_tgt*nvar)+ivar]/sumOfWeights[cls];
            Double_t rms = TMath::Sqrt( x2[cls][(var_tgt*nvar)+ivar]/sumOfWeights[cls] - mean*mean);
            AddStats(cls, (var_tgt*nvar)+ivar, mean, rms, varMin[cls][(var_tgt*nvar)+ivar], varMax[cls][(var_tgt*nvar)+ivar]);
         }
      }
   }

   // ------ pretty output of basic statistics -------------------------------
   // find maximum length in V (and column title)
   UInt_t maxL = 8, maxV = 0;
   std::vector<UInt_t> vLengths;
   for (UInt_t ivar=0; ivar<nvar+ntgt; ivar++) {
      if( ivar < nvar )
         maxL = TMath::Max( (UInt_t)Variable(ivar).GetLabel().Length(), maxL );
      else
         maxL = TMath::Max( (UInt_t)Target(ivar-nvar).GetLabel().Length(), maxL );
   }
   maxV = maxL + 2;
   // full column length
   UInt_t clen = maxL + 4*maxV + 11;
   Log() << kHEADER ;
   //for (UInt_t i=0; i<clen; i++) //Log() << "-";

   //Log() << Endl;
   // full column length
   Log() << std::setw(maxL) << "Variable";
   Log() << "  " << std::setw(maxV) << "Mean";
   Log() << " " << std::setw(maxV) << "RMS";
   Log() << "   " << std::setw(maxV) << "[        Min ";
   Log() << "  " << std::setw(maxV) << "    Max ]"<< Endl;;
   for (UInt_t i=0; i<clen; i++) Log() << "-";
   Log() << Endl;

   // the numbers
   TString format = "%#11.5g";
   for (UInt_t ivar=0; ivar<nvar+ntgt; ivar++) {
      if( ivar < nvar )
         Log() << std::setw(maxL) << Variable(ivar).GetLabel() << ":";
      else
         Log() << std::setw(maxL) << Target(ivar-nvar).GetLabel() << ":";
      Log() << std::setw(maxV) << Form( format.Data(), GetMean(ivar) );
      Log() << std::setw(maxV) << Form( format.Data(), GetRMS(ivar) );
      Log() << "   [" << std::setw(maxV) << Form( format.Data(), GetMin(ivar) );
      Log() << std::setw(maxV) << Form( format.Data(), GetMax(ivar) ) << " ]";
      Log() << Endl;
   }
   for (UInt_t i=0; i<clen; i++) Log() << "-";
   Log() << Endl;
   // ------------------------------------------------------------------------

   delete[] sumOfWeights;
   for (Int_t cls=0; cls<fNumC; cls++) {
      delete [] x2[cls];
      delete [] x0[cls];
      delete [] varMin[cls];
      delete [] varMax[cls];
   }
   delete [] x2;
   delete [] x0;
   delete [] varMin;
   delete [] varMax;
}

////////////////////////////////////////////////////////////////////////////////
/// create transformation function

void TMVA::TransformationHandler::MakeFunction( std::ostream& fout, const TString& fncName, Int_t part ) const
{
   TListIter trIt(&fTransformations);
   std::vector< Int_t >::const_iterator rClsIt = fTransformationsReferenceClasses.begin();
   UInt_t trCounter=1;
   while (VariableTransformBase *trf = (VariableTransformBase*) trIt() ) {
      trf->MakeFunction(fout, fncName, part, trCounter++, (*rClsIt) );
      ++rClsIt;
   }
   if (part==1) {
      for (Int_t i=0; i<fTransformations.GetSize(); i++) {
         fout << "   void InitTransform_"<<i+1<<"();" << std::endl;
         fout << "   void Transform_"<<i+1<<"( std::vector<double> & iv, int sigOrBgd ) const;" << std::endl;
      }
   }
   if (part==2) {
      fout << std::endl;
      fout << "//_______________________________________________________________________" << std::endl;
      fout << "inline void " << fncName << "::InitTransform()" << std::endl;
      fout << "{" << std::endl;
      for (Int_t i=0; i<fTransformations.GetSize(); i++)
         fout << "   InitTransform_"<<i+1<<"();" << std::endl;
      fout << "}" << std::endl;
      fout << std::endl;
      fout << "//_______________________________________________________________________" << std::endl;
      fout << "inline void " << fncName << "::Transform( std::vector<double>& iv, int sigOrBgd ) const" << std::endl;
      fout << "{" << std::endl;
      for (Int_t i=0; i<fTransformations.GetSize(); i++)
         fout << "   Transform_"<<i+1<<"( iv, sigOrBgd );" << std::endl;

      fout << "}" << std::endl;
   }
}

////////////////////////////////////////////////////////////////////////////////
/// return transformation name

TString TMVA::TransformationHandler::GetName() const
{
   TString name("Id");
   TListIter trIt(&fTransformations);
   VariableTransformBase *trf;
   if ((trf = (VariableTransformBase*) trIt())) {
      name = TString(trf->GetShortName());
      while ((trf = (VariableTransformBase*) trIt())) name += "_" + TString(trf->GetShortName());
   }
   return name;
}

////////////////////////////////////////////////////////////////////////////////
/// incorporates transformation type into title axis (usually for histograms)

TString TMVA::TransformationHandler::GetVariableAxisTitle( const VariableInfo& info ) const
{
   TString xtit = info.GetTitle();
   // indicate transformation, but not in case of single identity transform
   if (fTransformations.GetSize() >= 1) {
      if (fTransformations.GetSize() > 1 ||
          ((VariableTransformBase*)GetTransformationList().Last())->GetVariableTransform() != Types::kIdentity) {
         xtit += " (" + GetName() + ")";
      }
   }
   return xtit;
}


////////////////////////////////////////////////////////////////////////////////
/// create histograms from the input variables
///  - histograms for all input variables
///  - scatter plots for all pairs of input variables

void TMVA::TransformationHandler::PlotVariables (const std::vector<Event*>& events, TDirectory* theDirectory )
{
   if (fRootBaseDir==0 && theDirectory == 0) return;

   Log() << kDEBUG << "Plot event variables for ";
   if (theDirectory !=0) Log()<< TString(theDirectory->GetName()) << Endl;
   else Log() << GetName() << Endl;

   // extension for transformation type
   TString transfType = "";
   if (theDirectory == 0) {
      transfType += "_";
      transfType += GetName();
   }else{ // you plot for the individual classifiers. Note, here the "statistics" still need to be calculated as you are in the testing phase
      CalcStats(events);
   }

   const UInt_t nvar = fDataSetInfo.GetNVariables();
   const UInt_t ntgt = fDataSetInfo.GetNTargets();
   const Int_t  ncls = fDataSetInfo.GetNClasses();

   // Create all histograms
   // do both, scatter and profile plots
   std::vector<std::vector<TH1*> > hVars( ncls );  // histograms for variables
   std::vector<std::vector<std::vector<TH2F*> > >     mycorr( ncls ); // histograms for correlations
   std::vector<std::vector<std::vector<TProfile*> > > myprof( ncls ); // histograms for profiles

   for (Int_t cls = 0; cls < ncls; cls++) {
      hVars.at(cls).resize ( nvar+ntgt );
      hVars.at(cls).assign ( nvar+ntgt, 0 ); // fill with zeros
      mycorr.at(cls).resize( nvar+ntgt );
      myprof.at(cls).resize( nvar+ntgt );
      for (UInt_t ivar=0; ivar < nvar+ntgt; ivar++) {
         mycorr.at(cls).at(ivar).resize( nvar+ntgt );
         myprof.at(cls).at(ivar).resize( nvar+ntgt );
         mycorr.at(cls).at(ivar).assign( nvar+ntgt, 0 ); // fill with zeros
         myprof.at(cls).at(ivar).assign( nvar+ntgt, 0 ); // fill with zeros
      }
   }

   // if there are too many input variables, the creation of correlations plots blows up
   // memory and basically kills the TMVA execution
   // --> avoid above critical number (which can be user defined)
   if (nvar+ntgt > (UInt_t)gConfig().GetVariablePlotting().fMaxNumOfAllowedVariablesForScatterPlots) {
      Int_t nhists = (nvar+ntgt)*(nvar+ntgt - 1)/2;
      Log() << kINFO << gTools().Color("dgreen") << Endl;
      Log() << kINFO << "<PlotVariables> Will not produce scatter plots ==> " << Endl;
      Log() << kINFO
            << "|  The number of " << nvar << " input variables and " << ntgt << " target values would require "
            << nhists << " two-dimensional" << Endl;
      Log() << kINFO
            << "|  histograms, which would occupy the computer's memory. Note that this" << Endl;
      Log() << kINFO
            << "|  suppression does not have any consequences for your analysis, other" << Endl;
      Log() << kINFO
            << "|  than not disposing of these scatter plots. You can modify the maximum" << Endl;
      Log() << kINFO
            << "|  number of input variables allowed to generate scatter plots in your" << Endl;
      Log() << "|  script via the command line:" << Endl;
      Log() << kINFO
            << "|  \"(TMVA::gConfig().GetVariablePlotting()).fMaxNumOfAllowedVariablesForScatterPlots = <some int>;\""
            << gTools().Color("reset") << Endl;
      Log() << Endl;
      Log() << kINFO << "Some more output" << Endl;
   }

   Double_t timesRMS = gConfig().GetVariablePlotting().fTimesRMS;
   UInt_t   nbins1D  = gConfig().GetVariablePlotting().fNbins1D;
   UInt_t   nbins2D  = gConfig().GetVariablePlotting().fNbins2D;

   for (UInt_t var_tgt = 0; var_tgt < 2; var_tgt++) { // create the histos first for the variables, then for the targets
      UInt_t nloops = ( var_tgt == 0? nvar:ntgt );     // number of variables or number of targets
      for (UInt_t ivar=0; ivar<nloops; ivar++) {
         const VariableInfo& info = ( var_tgt == 0 ? Variable( ivar ) : Target(ivar) ); // choose the appropriate one (variable or target)
         TString myVari = info.GetInternalName();

         Double_t mean = fVariableStats.at(fNumC-1).at( ( var_tgt*nvar )+ivar).fMean;
         Double_t rms  = fVariableStats.at(fNumC-1).at( ( var_tgt*nvar )+ivar).fRMS;

         for (Int_t cls = 0; cls < ncls; cls++) {

            TString className = fDataSetInfo.GetClassInfo(cls)->GetName();

            // add "target" in case of target variable (required for plotting macros)
            className += (ntgt == 1 && var_tgt == 1 ? "_target" : "");

            // choose reasonable histogram ranges, by removing outliers
            TH1* h = 0;
            if (info.GetVarType() == 'I') {
               // special treatment for integer variables
               Int_t xmin = TMath::Nint( GetMin( ( var_tgt*nvar )+ivar) );
               Int_t xmax = TMath::Nint( GetMax( ( var_tgt*nvar )+ivar) + 1 );
               Int_t nbins = xmax - xmin;

               h = new TH1F( Form("%s__%s%s", myVari.Data(), className.Data(), transfType.Data()),
                             info.GetTitle(), nbins, xmin, xmax );
            }
            else {
               Double_t xmin = TMath::Max( GetMin( ( var_tgt*nvar )+ivar), mean - timesRMS*rms );
               Double_t xmax = TMath::Min( GetMax( ( var_tgt*nvar )+ivar), mean + timesRMS*rms );

               //std::cout << "Class="<<cls<<" xmin="<<xmin << " xmax="<<xmax<<" mean="<<mean<<" rms="<<rms<<" timesRMS="<<timesRMS<<std::endl;
               // protection
               if (xmin >= xmax) xmax = xmin*1.1; // try first...
               if (xmin >= xmax) xmax = xmin + 1; // this if xmin == xmax == 0
               // safety margin for values equal to the maximum within the histogram
               xmax += (xmax - xmin)/nbins1D;

               h = new TH1F( Form("%s__%s%s", myVari.Data(), className.Data(), transfType.Data()),
                             info.GetTitle(), nbins1D, xmin, xmax );
            }

            h->GetXaxis()->SetTitle( gTools().GetXTitleWithUnit( GetVariableAxisTitle( info ), info.GetUnit() ) );
            h->GetYaxis()->SetTitle( gTools().GetYTitleWithUnit( *h, info.GetUnit(), kFALSE ) );
            hVars.at(cls).at((var_tgt*nvar)+ivar) = h;

            // profile and scatter plots
            if (nvar+ntgt <= (UInt_t)gConfig().GetVariablePlotting().fMaxNumOfAllowedVariablesForScatterPlots) {

               for (UInt_t v_t = 0; v_t < 2; v_t++) {
                  UInt_t nl = ( v_t==0?nvar:ntgt );
                  UInt_t start = ( v_t==0? (var_tgt==0?ivar+1:0):(var_tgt==0?nl:ivar+1) );
                  for (UInt_t j=start; j<nl; j++) {
                     // choose the appropriate one (variable or target)
                     const VariableInfo& infoj = ( v_t == 0 ? Variable( j ) : Target(j) );
                     TString myVarj = infoj.GetInternalName();

                     Double_t rxmin = fVariableStats.at(fNumC-1).at( ( v_t*nvar )+ivar).fMin;
                     Double_t rxmax = fVariableStats.at(fNumC-1).at( ( v_t*nvar )+ivar).fMax;
                     Double_t rymin = fVariableStats.at(fNumC-1).at( ( v_t*nvar )+j).fMin;
                     Double_t rymax = fVariableStats.at(fNumC-1).at( ( v_t*nvar )+j).fMax;

                     // scatter plot
                     TH2F* h2 = new TH2F( Form( "scat_%s_vs_%s_%s%s" , myVarj.Data(), myVari.Data(),
                                                className.Data(), transfType.Data() ),
                                          Form( "%s versus %s (%s)%s", infoj.GetTitle(), info.GetTitle(),
                                                className.Data(), transfType.Data() ),
                                          nbins2D, rxmin , rxmax,
                                          nbins2D, rymin , rymax );

                     h2->GetXaxis()->SetTitle( gTools().GetXTitleWithUnit( GetVariableAxisTitle( info  ), info .GetUnit() ) );
                     h2->GetYaxis()->SetTitle( gTools().GetXTitleWithUnit( GetVariableAxisTitle( infoj ), infoj.GetUnit() ) );
                     mycorr.at(cls).at((var_tgt*nvar)+ivar).at((v_t*nvar)+j) = h2;

                     // profile plot
                     TProfile* p = new TProfile( Form( "prof_%s_vs_%s_%s%s", myVarj.Data(),
                                                       myVari.Data(), className.Data(),
                                                       transfType.Data() ),
                                                 Form( "profile %s versus %s (%s)%s",
                                                       infoj.GetTitle(), info.GetTitle(),
                                                       className.Data(), transfType.Data() ), nbins1D,
                                                 rxmin, rxmax );
                     //                                                 info.GetMin(), info.GetMax() );

                     p->GetXaxis()->SetTitle( gTools().GetXTitleWithUnit( GetVariableAxisTitle( info  ), info .GetUnit() ) );
                     p->GetYaxis()->SetTitle( gTools().GetXTitleWithUnit( GetVariableAxisTitle( infoj ), infoj.GetUnit() ) );
                     myprof.at(cls).at((var_tgt*nvar)+ivar).at((v_t*nvar)+j) = p;
                  }
               }
            }
         }
      }
   }

   UInt_t nevts = events.size();

   // compute correlation coefficient between target value and variables (regression only)
   std::vector<Double_t> xregmean ( nvar+1, 0 );
   std::vector<Double_t> x2regmean( nvar+1, 0 );
   std::vector<Double_t> xCregmean( nvar+1, 0 );

   // fill the histograms (this approach should be faster than individual projection
   for (UInt_t ievt=0; ievt<nevts; ievt++) {

      const Event* ev = events[ievt];

      Float_t weight = ev->GetWeight();
      Int_t   cls    = ev->GetClass();

      // average correlation between first target and variables (so far only for single-target regression)
      if (ntgt == 1) {
         Float_t valr = ev->GetTarget(0);
         xregmean[nvar]  += valr;
         x2regmean[nvar] += valr*valr;
         for (UInt_t ivar=0; ivar<nvar; ivar++) {
            Float_t vali = ev->GetValue(ivar);
            xregmean[ivar]  += vali;
            x2regmean[ivar] += vali*vali;
            xCregmean[ivar] += vali*valr;
         }
      }

      // fill correlation histograms
      for (UInt_t var_tgt = 0; var_tgt < 2; var_tgt++) { // create the histos first for the variables, then for the targets
         UInt_t nloops = ( var_tgt == 0? nvar:ntgt );    // number of variables or number of targets
         for (UInt_t ivar=0; ivar<nloops; ivar++) {
            Float_t vali = ( var_tgt == 0 ? ev->GetValue(ivar) : ev->GetTarget(ivar) );

            // variable histos
            hVars.at(cls).at( ( var_tgt*nvar )+ivar)->Fill( vali, weight );

            // correlation histos
            if (nvar+ntgt <= (UInt_t)gConfig().GetVariablePlotting().fMaxNumOfAllowedVariablesForScatterPlots) {

               for (UInt_t v_t = 0; v_t < 2; v_t++) {
                  UInt_t nl    = ( v_t==0 ? nvar : ntgt );
                  UInt_t start = ( v_t==0 ? (var_tgt==0?ivar+1:0) : (var_tgt==0?nl:ivar+1) );
                  for (UInt_t j=start; j<nl; j++) {
                     Float_t valj = ( v_t == 0 ? ev->GetValue(j) : ev->GetTarget(j) );
                     mycorr.at(cls).at( ( var_tgt*nvar )+ivar).at( ( v_t*nvar )+j)->Fill( vali, valj, weight );
                     myprof.at(cls).at( ( var_tgt*nvar )+ivar).at( ( v_t*nvar )+j)->Fill( vali, valj, weight );
                  }
               }
            }
         }
      }
   }

   // correlation analysis for ranking  (single-target regression only)
   if (ntgt == 1) {
      for (UInt_t ivar=0; ivar<=nvar; ivar++) {
         xregmean[ivar] /= nevts;
         x2regmean[ivar] = x2regmean[ivar]/nevts - xregmean[ivar]*xregmean[ivar];
      }
      for (UInt_t ivar=0; ivar<nvar; ivar++) {
         xCregmean[ivar] = xCregmean[ivar]/nevts - xregmean[ivar]*xregmean[nvar];
         xCregmean[ivar] /= TMath::Sqrt( x2regmean[ivar]*x2regmean[nvar] );
      }

      fRanking.push_back( new Ranking( GetName() + "Transformation", "|Correlation with target|" ) );
      for (UInt_t ivar=0; ivar<nvar; ivar++) {
         Double_t abscor = TMath::Abs( xCregmean[ivar] );
         fRanking.back()->AddRank( Rank( fDataSetInfo.GetVariableInfo(ivar).GetLabel(), abscor ) );
      }

      if (nvar+ntgt <= (UInt_t)gConfig().GetVariablePlotting().fMaxNumOfAllowedVariablesForScatterPlots) {

         // compute also mutual information (non-linear correlation measure)
         fRanking.push_back( new Ranking( GetName() + "Transformation", "Mutual information" ) );
         for (UInt_t ivar=0; ivar<nvar; ivar++) {
            TH2F* h1 = mycorr.at(0).at( nvar ).at( ivar );
            Double_t mi = gTools().GetMutualInformation( *h1 );
            fRanking.back()->AddRank( Rank( fDataSetInfo.GetVariableInfo(ivar).GetLabel(), mi ) );
         }

         // compute correlation ratio (functional correlations measure)
         fRanking.push_back( new Ranking( GetName() + "Transformation", "Correlation Ratio" ) );
         for (UInt_t ivar=0; ivar<nvar; ivar++) {
            TH2F*    h2 = mycorr.at(0).at( nvar ).at( ivar );
            Double_t cr = gTools().GetCorrelationRatio( *h2 );
            fRanking.back()->AddRank( Rank( fDataSetInfo.GetVariableInfo(ivar).GetLabel(), cr ) );
         }

         // additionally compute correlation ratio from transposed histograms since correlation ratio is asymmetric
         fRanking.push_back( new Ranking( GetName() + "Transformation", "Correlation Ratio (T)" ) );
         for (UInt_t ivar=0; ivar<nvar; ivar++) {
            TH2F*    h2T = gTools().TransposeHist( *mycorr.at(0).at( nvar ).at( ivar ) );
            Double_t cr  = gTools().GetCorrelationRatio( *h2T  );
            fRanking.back()->AddRank( Rank( fDataSetInfo.GetVariableInfo(ivar).GetLabel(), cr ) );
            delete h2T;
         }
      }
   }
   // computes ranking of input variables
   // separation for 2-class classification
   else if (fDataSetInfo.GetNClasses() == 2
            && fDataSetInfo.GetClassInfo("Signal") != NULL
            && fDataSetInfo.GetClassInfo("Background") != NULL
            ) { // TODO: ugly hack.. adapt to new framework
      fRanking.push_back( new Ranking( GetName() + "Transformation", "Separation" ) );
      for (UInt_t i=0; i<nvar; i++) {
         Double_t sep = gTools().GetSeparation( hVars.at(fDataSetInfo.GetClassInfo("Signal")    ->GetNumber()).at(i),
                                                hVars.at(fDataSetInfo.GetClassInfo("Background")->GetNumber()).at(i) );
         fRanking.back()->AddRank( Rank( hVars.at(fDataSetInfo.GetClassInfo("Signal")->GetNumber()).at(i)->GetTitle(),
                                         sep ) );
      }
   }

   // for regression compute performance from correlation with target value

   // write histograms

   TDirectory* localDir = theDirectory;
   if (theDirectory == 0) {
      // create directory in root dir
      fRootBaseDir->cd();
      TString outputDir = TString("InputVariables");
      TListIter trIt(&fTransformations);
      while (VariableTransformBase *trf = (VariableTransformBase*) trIt())
         outputDir += "_" + TString(trf->GetShortName());

      TString uniqueOutputDir = outputDir;
      Int_t counter = 0;
      TObject* o = NULL;
      while( (o = fRootBaseDir->FindObject(uniqueOutputDir)) != 0 ){
         uniqueOutputDir = outputDir+Form("_%d",counter);
         Log() << kINFO << "A " << o->ClassName() << " with name " << o->GetName() << " already exists in "
               << fRootBaseDir->GetPath() << ", I will try with "<<uniqueOutputDir<<"." << Endl;
         ++counter;
      }

      //       TObject* o = fRootBaseDir->FindObject(outputDir);
      //       if (o != 0) {
      //          Log() << kFATAL << "A " << o->ClassName() << " with name " << o->GetName() << " already exists in "
      //                << fRootBaseDir->GetPath() << "("<<outputDir<<")" << Endl;
      //       }
      localDir = fRootBaseDir->mkdir( uniqueOutputDir );
      localDir->cd();

      Log() << kVERBOSE << "Create and switch to directory " << localDir->GetPath() << Endl;
   }
   else {
      theDirectory->cd();
   }

   for (UInt_t i=0; i<nvar+ntgt; i++) {
      for (Int_t cls = 0; cls < ncls; cls++) {
         if (hVars.at(cls).at(i) != 0) {
            hVars.at(cls).at(i)->Write();
            hVars.at(cls).at(i)->SetDirectory(0);
            delete hVars.at(cls).at(i);
         }
      }
   }

   // correlation plots have dedicated directory
   if (nvar+ntgt <= (UInt_t)gConfig().GetVariablePlotting().fMaxNumOfAllowedVariablesForScatterPlots) {

      localDir = localDir->mkdir( "CorrelationPlots" );
      localDir ->cd();
      Log() << kDEBUG << "Create scatter and profile plots in target-file directory: " << Endl;
      Log() << kDEBUG << localDir->GetPath() << Endl;


      for (UInt_t i=0; i<nvar+ntgt; i++) {
         for (UInt_t j=i+1; j<nvar+ntgt; j++) {
            for (Int_t cls = 0; cls < ncls; cls++) {
               if (mycorr.at(cls).at(i).at(j) != 0 ) {
                  mycorr.at(cls).at(i).at(j)->Write();
                  mycorr.at(cls).at(i).at(j)->SetDirectory(0);
                  delete mycorr.at(cls).at(i).at(j);
               }
               if (myprof.at(cls).at(i).at(j) != 0) {
                  myprof.at(cls).at(i).at(j)->Write();
                  myprof.at(cls).at(i).at(j)->SetDirectory(0);
                  delete myprof.at(cls).at(i).at(j);
               }
            }
         }
      }
   }
   if (theDirectory != 0 ) theDirectory->cd();
   else                    fRootBaseDir->cd();
}

////////////////////////////////////////////////////////////////////////////////
/// returns string for transformation

std::vector<TString>* TMVA::TransformationHandler::GetTransformationStringsOfLastTransform() const
{
   VariableTransformBase* trf = ((VariableTransformBase*)GetTransformationList().Last());
   if (!trf) return 0;
   else      return trf->GetTransformationStrings( fTransformationsReferenceClasses.back() );
}

////////////////////////////////////////////////////////////////////////////////
/// returns string for transformation

const char* TMVA::TransformationHandler::GetNameOfLastTransform() const
{
   VariableTransformBase* trf = ((VariableTransformBase*)GetTransformationList().Last());
   if (!trf) return 0;
   else      return trf->GetName();
}

////////////////////////////////////////////////////////////////////////////////
/// write transformation to stream

void TMVA::TransformationHandler::WriteToStream( std::ostream& o ) const
{
   TListIter trIt(&fTransformations);
   std::vector< Int_t >::const_iterator rClsIt = fTransformationsReferenceClasses.begin();

   o << "NTransformtations " << fTransformations.GetSize() << std::endl << std::endl;

   ClassInfo* ci;
   UInt_t i = 1;
   while (VariableTransformBase *trf = (VariableTransformBase*) trIt()) {
      o << "#TR -*-*-*-*-*-*-* transformation " << i++ << ": " << trf->GetName() << " -*-*-*-*-*-*-*-" << std::endl;
      trf->WriteTransformationToStream(o);
      ci = fDataSetInfo.GetClassInfo( (*rClsIt) );
      TString clsName;
      if (ci == 0 ) clsName = "AllClasses";
      else clsName = ci->GetName();
      o << "ReferenceClass " << clsName << std::endl;
      ++rClsIt;
   }
}


////////////////////////////////////////////////////////////////////////////////
/// XML node describing the transformation

void TMVA::TransformationHandler::AddXMLTo( void* parent ) const
{
   if(!parent) return;
   void* trfs = gTools().AddChild(parent, "Transformations");
   gTools().AddAttr( trfs, "NTransformations", fTransformations.GetSize() );
   TListIter trIt(&fTransformations);
   while (VariableTransformBase *trf = (VariableTransformBase*) trIt()) trf->AttachXMLTo(trfs);
}

////////////////////////////////////////////////////////////////////////////////

void TMVA::TransformationHandler::ReadFromStream( std::istream& )
{
   Log() << kFATAL << "Read transformations not implemented" << Endl;
   // TODO
}

////////////////////////////////////////////////////////////////////////////////

void TMVA::TransformationHandler::ReadFromXML( void* trfsnode )
{
   void* ch = gTools().GetChild( trfsnode );
   while(ch) {
      Int_t idxCls = -1;
      TString trfname;
      gTools().ReadAttr(ch, "Name", trfname);

      VariableTransformBase* newtrf = 0;

      if (trfname == "Decorrelation" ) {
         newtrf = new VariableDecorrTransform(fDataSetInfo);
      }
      else if (trfname == "PCA" ) {
         newtrf = new VariablePCATransform(fDataSetInfo);
      }
      else if (trfname == "Gauss" ) {
         newtrf = new VariableGaussTransform(fDataSetInfo);
      }
      else if (trfname == "Uniform" ) {
         newtrf = new VariableGaussTransform(fDataSetInfo, "Uniform");
      }
      else if (trfname == "Normalize" ) {
         newtrf = new VariableNormalizeTransform(fDataSetInfo);
      }
      else if (trfname == "Rearrange" ) {
         newtrf = new VariableRearrangeTransform(fDataSetInfo);
      }
      else if (trfname != "None") {
      }
      else {
         Log() << kFATAL << "<ReadFromXML> Variable transform '"
               << trfname << "' unknown." << Endl;
      }
      newtrf->ReadFromXML( ch );
      AddTransformation( newtrf, idxCls );
      ch = gTools().GetNextChild(ch);
   }
}

////////////////////////////////////////////////////////////////////////////////
/// prints ranking of input variables

void TMVA::TransformationHandler::PrintVariableRanking() const
{
  //Log() << kINFO << " " << Endl;
   Log() << kINFO << "Ranking input variables (method unspecific)..." << Endl;
   std::vector<Ranking*>::const_iterator it = fRanking.begin();
   for (; it != fRanking.end(); ++it) (*it)->Print();
}

////////////////////////////////////////////////////////////////////////////////

Double_t TMVA::TransformationHandler::GetMean( Int_t ivar, Int_t cls ) const
{
   try {
      return fVariableStats.at(cls).at(ivar).fMean;
   }
   catch(...) {
      try {
         return fVariableStats.at(fNumC-1).at(ivar).fMean;
      }
      catch(...) {
         Log() << kWARNING << "Inconsistent variable state when reading the mean value. " << Endl;
      }
   }
   Log() << kWARNING << "Inconsistent variable state when reading the mean value. Value 0 given back" << Endl;
   return 0;
}

////////////////////////////////////////////////////////////////////////////////

Double_t TMVA::TransformationHandler::GetRMS( Int_t ivar, Int_t cls ) const
{
   try {
      return fVariableStats.at(cls).at(ivar).fRMS;
   }
   catch(...) {
      try {
         return fVariableStats.at(fNumC-1).at(ivar).fRMS;
      }
      catch(...) {
         Log() << kWARNING << "Inconsistent variable state when reading the RMS value. " << Endl;
      }
   }
   Log() << kWARNING << "Inconsistent variable state when reading the RMS value. Value 0 given back" << Endl;
   return 0;
}

////////////////////////////////////////////////////////////////////////////////

Double_t TMVA::TransformationHandler::GetMin( Int_t ivar, Int_t cls ) const
{
   try {
      return fVariableStats.at(cls).at(ivar).fMin;
   }
   catch(...) {
      try {
         return fVariableStats.at(fNumC-1).at(ivar).fMin;
      }
      catch(...) {
         Log() << kWARNING << "Inconsistent variable state when reading the minimum value. " << Endl;
      }
   }
   Log() << kWARNING << "Inconsistent variable state when reading the minimum value. Value 0 given back" << Endl;
   return 0;
}

////////////////////////////////////////////////////////////////////////////////

Double_t TMVA::TransformationHandler::GetMax( Int_t ivar, Int_t cls ) const
{
   try {
      return fVariableStats.at(cls).at(ivar).fMax;
   }
   catch(...) {
      try {
         return fVariableStats.at(fNumC-1).at(ivar).fMax;
      }
      catch(...) {
         Log() << kWARNING << "Inconsistent variable state when reading the maximum value. " << Endl;
      }
   }
   Log() << kWARNING << "Inconsistent variable state when reading the maximum value. Value 0 given back" << Endl;
   return 0;
}