THnSparse feature request

  1. Would it be possible to add a THnSparse constructor to construct on a TH1* (and derivatives)?

  2. It would be useful to have {Multiply,Divide,Add} for TF1’s, and to be able to specify which dimension to apply it on to? (i.e. the x of TF1 should be applied to Axis(i))

Cheers,

  • Pete

Pete,

We prefer not to connect very basic classes like TH1 to more elaborated classes like THnSparse. This prevents possible future reorganization of classes in different libraries.
You can achieve the same result with THnSparse::Projection.

Rene

Rene, I was referring to doing

TH2F *MyHist…
THnSparseF newHist(MyHist);

Rather than newHist -> MyHist.

  • Pete

OK, sorry if I misread your original mail.
Would you be willing to provide the code? ::slight_smile:

Rene

Sure. Here is something that does 2) for Multiply, although I wasn’t sure exactly what to do for the SetError() because I couldn’t access the sumw2, since that is buried behind private data members…

It should be trivial to modify for Add/Divide? I wanted to post this to make sure I had things along the right lines…

void THnSparse__Multiply( THnSparse *h, TF1 *f )
{
    Long64_t TotalEntries = h->GetNbins();
    UInt_t nDims = h->GetNdimensions();
    
    Int_t   *index = new Int_t[ nDims ];
    Double_t *coord = new Double_t[ nDims ];
    memset( coord, 0, sizeof( Float_t ) * nDims );
    
    Float_t content = 0, error = 0;
    
    for ( Long64_t i = 0; i < TotalEntries; i++ )
    {
        content = h->GetBinContent( i, index );
        error   = h->GetBinError( i );
        
        for ( UInt_t j = 0; j < nDims; j++ )
            coord[j] = h->GetAxis( j )->GetBinCenter( index[j] );
        
        Float_t cu = f->EvalPar( coord, NULL );
        h->SetBinContent( index, cu * content );
        h->SetBinError  ( index, cu * error );
    }
    
    delete [] coord;
    delete [] index;
}

void THnSparse__MultiplyDim( THnSparse *h, TF1 *f, UInt_t Dim = 0 )
{
    Long64_t TotalEntries = h->GetNbins();
    UInt_t nDims = h->GetNdimensions();
    
    Int_t   *index = new Int_t[ nDims ];
    Double_t *coord = new Double_t[ nDims ];
    memset( coord, 0, sizeof( Float_t ) * nDims );
    
    UInt_t nBins = h->GetAxis( Dim )->GetNbins();
    TAxis *ax = h->GetAxis(Dim);
    Double_t *CachedValues = new Double_t[ nBins ];
    for ( UInt_t i = 0; i < nBins; i++ )
        CachedValues[i] = f->Eval( ax->GetBinCenter( i ) );
    
    Double_t content = 0, error = 0;
    
    for ( Long64_t i = 0; i < TotalEntries; i++ )
    {
        content = h->GetBinContent( i, index );
        error   = h->GetBinError( i );
        Double_t cu = CachedValues[ index[ Dim ] ];
        h->SetBinContent( index, cu * content );
        h->SetBinError  ( index, cu * error );
    }
    
    delete [] CachedValues;
    delete [] coord;
    delete [] index;
}

Pete,

Sorry to ask too much ::slight_smile:
Could you document these functions and fix problems like

memset( coord, 0, sizeof( Float_t ) * nDims ); ==> memset( coord, 0, sizeof( Double_t ) * nDims );

and use Double_t internally for intermediate vars?

Rene

Is this looking closer?

[code]#include <TAxis.h>
#include <TF1.h>

void THnSparse::Multiply( TF1 f1 )
{
// Performs the operation: this = this
c1*f1
// if errors are defined, errors are also recalculated.
//
// Only bins inside the function range are recomputed.
// IMPORTANT NOTE: If you intend to use the errors of this histogram later
// you should call Sumw2 before making this operation.
// This is particularly important if you fit the histogram after THnSparse::Multiply

UInt_t nDims = GetNdimensions();

Int_t    *index = new Int_t[ nDims ];
Double_t *coord = new Double_t[ nDims ];
memset( coord, 0, sizeof( Double_t ) * nDims );

Double_t content = 0, error = 0;

ULong64_t TotalEntries = GetNbins();
for ( ULong64_t i = 0; i < TotalEntries; i++ )
{
    content = GetBinContent( i, index );
    #error Don't know how to determine whether sumw2 is on
    if ( Sum of weights is on? )
        error   = GetBinError( i );
    
    // Get the bin co-ordinates given an index
    for ( UInt_t j = 0; j < nDims; j++ )
        coord[j] = GetAxis( j )->GetBinCenter( index[j] );
    
    if ( !f1->IsInside( coord ) ) 
        continue;
    TF1::RejectPoint(kFALSE);

    // Evaulate function at coord
    Double_t cu = f1->EvalPar( coord, NULL );
    
    SetBinContent( index, cu * content );
    // *** FIXME: Not sure what is correct to do here for Sumw2 and
    //            related variables
    #error Don't know how to determine whether sumw2 is on
    if ( Sum of weights is on? )
        SetBinError  ( index, cu * error );
}

delete [] coord;
delete [] index;

}

void THnSparse::MultiplyDim( TF1 f1, UInt_t Dim = 0, Bool_t UseCache = false )
{
// Performs the operation: this = this
c1*f1
// if errors are defined, errors are also recalculated.
//
// Only bins inside the function range are recomputed.
// IMPORTANT NOTE: If you intend to use the errors of this histogram later
// you should call Sumw2 before making this operation.
// This is particularly important if you fit the histogram after THnSparse::Multiply

UInt_t nDims = GetNdimensions();
if ( Dim >= nDims )
{
    // Raise an error here
    #error Need some error handling code
    return;
}

Int_t *index = new Int_t[ nDims ];
Double_t *CachedValues = NULL;

TAxis *multiplyAxis = GetAxis( Dim );

if ( UseCache )
{
    // Cache all of the values along our axis, if the function is slow
    UInt_t nBins = GetAxis( Dim )->GetNbins();
    CachedValues = new Double_t[ nBins ];
    for ( UInt_t i = 0; i < nBins; i++ )
        CachedValues[i] = f1->Eval( multiplyAxis->GetBinCenter( i ) );
}

Double_t content = 0, error = 0, cu = 0;

Long64_t TotalEntries = GetNbins();
for ( Long64_t i = 0; i < TotalEntries; i++ )
{
    content = GetBinContent( i, index );
    #error Don't know how to determine whether sumw2 is on
    if ( Sum of weights is on? )
        error = GetBinError( i );
    
    // FIXME: 
    // Maybe refactor this outside of the loop (Have two seperate loops?)
    if ( UseCache )
        cu = CachedValues[ index[ Dim ] ];
    else
        cu = f1->Eval( multiplyAxis->GetBinCenter( i ) );
    
    SetBinContent( index, cu * content );
    // *** FIXME: Not sure what is correct to do here for Sumw2 and
    //            related variables
    #error Don't know how to determine whether sumw2 is on
    if ( Sum of weights is on? )
        SetBinError  ( index, cu * error );
}

if ( CachedValues ) delete [] CachedValues;
delete [] index;

}
[/code]

Here is some sample code for 1).

Untested - sorry.

[code]#include <TH1.h>

THnSparse::THnSparse( TH1 *originalHist, Int_t chunksize ):
TNamed(name, title), fNdimensions( originalHist->GetDimension() ),
fChunkSize(chunksize), fFilledBins(0),
fAxes( originalHist->GetDimension() ), fEntries(0), fTsumw(0), fTsumw2(-1.),
fTsumwx( originalHist->GetDimension() ),
fTsumwx2( originalHist->GetDimension() ),
fCompactCoord(0), fIntegral(0), fIntegralStatus(kNoInt)
{
// Construct a THnSparse from a dense histogram
Int_t nbins[3] = { -1, -1, -1 };
TAxis *axes[] = { originalHist->GetXaxis(),
originalHist->GetYaxis(),
originalHist->GetZaxis() };

for ( Int_t i = 0; i < fNdimensions; ++i )
{
TAxis *axis = axes[i]->Clone();
fAxes.AddAtAndExpand( axis, i );
nbins[i] = axis->GetNbins();
}
fAxes.SetOwner();

fCompactCoord = new THnSparseCompactBinCoord( dim, nbins );
fBinContent.SetOwner();

Double_t value = 0;

Int_t coord[3] = { 0, 0, 0 };
UInt_t origBin = 0;
for ( coord[2] = 0; coord[2] <= nbins[2] + 1; coord[2]++ )
for ( coord[1] = 0; coord[1] <= nbins[1] + 1; coord[1]++ )
for ( coord[0] = 0; coord[0] <= nbins[0] + 1; coord[0]++ )
{
origBin = originalHist->FindBin( coord[0], coord[1], coord[2] );
value = originalHist->GetBinContent( origBin );
if ( value )
{
SetBinContent( coord, value );
#error Don’t know how to determine whether sumw2 is on
if ( sumw2 is on )
{
GetCompactCoord()->SetCoord( coord );
Long_t bin = GetBinIndexForCurrentBin(kTRUE);

        THnSparseArrayChunk* chunk = GetChunk( bin / fChunkSize );
        chunk->fSumw2->SetAt( originalHist->GetSumw2()[origBin], bin % fChunkSize);
     }
  }

}

fEntries = originalHist->GetEntries();

// Not sure what to do for these (THnSparse stores sumwx{,2} for each dim)
//~ fTsumw = originalHist->fTsumw;
//~ fTsumw2 = originalHist->fTsumw2;
//~ fTsumwx = originalHist->fTsumwx;
//~ fTsumwx2 = originalHist->fTsumwx2;

delete [] nbins;
}[/code]

Thanks Pete for the code. We will see how to adapt it to THnSparse for the next release in December. Lorenzo will get in touch with you if necessary ::slight_smile:

Rene

Thanks Rene :slight_smile:

  • Pete

Thank you for the code.

We are in the process of improving the Q/A for the next December release, and adding a test suite for histograms including the THnSparse.
We will test also this new code and when it will be validated we will include it in the root trunk.
I’ll let you know when the code will be committed or before if needed

Best Regards

Lorenzo

Another question that comes to mind: Are we far off being able to TTree::Draw to THnSparse?

  • Pete

This will not come before Christmas. We also need better multi-D visualization techniques. Filling a THnSparse should not be too difficult as the logic in TTree::Draw is already in place to process queries like
T.Draw(“x1:x2:x3:…:angry:n”) where n can be anything.

Rene

Did anything come of this?

Hi,

These two features are now implemented in the math branch. They will probably be out in the next ROOT release.

David

Out of curiosity, what of TTree::Draw >> THnSparse ?

1 Like