Histogram operation optimization

Hi,

I have implemented some histogram operations (flip x axis, fourier transform, convolution,…)
This is done bin by bin and I would like to speed up the process. How would you proceed ?

I was thinking about getting the bin content array and loop over vectors with some -O2 compilation optimizations, but this is a bit more complex with 2D histograms…

Sorry for the silly question and confusions :slight_smile:

This is not a silly question. Your message suggests you managed to do what you want for 1D histograms.
Maybe you can post here what you did for 1D. In any case, the best person to help you is @moneta .

Hi

Sure, here is a typical naive example that I would like to optimize:

TH1* Flip(TH1* h, Option_t* option)
{
        if( h == NULL) return NULL;

        TString opt = option;
        opt.ToLower();

        TString hname = (TString) h->GetName() + ":" + opt + "-flip";
        TH1 *h0 = (TH1*) h->Clone();

        TAxis *x1, *x2, *x3;
        if(opt == "x") {

                x3 = h->GetXaxis();
                x1 = h->GetYaxis();
                x2 = h->GetZaxis();

        } else if(opt == "y") {

                x3 = h->GetYaxis();
                x1 = h->GetXaxis();
                x2 = h->GetZaxis();

        } else if(opt == "z") {

                x3 = h->GetZaxis();
                x1 = h->GetXaxis();
                x2 = h->GetYaxis();

        } else {

                
                return NULL;
        }

        int N = h->GetEntries();
        for(int x1_i = 1, x1_N = x1->GetNbins(); x1_i <= x1_N; x1_i++) {

                for(int x2_i = 1, x2_N = x2->GetNbins(); x2_i <= x2_N; x2_i++) {

                        double c1 = 0, c2 = 0;
                        double c1_err = 0, c2_err = 0;

                        for(int x3_i = 1, x3_N = x3->GetNbins(); x3_i <= x3_N; x3_i++) {

                                if(opt == "x") {

                                        c1     = h->GetBinContent(x3_i, x1_i, x2_i);
                                        c1_err = h->GetBinError(x3_i, x1_i, x2_i);
                                        c2     = h->GetBinContent(x3->GetNbins()+1-x3_i, x1_i, x2_i);
                                        c2_err = h->GetBinError(x3->GetNbins()+1-x3_i, x1_i, x2_i);

                                        h0->SetBinContent(x3_i, x1_i, x2_i, c2);
                                        h0->SetBinError(x3_i, x1_i, x2_i, c2_err);
                                        h0->SetBinContent(x3->GetNbins()+1-x3_i, x1_i, x2_i, c1);
                                        h0->SetBinError(x3->GetNbins()+1-x3_i, x1_i, x2_i, c1_err);
                                        
                                } else if (opt == "y") {

                                        c1     = h->GetBinContent(x1_i, x3_i, x2_i);
                                        c1_err = h->GetBinError(x1_i, x3_i, x2_i);
                                        c2     = h->GetBinContent(x1_i, x3->GetNbins()+1-x3_i, x2_i);
                                        c2_err = h->GetBinError(x1_i, x3->GetNbins()+1-x3_i, x2_i);

                                        h0->SetBinContent(x1_i, x3_i, x2_i, c2);
                                        h0->SetBinError(x1_i, x3_i, x2_i, c2_err);
                                        h0->SetBinContent(x1_i, x3->GetNbins()+1-x3_i, x2_i, c1);
                                        h0->SetBinError(x1_i, x3->GetNbins()+1-x3_i, x2_i, c1_err);
                                        
                                } else if (opt == "z") {

                                        c1     = h->GetBinContent(x1_i, x2_i, x3_i);
                                        c1_err = h->GetBinError(x1_i, x2_i, x3_i);
                                        c2     = h->GetBinContent(x1_i, x2_i, x3->GetNbins()+1-x3_i);
                                        c2_err = h->GetBinError(x1_i, x2_i, x3->GetNbins()+1-x3_i);

                                        h0->SetBinContent(x1_i, x2_i, x3_i, c2);
                                        h0->SetBinError(x1_i, x2_i, x3_i, c2_err);
                                        h0->SetBinContent(x1_i, x2_i, x3->GetNbins()+1-x3_i, c1);
                                        h0->SetBinError(x1_i, x2_i, x3->GetNbins()+1-x3_i, c1_err);
                                }
                        }
                }
        }

        h0->SetEntries(N);
        
        return h0;
}

So it seems you implement your code for all cases ? X Y and Z ? what is the issue with the example ?

This code is handling one by one the histogram bin flip.
So I was wondering whether it is possible le to take advantage of compiler optimizations, vectorization or parallelization of loops if possible.

I am not sure if there is a better way or not. I let @moneta reply.

Hi,
At the moment we don’t have implemented parallelisation of loops of histogram contents. For taking advantage of vectorisation you should probably return the content in a the histogram array.
Otherwise for convenience, and faster than looping on bins, is using the class THistRange and TBinIterator, see ROOT: THistRange Class Reference

Lorenzo

1 Like