ROOT logo
// @(#)root/roostats:$Id: MarkovChain.cxx 31276 2009-11-18 15:06:42Z moneta $
// Authors: Kevin Belasco        17/06/2009
// Authors: Kyle Cranmer         17/06/2009
/*************************************************************************
 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

//_________________________________________________
/*
BEGIN_HTML
<p>
Stores the steps in a Markov Chain of points.  Allows user to access the
weight and NLL value (if applicable) with which a point was added to the
MarkovChain.
</p>
END_HTML
*/
//_________________________________________________

#ifndef ROOT_Rtypes
#include "Rtypes.h"
#endif

#ifndef ROOT_TNamed
#include "TNamed.h"
#endif
#ifndef ROOSTATS_MarkovChain
#include "RooStats/MarkovChain.h"
#endif
#ifndef ROO_DATA_SET
#include "RooDataSet.h"
#endif
#ifndef ROO_ARG_SET
#include "RooArgSet.h"
#endif
#ifndef ROO_REAL_VAR
#include "RooRealVar.h"
#endif
#ifndef RooStats_RooStatsUtils
#include "RooStats/RooStatsUtils.h"
#endif
#ifndef ROO_DATA_HIST
#include "RooDataHist.h"
#endif
#ifndef ROOT_THnSparse
#include "THnSparse.h"
#endif

ClassImp(RooStats::MarkovChain);

using namespace RooFit;
using namespace RooStats;

static const char* WEIGHT_NAME = "weight_MarkovChain_local_";
static const char* NLL_NAME = "nll_MarkovChain_local_";
static const char* DATASET_NAME = "dataset_MarkovChain_local_";
static const char* DEFAULT_NAME = "_markov_chain";
static const char* DEFAULT_TITLE = "Markov Chain";

MarkovChain::MarkovChain() :
   TNamed(DEFAULT_NAME, DEFAULT_TITLE)
{
   fParameters = NULL;
   fDataEntry = NULL;
   fChain = NULL;
   fNLL = NULL;
   fWeight = NULL;
}

MarkovChain::MarkovChain(RooArgSet& parameters) :
   TNamed(DEFAULT_NAME, DEFAULT_TITLE)
{
   fParameters = NULL;
   fDataEntry = NULL;
   fChain = NULL;
   fNLL = NULL;
   fWeight = NULL;
   SetParameters(parameters);
}

MarkovChain::MarkovChain(const char* name, const char* title,
      RooArgSet& parameters) : TNamed(name, title)
{
   fParameters = NULL;
   fDataEntry = NULL;
   fChain = NULL;
   fNLL = NULL;
   fWeight = NULL;
   SetParameters(parameters);
}

void MarkovChain::SetParameters(RooArgSet& parameters)
{
   delete fChain;
   delete fParameters;
   delete fDataEntry;

   fParameters = new RooArgSet();
   fParameters->addClone(parameters);

   // kbelasco: consider setting fDataEntry = fChain->get()
   // to see if that makes it possible to get values of variables without
   // doing string comparison
   RooRealVar nll(NLL_NAME, "-log Likelihood", 0);
   RooRealVar weight(WEIGHT_NAME, "weight", 0);

   fDataEntry = new RooArgSet();
   fDataEntry->addClone(parameters);
   fDataEntry->addClone(nll);
   fDataEntry->addClone(weight);
   fNLL = (RooRealVar*)fDataEntry->find(NLL_NAME);
   fWeight = (RooRealVar*)fDataEntry->find(WEIGHT_NAME);

   fChain = new RooDataSet(DATASET_NAME, "Markov Chain", *fDataEntry,WEIGHT_NAME);
}

void MarkovChain::Add(RooArgSet& entry, Double_t nllValue, Double_t weight)
{
   if (fParameters == NULL)
      SetParameters(entry);
   RooStats::SetParameters(&entry, fDataEntry);
   fNLL->setVal(nllValue);
   //kbelasco: this is stupid, but some things might require it, so be doubly sure
   fWeight->setVal(weight);
   fChain->add(*fDataEntry, weight);
   //fChain->add(*fDataEntry);
}

void MarkovChain::AddFast(RooArgSet& entry, Double_t nllValue, Double_t weight)
{
   RooStats::SetParameters(&entry, fDataEntry);
   fNLL->setVal(nllValue);
   //kbelasco: this is stupid, but some things might require it, so be doubly sure
   fWeight->setVal(weight);
   fChain->addFast(*fDataEntry, weight);
   //fChain->addFast(*fDataEntry);
}

RooDataSet* MarkovChain::GetAsDataSet(RooArgSet* whichVars) const
{
   RooArgSet args;
   if (whichVars == NULL) {
      //args.add(*fParameters);
      //args.add(*fNLL);
      args.add(*fDataEntry);
   } else {
      args.add(*whichVars);
   }

   RooDataSet* data;
   //data = dynamic_cast<RooDataSet*>(fChain->reduce(args));
   data = (RooDataSet*)fChain->reduce(args);

   return data;
}

RooDataSet* MarkovChain::GetAsDataSet(RooCmdArg arg1, RooCmdArg arg2,
      RooCmdArg arg3, RooCmdArg arg4, RooCmdArg arg5, RooCmdArg arg6,
      RooCmdArg arg7, RooCmdArg arg8) const
{
   RooDataSet* data;
   data = (RooDataSet*)fChain->reduce(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8);
   return data;
}

RooDataHist* MarkovChain::GetAsDataHist(RooArgSet* whichVars) const
{
   RooArgSet args;
   if (whichVars == NULL) {
      args.add(*fParameters);
      //args.add(*fNLL);
      //args.add(*fDataEntry);
   } else {
      args.add(*whichVars);
   }

   RooDataSet* data = (RooDataSet*)fChain->reduce(args);
   RooDataHist* hist = data->binnedClone();
   delete data;

   return hist;
}

RooDataHist* MarkovChain::GetAsDataHist(RooCmdArg arg1, RooCmdArg arg2,
      RooCmdArg arg3, RooCmdArg arg4, RooCmdArg arg5, RooCmdArg arg6,
      RooCmdArg arg7, RooCmdArg arg8) const
{
   RooDataSet* data;
   RooDataHist* hist;
   data = (RooDataSet*)fChain->reduce(arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8);
   hist = data->binnedClone();
   delete data;

   return hist;
}

THnSparse* MarkovChain::GetAsSparseHist(RooAbsCollection* whichVars) const
{
   RooArgList axes;
   if (whichVars == NULL)
      axes.add(*fParameters);
   else 
      axes.add(*whichVars);

   Int_t dim = axes.getSize();
   Double_t* min = new Double_t[dim];
   Double_t* max = new Double_t[dim];
   Int_t* bins = new Int_t[dim];
   TIterator* it = axes.createIterator();
   for (Int_t i = 0; i < dim; i++) {
      min[i] = ((RooRealVar*)it->Next())->getMin();
      max[i] = ((RooRealVar*)it->Next())->getMax();
      bins[i] = ((RooRealVar*)it->Next())->numBins();
   }
   THnSparseF* sparseHist = new THnSparseF("posterior", "MCMC Posterior Histogram",
         dim, bins, min, max);
   // kbelasco: check here for memory leaks: does THnSparse copy min, max, and bins
   // or can we not delete them?
   //delete[] min;
   //delete[] max;
   //delete[] bins;

   // kbelasco: it appears we need to call Sumw2() just to get the
   // histogram to keep a running total of the weight so that Getsumw doesn't
   // just return 0
   sparseHist->Sumw2();

   // Fill histogram
   Int_t size = fChain->numEntries();
   const RooArgSet* entry;
   Double_t* x = new Double_t[dim];
   for (Int_t i = 0; i < size; i++) {
      entry = fChain->get(i);
      it->Reset();
      for (Int_t ii = 0; ii < dim; ii++)
         x[ii] = entry->getRealValue(it->Next()->GetName());
      sparseHist->Fill(x, fChain->weight());
   }
   delete[] x;
   delete it;

   return sparseHist;
}

Double_t MarkovChain::NLL(Int_t i) const
{
   // kbelasco: how to do this?
   //fChain->get(i);
   //return fNLL->getVal();
   return fChain->get(i)->getRealValue(NLL_NAME);
}

Double_t MarkovChain::NLL() const
{
   // kbelasco: how to do this?
   //fChain->get();
   //return fNLL->getVal();
   return fChain->get()->getRealValue(NLL_NAME);
}

Double_t MarkovChain::Weight() const
{
   return fChain->weight();
}

Double_t MarkovChain::Weight(Int_t i) const
{
   fChain->get(i);
   return fChain->weight();
}
 MarkovChain.cxx:1
 MarkovChain.cxx:2
 MarkovChain.cxx:3
 MarkovChain.cxx:4
 MarkovChain.cxx:5
 MarkovChain.cxx:6
 MarkovChain.cxx:7
 MarkovChain.cxx:8
 MarkovChain.cxx:9
 MarkovChain.cxx:10
 MarkovChain.cxx:11
 MarkovChain.cxx:12
 MarkovChain.cxx:13
 MarkovChain.cxx:14
 MarkovChain.cxx:15
 MarkovChain.cxx:16
 MarkovChain.cxx:17
 MarkovChain.cxx:18
 MarkovChain.cxx:19
 MarkovChain.cxx:20
 MarkovChain.cxx:21
 MarkovChain.cxx:22
 MarkovChain.cxx:23
 MarkovChain.cxx:24
 MarkovChain.cxx:25
 MarkovChain.cxx:26
 MarkovChain.cxx:27
 MarkovChain.cxx:28
 MarkovChain.cxx:29
 MarkovChain.cxx:30
 MarkovChain.cxx:31
 MarkovChain.cxx:32
 MarkovChain.cxx:33
 MarkovChain.cxx:34
 MarkovChain.cxx:35
 MarkovChain.cxx:36
 MarkovChain.cxx:37
 MarkovChain.cxx:38
 MarkovChain.cxx:39
 MarkovChain.cxx:40
 MarkovChain.cxx:41
 MarkovChain.cxx:42
 MarkovChain.cxx:43
 MarkovChain.cxx:44
 MarkovChain.cxx:45
 MarkovChain.cxx:46
 MarkovChain.cxx:47
 MarkovChain.cxx:48
 MarkovChain.cxx:49
 MarkovChain.cxx:50
 MarkovChain.cxx:51
 MarkovChain.cxx:52
 MarkovChain.cxx:53
 MarkovChain.cxx:54
 MarkovChain.cxx:55
 MarkovChain.cxx:56
 MarkovChain.cxx:57
 MarkovChain.cxx:58
 MarkovChain.cxx:59
 MarkovChain.cxx:60
 MarkovChain.cxx:61
 MarkovChain.cxx:62
 MarkovChain.cxx:63
 MarkovChain.cxx:64
 MarkovChain.cxx:65
 MarkovChain.cxx:66
 MarkovChain.cxx:67
 MarkovChain.cxx:68
 MarkovChain.cxx:69
 MarkovChain.cxx:70
 MarkovChain.cxx:71
 MarkovChain.cxx:72
 MarkovChain.cxx:73
 MarkovChain.cxx:74
 MarkovChain.cxx:75
 MarkovChain.cxx:76
 MarkovChain.cxx:77
 MarkovChain.cxx:78
 MarkovChain.cxx:79
 MarkovChain.cxx:80
 MarkovChain.cxx:81
 MarkovChain.cxx:82
 MarkovChain.cxx:83
 MarkovChain.cxx:84
 MarkovChain.cxx:85
 MarkovChain.cxx:86
 MarkovChain.cxx:87
 MarkovChain.cxx:88
 MarkovChain.cxx:89
 MarkovChain.cxx:90
 MarkovChain.cxx:91
 MarkovChain.cxx:92
 MarkovChain.cxx:93
 MarkovChain.cxx:94
 MarkovChain.cxx:95
 MarkovChain.cxx:96
 MarkovChain.cxx:97
 MarkovChain.cxx:98
 MarkovChain.cxx:99
 MarkovChain.cxx:100
 MarkovChain.cxx:101
 MarkovChain.cxx:102
 MarkovChain.cxx:103
 MarkovChain.cxx:104
 MarkovChain.cxx:105
 MarkovChain.cxx:106
 MarkovChain.cxx:107
 MarkovChain.cxx:108
 MarkovChain.cxx:109
 MarkovChain.cxx:110
 MarkovChain.cxx:111
 MarkovChain.cxx:112
 MarkovChain.cxx:113
 MarkovChain.cxx:114
 MarkovChain.cxx:115
 MarkovChain.cxx:116
 MarkovChain.cxx:117
 MarkovChain.cxx:118
 MarkovChain.cxx:119
 MarkovChain.cxx:120
 MarkovChain.cxx:121
 MarkovChain.cxx:122
 MarkovChain.cxx:123
 MarkovChain.cxx:124
 MarkovChain.cxx:125
 MarkovChain.cxx:126
 MarkovChain.cxx:127
 MarkovChain.cxx:128
 MarkovChain.cxx:129
 MarkovChain.cxx:130
 MarkovChain.cxx:131
 MarkovChain.cxx:132
 MarkovChain.cxx:133
 MarkovChain.cxx:134
 MarkovChain.cxx:135
 MarkovChain.cxx:136
 MarkovChain.cxx:137
 MarkovChain.cxx:138
 MarkovChain.cxx:139
 MarkovChain.cxx:140
 MarkovChain.cxx:141
 MarkovChain.cxx:142
 MarkovChain.cxx:143
 MarkovChain.cxx:144
 MarkovChain.cxx:145
 MarkovChain.cxx:146
 MarkovChain.cxx:147
 MarkovChain.cxx:148
 MarkovChain.cxx:149
 MarkovChain.cxx:150
 MarkovChain.cxx:151
 MarkovChain.cxx:152
 MarkovChain.cxx:153
 MarkovChain.cxx:154
 MarkovChain.cxx:155
 MarkovChain.cxx:156
 MarkovChain.cxx:157
 MarkovChain.cxx:158
 MarkovChain.cxx:159
 MarkovChain.cxx:160
 MarkovChain.cxx:161
 MarkovChain.cxx:162
 MarkovChain.cxx:163
 MarkovChain.cxx:164
 MarkovChain.cxx:165
 MarkovChain.cxx:166
 MarkovChain.cxx:167
 MarkovChain.cxx:168
 MarkovChain.cxx:169
 MarkovChain.cxx:170
 MarkovChain.cxx:171
 MarkovChain.cxx:172
 MarkovChain.cxx:173
 MarkovChain.cxx:174
 MarkovChain.cxx:175
 MarkovChain.cxx:176
 MarkovChain.cxx:177
 MarkovChain.cxx:178
 MarkovChain.cxx:179
 MarkovChain.cxx:180
 MarkovChain.cxx:181
 MarkovChain.cxx:182
 MarkovChain.cxx:183
 MarkovChain.cxx:184
 MarkovChain.cxx:185
 MarkovChain.cxx:186
 MarkovChain.cxx:187
 MarkovChain.cxx:188
 MarkovChain.cxx:189
 MarkovChain.cxx:190
 MarkovChain.cxx:191
 MarkovChain.cxx:192
 MarkovChain.cxx:193
 MarkovChain.cxx:194
 MarkovChain.cxx:195
 MarkovChain.cxx:196
 MarkovChain.cxx:197
 MarkovChain.cxx:198
 MarkovChain.cxx:199
 MarkovChain.cxx:200
 MarkovChain.cxx:201
 MarkovChain.cxx:202
 MarkovChain.cxx:203
 MarkovChain.cxx:204
 MarkovChain.cxx:205
 MarkovChain.cxx:206
 MarkovChain.cxx:207
 MarkovChain.cxx:208
 MarkovChain.cxx:209
 MarkovChain.cxx:210
 MarkovChain.cxx:211
 MarkovChain.cxx:212
 MarkovChain.cxx:213
 MarkovChain.cxx:214
 MarkovChain.cxx:215
 MarkovChain.cxx:216
 MarkovChain.cxx:217
 MarkovChain.cxx:218
 MarkovChain.cxx:219
 MarkovChain.cxx:220
 MarkovChain.cxx:221
 MarkovChain.cxx:222
 MarkovChain.cxx:223
 MarkovChain.cxx:224
 MarkovChain.cxx:225
 MarkovChain.cxx:226
 MarkovChain.cxx:227
 MarkovChain.cxx:228
 MarkovChain.cxx:229
 MarkovChain.cxx:230
 MarkovChain.cxx:231
 MarkovChain.cxx:232
 MarkovChain.cxx:233
 MarkovChain.cxx:234
 MarkovChain.cxx:235
 MarkovChain.cxx:236
 MarkovChain.cxx:237
 MarkovChain.cxx:238
 MarkovChain.cxx:239
 MarkovChain.cxx:240
 MarkovChain.cxx:241
 MarkovChain.cxx:242
 MarkovChain.cxx:243
 MarkovChain.cxx:244
 MarkovChain.cxx:245
 MarkovChain.cxx:246
 MarkovChain.cxx:247
 MarkovChain.cxx:248
 MarkovChain.cxx:249
 MarkovChain.cxx:250
 MarkovChain.cxx:251
 MarkovChain.cxx:252
 MarkovChain.cxx:253
 MarkovChain.cxx:254
 MarkovChain.cxx:255
 MarkovChain.cxx:256
 MarkovChain.cxx:257
 MarkovChain.cxx:258
 MarkovChain.cxx:259
 MarkovChain.cxx:260
 MarkovChain.cxx:261
 MarkovChain.cxx:262
 MarkovChain.cxx:263
 MarkovChain.cxx:264
 MarkovChain.cxx:265
 MarkovChain.cxx:266
 MarkovChain.cxx:267
 MarkovChain.cxx:268
 MarkovChain.cxx:269
 MarkovChain.cxx:270
 MarkovChain.cxx:271
 MarkovChain.cxx:272
 MarkovChain.cxx:273
 MarkovChain.cxx:274