From 59f987c5366e4eb4f69714869ebaf6dd4de444aa Mon Sep 17 00:00:00 2001 From: Mathieu Guigue Date: Thu, 28 Mar 2019 10:04:49 +0900 Subject: [PATCH 01/12] First non working version --- .../sampling/PyBindRooFitProcessor.py | 252 ++++++++++++++++++ morpho/processors/sampling/__init__.py | 1 + tests/sampling/sampling_test.py | 23 ++ 3 files changed, 276 insertions(+) create mode 100644 morpho/processors/sampling/PyBindRooFitProcessor.py diff --git a/morpho/processors/sampling/PyBindRooFitProcessor.py b/morpho/processors/sampling/PyBindRooFitProcessor.py new file mode 100644 index 00000000..219f055a --- /dev/null +++ b/morpho/processors/sampling/PyBindRooFitProcessor.py @@ -0,0 +1,252 @@ +''' +Processor for linear fitting +Authors: M. Guigue +Date: 06/26/18 +''' + +try: + import ROOT +except ImportError: + pass + +from morpho.utilities import morphologging, reader +from morpho.processors.sampling.RooFitInterfaceProcessor import RooFitInterfaceProcessor +from morpho.processors.BaseProcessor import BaseProcessor +logger = morphologging.getLogger(__name__) + +__all__ = [] +__all__.append(__name__) + + + +import numpy as np +import scipy.special as scs + +c = 299792458 +m_kg = 9.10938291*1e-31 +me_keV = 510.998 +q_C = 1.60217657*1e-19 +Z0 = 119.917 * np.pi +Rc = 0.06e-2 +a = 1.006e-2/2 +B0 = 1 + +def sinc(x): + if x == 0 : + return 1 + return np.sin(x)/x + +def cot(x): + return np.cos(x)/np.sin(x) + +def beta_function(n, L0, E, H, theta ): + + fc11 = 1.841 * c / (2 * np.pi * a) + gamma = 1 + E / me_keV + wc = q_C * H / ( gamma * m_kg ) + v0 = c * ( 1 - 1 / gamma ** 2 )**0.5 + vz0 = v0 * np.cos(theta) + wa = v0 / L0 * np.sin(theta) + f = wc / ( 2 * np.pi ) + Dw = 0.5 * wc * cot(theta)**2 + vp = c/(1-(2*np.pi*fc11/(wc+Dw))**2)**0.5 + k = (wc + Dw) / vp + zmax = L0 * cot( theta ) + + return scs.jv( n , k * zmax ) + +class start_freq_func(ROOT.TPyMultiGenFunction): + def __init__( self ): + print("f CREATED") + ROOT.TPyMultiGenFunction.__init__( self, self ) + + def NDim( self ): + print('PYTHON NDim Freq. called') + return 2 + def DoEval(self,args): + #(E, theta) + E = args[0] + theta = args[1] + gamma = 1 + E / me_keV + wc = q_C * B0 / ( gamma * m_kg ) + #Dw = 0.5 * wc * cot(theta)**2 + return ( wc ) / ( 2e6 * np.pi ) + +class power_function(ROOT.TPyMultiGenFunction): + def __init__( self ): + print("p CREATED") + ROOT.TPyMultiGenFunction.__init__( self, self ) + + def NDim( self ): + print('PYTHON NDim power called') + return 3 + def DoEval(self,args): + #(E, theta, L0) + E = args[0] + theta = args[1] + L0 = args[2] + return beta_function(n = 0, L0 = L0, E = E, H = B0, theta = theta) + +c1 = ROOT.TCanvas("c", "c", 800, 600) +L0 = ROOT.RooRealVar('L0', 'L0', 0.3, 0, 1) + +E = ROOT.RooRealVar('E','E',17,18) +theta = ROOT.RooRealVar('theta','theta',89*np.pi/180,np.pi/2) + + +startfreq_obj = start_freq_func() +ff01 = ROOT.TPyMultiGenFunction(startfreq_obj) +# ff0 = ROOT.RooFit.bindFunction("ff0",ff01, ROOT.RooArgList(E, theta)) + +power_obj = power_function() +fp1 = ROOT.TPyMultiGenFunction(power_obj) +# fp = ROOT.RooFit.bindFunction("fp",fp1, ROOT.RooArgList(E, theta, L0)) + + + + + + + +class PyFunctionObject(ROOT.TPyMultiGenFunction): + def __init__( self, pythonFunction, dimension=2 ): + logger.info("Created PyFunctionObject") + ROOT.TPyMultiGenFunction.__init__( self, self ) + self.pythonFunction = pythonFunction + self.dimension = dimension + + def NDim( self ): + logger.info('PYTHON NDim called') + return self.dimension + + def __call__(self,args): + # x = args[0] + # y = args[1]; + # tmp1 = y-x*x; + # tmp2 = 1-x; + return self.pythonFunction(args); + + +class PyBindRooFitProcessor(RooFitInterfaceProcessor): + ''' + Linear fit of data using RootFit Likelihood sampler. + We redefine the _defineDataset method as this analysis requires datapoints in a 2D space. + Users should feel free to change this method as they see fit. + + Parameters: + varName (required): name(s) of the variable in the data + nuisanceParams (required): parameters to be discarded at end of sampling + interestParams (required): parameters to be saved in the results variable + iter (required): total number of iterations (warmup and sampling) + warmup: number of warmup iterations (default=iter/2) + chain: number of chains (default=1) + n_jobs: number of parallel cores running (default=1) + binned: should do binned analysis (default=false) + options: other options + a (required): range of slopes (list) + b (required): range of intercepts (list) + x (required): range of x (list) + y (required): range of y (list) + witdh (required): range of width (list) + + Input: + data: dictionary containing model input data + + Results: + results: dictionary containing the result of the sampling of the parameters of interest + ''' + + def InternalConfigure(self, config_dict): + super().InternalConfigure(config_dict) + self.ranges = reader.read_param(config_dict, "paramRange", "required") + self.module_name = reader.read_param(config_dict, "module_name", "required") + self.function_name = reader.read_param(config_dict, "function_name", "required") + # Test if the module exists + try: + import imp + self.module = imp.load_source( + self.module_name, self.module_name+'.py') + except Exception as err: + logger.critical(err) + return 0 + # Test if the function exists in the file + if hasattr(self.module, self.function_name): + logger.info("Found {} using {}".format( + self.function_name, self.module_name)) + else: + logger.critical("Couldn't find {} using {}".format( + self.function_name, self.module_name)) + return False + return True + # for varName in self.varName: + # if varName not in self.ranges.keys(): + # logger.error("Missing range for {}".format(varName)) + # return False + # self.x_min, self.x_max = reader.read_param(config_dict, "paramRange", "required")["x"] + # self.mean_min, self.mean_max = reader.read_param(config_dict, "paramRange", "required")["mean"] + # self.width_min, self.width_max = reader.read_param(config_dict, "paramRange", "required")["width"] + + return True + + def _defineDataset(self, wspace): + rooVarSet = set() + + for aVarName in self.ranges.keys(): + rooVarSet.add(ROOT.RooRealVar(str(aVarName), str(aVarName), min(self._data[aVarName]), max(self._data[aVarName]))) + print(rooVarSet) + # data = ROOT.RooDataSet(self.datasetName, self.datasetName, ROOT.RooArgSet(*rooVarSet)) + # for x in self._data["x"]: + # varX.setVal(x) + # data.add(ROOT.RooArgSet(varX)) + # getattr(wspace, 'import')(data) + return wspace + + def definePdf(self, wspace): + ''' + Define the model which is that the residual of the linear fit should be normally distributed. + ''' + logger.debug("Defining pdf") + # mean = ROOT.RooRealVar("mean", "mean", 0, self.mean_min, self.mean_max) + # width = ROOT.RooRealVar("width", "width", 1., self.width_min, self.width_max) + # x = ROOT.RooRealVar("x", "x", 0, self.x_min, self.x_max) + + rooVarSet = set() + + for aVarName in self.ranges.keys(): + rooVarSet.add(ROOT.RooRealVar(str(aVarName), str(aVarName), self.ranges[aVarName][0], self.ranges[aVarName][1])) + + print(rooVarSet) + # pdf = ROOT.RooGaussian("pdf", "pdf", x, mean, width) + # try: + self.func = getattr(self.module, self.function_name) + print(self.func) + # self.pyFuncObj = PyFunctionObject(self.func) + # print(self.pyFuncObj) + self.f = ROOT.TPyMultiGenFunction(self.func) + print(self.f) + + # except: + # logger.critical("I failed") + + + self.pdf = ROOT.RooFit.bindFunction("test",self.f, ROOT.RooArgList(*rooVarSet)) + print(self.pdf) + # try: + # self.results = getattr(self.module, self.function_name)(self.config_dict) + # return True + # except Exception as err: + # logger.critical(err) + # return False + + # Save pdf: this will save all required variables and functions + getattr(wspace, 'import')(self.pdf) + + return wspace + + +if __name__ == "__main__": + rose = RosenBrock() + f = ROOT.TPyMultiGenFunction(rose) + x = ROOT.RooRealVar("x","x",0,10) + a = ROOT.RooRealVar("a","a",1,2) + fx = ROOT.RooFit.bindFunction("test",f, ROOT.RooArgList(x, a)) \ No newline at end of file diff --git a/morpho/processors/sampling/__init__.py b/morpho/processors/sampling/__init__.py index d6acd886..e64d43e1 100644 --- a/morpho/processors/sampling/__init__.py +++ b/morpho/processors/sampling/__init__.py @@ -8,4 +8,5 @@ from .PyStanSamplingProcessor import PyStanSamplingProcessor from .RooFitInterfaceProcessor import RooFitInterfaceProcessor +from .PyBindRooFitProcessor import PyBindRooFitProcessor from .LinearFitRooFitProcessor import LinearFitRooFitProcessor diff --git a/tests/sampling/sampling_test.py b/tests/sampling/sampling_test.py index 3a86038a..602079ce 100644 --- a/tests/sampling/sampling_test.py +++ b/tests/sampling/sampling_test.py @@ -93,6 +93,29 @@ def mean(numbers): return float(sum(numbers)) / max(len(numbers), 1) self.assertTrue(mean(fitterProcessor.results["a"]) > 0.5) + + def test_PyBind(self): + logger.info("PyBind tester") + from morpho.processors.sampling.PyBindRooFitProcessor import PyBindRooFitProcessor + pybind_gene_config = { + "varName": "XY", + "paramRange": { + "x": [-10, 10], + "y": [-10, 50], + "a": [-10, 10], + "b": [-10, 10] + }, + "iter": 10000, + "fixedParams": {}, + "interestParams": ['a', 'b'], + "module_name": "myModule", + "function_name": "myFunction" + } + sampler = PyBindRooFitProcessor("sampler") + self.assertTrue(sampler.Configure(pybind_gene_config)) + self.assertTrue(sampler.Run()) + self.assertTrue(1) + def test_GaussianSampler(self): logger.info("GaussianSampler test") from morpho.processors.sampling.GaussianSamplingProcessor import GaussianSamplingProcessor From f703dff5ab6c4beb3769ceed11e2d23bc077ba85 Mon Sep 17 00:00:00 2001 From: Mathieu Guigue Date: Thu, 28 Mar 2019 10:05:15 +0900 Subject: [PATCH 02/12] Adding simple model --- tests/sampling/myModule.py | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 tests/sampling/myModule.py diff --git a/tests/sampling/myModule.py b/tests/sampling/myModule.py new file mode 100644 index 00000000..b06a7cb4 --- /dev/null +++ b/tests/sampling/myModule.py @@ -0,0 +1,11 @@ +''' +A super simple module for unittest +''' + +from math import sin +from morpho.utilities import morphologging +logger = morphologging.getLogger(__name__) + +def myFunction(x,a,b): + logger.info("This is my function") + return a*sin(b*x) \ No newline at end of file From 62d5ccc0a0d607914fdfd4c3ae30e5b7ed0238cb Mon Sep 17 00:00:00 2001 From: Mathieu GUIGUE Date: Tue, 28 May 2019 12:05:15 +0200 Subject: [PATCH 03/12] Disgusting but working... --- .../sampling/PyBindRooFitProcessor.py | 198 +++++++++++++----- .../sampling/RooFitInterfaceProcessor.py | 2 + tests/sampling/myModule.py | 7 +- tests/sampling/sampling_test.py | 23 +- 4 files changed, 166 insertions(+), 64 deletions(-) diff --git a/morpho/processors/sampling/PyBindRooFitProcessor.py b/morpho/processors/sampling/PyBindRooFitProcessor.py index 219f055a..b5c0ce70 100644 --- a/morpho/processors/sampling/PyBindRooFitProcessor.py +++ b/morpho/processors/sampling/PyBindRooFitProcessor.py @@ -1,3 +1,8 @@ +import scipy.special as scs +import numpy as np +from morpho.processors.BaseProcessor import BaseProcessor +from morpho.processors.sampling.RooFitInterfaceProcessor import RooFitInterfaceProcessor +from morpho.utilities import morphologging, reader ''' Processor for linear fitting Authors: M. Guigue @@ -8,20 +13,17 @@ import ROOT except ImportError: pass +value = ROOT.gSystem.Load("libRooFit") +if value < 0: + print("Failed loading", value) + exit() -from morpho.utilities import morphologging, reader -from morpho.processors.sampling.RooFitInterfaceProcessor import RooFitInterfaceProcessor -from morpho.processors.BaseProcessor import BaseProcessor logger = morphologging.getLogger(__name__) __all__ = [] __all__.append(__name__) - -import numpy as np -import scipy.special as scs - c = 299792458 m_kg = 9.10938291*1e-31 me_keV = 510.998 @@ -31,67 +33,75 @@ a = 1.006e-2/2 B0 = 1 + def sinc(x): - if x == 0 : + if x == 0: return 1 return np.sin(x)/x + def cot(x): return np.cos(x)/np.sin(x) -def beta_function(n, L0, E, H, theta ): + +def beta_function(n, L0, E, H, theta): fc11 = 1.841 * c / (2 * np.pi * a) gamma = 1 + E / me_keV - wc = q_C * H / ( gamma * m_kg ) - v0 = c * ( 1 - 1 / gamma ** 2 )**0.5 + wc = q_C * H / (gamma * m_kg) + v0 = c * (1 - 1 / gamma ** 2)**0.5 vz0 = v0 * np.cos(theta) wa = v0 / L0 * np.sin(theta) - f = wc / ( 2 * np.pi ) + f = wc / (2 * np.pi) Dw = 0.5 * wc * cot(theta)**2 vp = c/(1-(2*np.pi*fc11/(wc+Dw))**2)**0.5 k = (wc + Dw) / vp - zmax = L0 * cot( theta ) + zmax = L0 * cot(theta) + + return scs.jv(n, k * zmax) - return scs.jv( n , k * zmax ) class start_freq_func(ROOT.TPyMultiGenFunction): - def __init__( self ): + def __init__(self): print("f CREATED") - ROOT.TPyMultiGenFunction.__init__( self, self ) + ROOT.TPyMultiGenFunction.__init__(self, self) - def NDim( self ): + def NDim(self): print('PYTHON NDim Freq. called') return 2 - def DoEval(self,args): + + def DoEval(self, args): #(E, theta) E = args[0] theta = args[1] gamma = 1 + E / me_keV - wc = q_C * B0 / ( gamma * m_kg ) + wc = q_C * B0 / (gamma * m_kg) #Dw = 0.5 * wc * cot(theta)**2 - return ( wc ) / ( 2e6 * np.pi ) + return (wc) / (2e6 * np.pi) + class power_function(ROOT.TPyMultiGenFunction): - def __init__( self ): + def __init__(self): print("p CREATED") - ROOT.TPyMultiGenFunction.__init__( self, self ) + ROOT.TPyMultiGenFunction.__init__(self, self) - def NDim( self ): + def NDim(self): print('PYTHON NDim power called') return 3 - def DoEval(self,args): + + def DoEval(self, args): #(E, theta, L0) E = args[0] theta = args[1] L0 = args[2] - return beta_function(n = 0, L0 = L0, E = E, H = B0, theta = theta) + return beta_function(n=0, L0=L0, E=E, H=B0, theta=theta) + c1 = ROOT.TCanvas("c", "c", 800, 600) L0 = ROOT.RooRealVar('L0', 'L0', 0.3, 0, 1) -E = ROOT.RooRealVar('E','E',17,18) -theta = ROOT.RooRealVar('theta','theta',89*np.pi/180,np.pi/2) +E = ROOT.RooRealVar('E', 'E', 17, 18) +theta = ROOT.RooRealVar('theta', 'theta', 89*np.pi/180, np.pi/2) startfreq_obj = start_freq_func() @@ -103,28 +113,34 @@ def DoEval(self,args): # fp = ROOT.RooFit.bindFunction("fp",fp1, ROOT.RooArgList(E, theta, L0)) - - - - - class PyFunctionObject(ROOT.TPyMultiGenFunction): - def __init__( self, pythonFunction, dimension=2 ): + def __init__(self, pythonFunction, dimension=2): logger.info("Created PyFunctionObject") - ROOT.TPyMultiGenFunction.__init__( self, self ) + ROOT.TPyMultiGenFunction.__init__(self, self) self.pythonFunction = pythonFunction self.dimension = dimension - def NDim( self ): - logger.info('PYTHON NDim called') + def NDim(self): + logger.info('PYTHON NDim called: {}'.format(self.dimension)) return self.dimension - def __call__(self,args): + def __call__(self, args): + E = args[0] + theta = args[1] + L0 = args[2] + return self.pythonFunction(E, theta, L0) + + def DoEval(self, args): + # print(args) # x = args[0] # y = args[1]; # tmp1 = y-x*x; - # tmp2 = 1-x; - return self.pythonFunction(args); + # # tmp2 = 1-x; + # return self.pythonFunction(args) + E = args[0] + theta = args[1] + L0 = args[2] + return self.pythonFunction(E, theta, L0) class PyBindRooFitProcessor(RooFitInterfaceProcessor): @@ -159,8 +175,10 @@ class PyBindRooFitProcessor(RooFitInterfaceProcessor): def InternalConfigure(self, config_dict): super().InternalConfigure(config_dict) self.ranges = reader.read_param(config_dict, "paramRange", "required") - self.module_name = reader.read_param(config_dict, "module_name", "required") - self.function_name = reader.read_param(config_dict, "function_name", "required") + self.module_name = reader.read_param( + config_dict, "module_name", "required") + self.function_name = reader.read_param( + config_dict, "function_name", "required") # Test if the module exists try: import imp @@ -190,9 +208,10 @@ def InternalConfigure(self, config_dict): def _defineDataset(self, wspace): rooVarSet = set() - + for aVarName in self.ranges.keys(): - rooVarSet.add(ROOT.RooRealVar(str(aVarName), str(aVarName), min(self._data[aVarName]), max(self._data[aVarName]))) + rooVarSet.add(ROOT.RooRealVar(str(aVarName), str(aVarName), min( + self._data[aVarName]), max(self._data[aVarName]))) print(rooVarSet) # data = ROOT.RooDataSet(self.datasetName, self.datasetName, ROOT.RooArgSet(*rooVarSet)) # for x in self._data["x"]: @@ -211,26 +230,88 @@ def definePdf(self, wspace): # x = ROOT.RooRealVar("x", "x", 0, self.x_min, self.x_max) rooVarSet = set() - + for aVarName in self.ranges.keys(): - rooVarSet.add(ROOT.RooRealVar(str(aVarName), str(aVarName), self.ranges[aVarName][0], self.ranges[aVarName][1])) + if str(aVarName) == "x": + rooVarSet.add(ROOT.RooRealVar(str(aVarName), str(aVarName), self.ranges[aVarName][0], self.ranges[aVarName][1])) + else: + rooVarSet.add(ROOT.RooRealVar(str(aVarName), str(aVarName), 1, self.ranges[aVarName][0], self.ranges[aVarName][1])) print(rooVarSet) # pdf = ROOT.RooGaussian("pdf", "pdf", x, mean, width) # try: self.func = getattr(self.module, self.function_name) print(self.func) - # self.pyFuncObj = PyFunctionObject(self.func) + # self.f = PyFunctionObject(self.func,len(self.ranges)) + self.f = PyFunctionObject(self.func, 3) # print(self.pyFuncObj) - self.f = ROOT.TPyMultiGenFunction(self.func) + # self.f = ROOT.TPyMultiGenFunction(self.func) print(self.f) + self.bindFunc = ROOT.RooFit.bindFunction("test", self.f, ROOT.RooArgList(*rooVarSet)) + # print(self.f.DoEval([1,1,1])) + + + # // Construct parameter mean2 and sigma + # mean = ROOT.RooRealVar("mean","mean",10,0,200) ; + # sigma = ROOT.RooRealVar("sigma","sigma",3,0.1,10) ; + + # # // Construct interpreted function mean = sqrt(mean^2) + # def testfunc(x, mean, sigma): + # return mean*sin(x*sigma) #1/sqrt(2*3.14159)*exp(-0.5*(x-mean)*(x-mean)) + # g2 = ROOT.RooFormulaVar("mean","mean","testfunc(x,mean,sigma)",ROOT.RooArgList(x,mean,sigma)) ; + # gen = ROOT.RooGenericPdf("gen","@0",ROOT.RooArgList(g2)) + + # // Construct a gaussian g2(x,sqrt(mean2),sigma) ; + # g2 = ROOT.RooGaussian("g2","h2",x,mean,sigma) ; + + + # // G e n e r a t e t o y d a t a + # // --------------------------------- + + # // Construct a separate gaussian g1(x,10,3) to generate a toy Gaussian dataset with mean 10 and width 3 + # g1 = ROOT.RooGaussian("g1","g1",x,ROOT.RooFit.RooConst(10),ROOT.RooFit.RooConst(3)) ; + # data2 = gen.generate(ROOT.RooArgSet(x),1000) ; + # data2.Print() + # xframe2 = x.frame(ROOT.RooFit.Title("Tailored Gaussian pdf")) ; + # data2.plotOn(xframe2) + # can = ROOT.TCanvas("can", "can", 600, 400) + # xframe2.Draw() + # can.SaveAs("test.pdf") # except: # logger.critical("I failed") + x = ROOT.RooRealVar("x", "x", 0, -5, 5) + a0 = ROOT.RooRealVar("a0","a0",0); + a0.setConstant(); + # RooChebychev will make a first order polynomial, set to a constant + bkg = ROOT.RooChebychev("a0_bkg","a0_bkg",x, ROOT.RooArgList(a0)) + # RooAbsReal *LSfcn = bindFunction(SH,m, RooArgList(deltaM,cw));//deltaM and cw are parameters + self.pdf = ROOT.RooRealSumPdf("pdf","pdf",self.bindFunc,bkg,ROOT.RooFit.RooConst(1.)) #; //combine the constant term (bkg) + + print("pdf:", self.pdf) + getattr(wspace, 'import')(self.pdf) + + pdf = wspace.pdf("pdf") + x = wspace.var("x") + a = wspace.var("a") + b = wspace.var("b") + wspace.Print() + x.Print() + a.Print() + # mean = ROOT.RooRealVar("mean", "mean", 0, -1, 1) + # width = ROOT.RooRealVar("width", "width", 1., 0, 2) - self.pdf = ROOT.RooFit.bindFunction("test",self.f, ROOT.RooArgList(*rooVarSet)) - print(self.pdf) + print(ROOT.RooArgSet(*rooVarSet)) + # pdf = ROOT.RooGaussian("pdf", "pdf", x, mean, width) + data2 = pdf.generate(ROOT.RooArgSet(*rooVarSet),10000) + data2.Print() + xframe2 = x.frame(ROOT.RooFit.Title("Tailored Gaussian pdf")) ; + data2.plotOn(xframe2) + can = ROOT.TCanvas("can", "can", 600, 400) + xframe2.Draw() + can.SaveAs("test.pdf") + # self.pdf.generate(100) # try: # self.results = getattr(self.module, self.function_name)(self.config_dict) # return True @@ -239,14 +320,23 @@ def definePdf(self, wspace): # return False # Save pdf: this will save all required variables and functions - getattr(wspace, 'import')(self.pdf) - + # getattr(wspace, 'import')(self.pdf) + # print(wspace) + # TF1 SH = new TF1("SH",signal_shape,x1,x2,nPar); // Define TF1 SH->SetParameters(a,b,c...); // set parameters - my function has 2 + # 3.) Bind the TF1 to a RooFit RooAbsReal, create a constant term, then use RooRealSumPdf (root.cern.ch/root/html/RooRealSumPdf.html 17) to create a PDF out of the function and the constant: + # RooRealVar a0("a0","a0",0); + # a0.setConstant(kTRUE); + # //RooChebychev will make a first order polynomial, set to a constant + # RooChebychev bkg("bkg","bkg",m, RooArgList(a0)); //this is now a constant + # RooAbsReal *LSfcn = bindFunction(SH,m, RooArgList(deltaM,cw));//deltaM and cw are parameters + # RooAbsPdf *LSPdf = new RooRealSumPdf("LSPdf","LSPdf",*LSfcn,bkg,c1); //combine the constant term (bkg) and the term RooAbsReal (which is a function) into a PDF. + # exit() return wspace if __name__ == "__main__": rose = RosenBrock() f = ROOT.TPyMultiGenFunction(rose) - x = ROOT.RooRealVar("x","x",0,10) - a = ROOT.RooRealVar("a","a",1,2) - fx = ROOT.RooFit.bindFunction("test",f, ROOT.RooArgList(x, a)) \ No newline at end of file + x = ROOT.RooRealVar("x", "x", 0, 10) + a = ROOT.RooRealVar("a", "a", 1, 2) + fx = ROOT.RooFit.bindFunction("test", f, ROOT.RooArgList(x, a)) diff --git a/morpho/processors/sampling/RooFitInterfaceProcessor.py b/morpho/processors/sampling/RooFitInterfaceProcessor.py index ab91d55e..986e8786 100644 --- a/morpho/processors/sampling/RooFitInterfaceProcessor.py +++ b/morpho/processors/sampling/RooFitInterfaceProcessor.py @@ -164,6 +164,7 @@ def _Generator(self): wspace = self._FixParams(wspace) pdf = wspace.pdf("pdf") paramOfInterest = self._getArgSet(wspace, self.paramOfInterestNames) + paramOfInterest.Print() data = pdf.generate(paramOfInterest, self.iter) data.Print() @@ -176,6 +177,7 @@ def _Generator(self): self.data[item].append( data.get(i).getRealValue(item)) self.data.update({"is_sample": [1]*(self.iter)}) + return True def _LikelihoodSampling(self): diff --git a/tests/sampling/myModule.py b/tests/sampling/myModule.py index b06a7cb4..57f75318 100644 --- a/tests/sampling/myModule.py +++ b/tests/sampling/myModule.py @@ -6,6 +6,7 @@ from morpho.utilities import morphologging logger = morphologging.getLogger(__name__) -def myFunction(x,a,b): - logger.info("This is my function") - return a*sin(b*x) \ No newline at end of file + +def myFunction(x, a, b): + # logger.info("This is my function: {},{},{} -> {}".format(x,a,b,abs(a*sin(b*x)))) + return abs(a*sin(b*x)) diff --git a/tests/sampling/sampling_test.py b/tests/sampling/sampling_test.py index 602079ce..49ee36f2 100644 --- a/tests/sampling/sampling_test.py +++ b/tests/sampling/sampling_test.py @@ -97,24 +97,33 @@ def mean(numbers): def test_PyBind(self): logger.info("PyBind tester") from morpho.processors.sampling.PyBindRooFitProcessor import PyBindRooFitProcessor + from morpho.processors.plots import Histogram pybind_gene_config = { "varName": "XY", "paramRange": { - "x": [-10, 10], - "y": [-10, 50], - "a": [-10, 10], - "b": [-10, 10] + "x": [-5, 5], + # "y": [-10, 50], + "a": [1, 10], + "b": [1, 10] }, "iter": 10000, - "fixedParams": {}, - "interestParams": ['a', 'b'], + "fixedParams": {'a':1, 'b':2}, + "interestParams": ['x'], "module_name": "myModule", "function_name": "myFunction" } + histo_config = { + "variables": "x", + "n_bins_x": 300, + "output_path": "plots" + } sampler = PyBindRooFitProcessor("sampler") + myhisto = Histogram("histo") self.assertTrue(sampler.Configure(pybind_gene_config)) + self.assertTrue(myhisto.Configure(histo_config)) self.assertTrue(sampler.Run()) - self.assertTrue(1) + myhisto.data = sampler.data + self.assertTrue(myhisto.Run()) def test_GaussianSampler(self): logger.info("GaussianSampler test") From ef4fc1ea4f8cf211e26e6bdd508d6778bfac441a Mon Sep 17 00:00:00 2001 From: Mathieu GUIGUE Date: Tue, 28 May 2019 12:45:53 +0200 Subject: [PATCH 04/12] Now variables are produced in the right order... --- .../sampling/PyBindRooFitProcessor.py | 84 +++++++++++-------- tests/sampling/myModule.py | 4 +- 2 files changed, 49 insertions(+), 39 deletions(-) diff --git a/morpho/processors/sampling/PyBindRooFitProcessor.py b/morpho/processors/sampling/PyBindRooFitProcessor.py index b5c0ce70..130d6b1c 100644 --- a/morpho/processors/sampling/PyBindRooFitProcessor.py +++ b/morpho/processors/sampling/PyBindRooFitProcessor.py @@ -124,11 +124,11 @@ def NDim(self): logger.info('PYTHON NDim called: {}'.format(self.dimension)) return self.dimension - def __call__(self, args): - E = args[0] - theta = args[1] - L0 = args[2] - return self.pythonFunction(E, theta, L0) + # def __call__(self, args): + # E = args[0] + # theta = args[1] + # L0 = args[2] + # return self.pythonFunction(E, theta, L0) def DoEval(self, args): # print(args) @@ -137,10 +137,17 @@ def DoEval(self, args): # tmp1 = y-x*x; # # tmp2 = 1-x; # return self.pythonFunction(args) - E = args[0] - theta = args[1] - L0 = args[2] - return self.pythonFunction(E, theta, L0) + # E = args[0] + # theta = args[1] + # L0 = args[2] + test_argv = list() + for i in range(self.dimension): + value = args[i] + test_argv.append(value) + # print("argv",*argv) + # print("normal",E,theta,L0) + return self.pythonFunction(*test_argv) + # return self.pythonFunction(E, theta, L0) class PyBindRooFitProcessor(RooFitInterfaceProcessor): @@ -204,7 +211,7 @@ def InternalConfigure(self, config_dict): # self.mean_min, self.mean_max = reader.read_param(config_dict, "paramRange", "required")["mean"] # self.width_min, self.width_max = reader.read_param(config_dict, "paramRange", "required")["width"] - return True + # return True def _defineDataset(self, wspace): rooVarSet = set() @@ -229,24 +236,27 @@ def definePdf(self, wspace): # width = ROOT.RooRealVar("width", "width", 1., self.width_min, self.width_max) # x = ROOT.RooRealVar("x", "x", 0, self.x_min, self.x_max) - rooVarSet = set() - + rooVarSet = list() + aVarSampling = 0 for aVarName in self.ranges.keys(): - if str(aVarName) == "x": - rooVarSet.add(ROOT.RooRealVar(str(aVarName), str(aVarName), self.ranges[aVarName][0], self.ranges[aVarName][1])) + print(aVarName) + + if aVarName in self.fixedParameters.keys(): + logger.debug("{} is fixed".format(aVarName)) + rooVarSet.append(ROOT.RooRealVar(str(aVarName), str(aVarName), self.fixedParameters[aVarName])) else: - rooVarSet.add(ROOT.RooRealVar(str(aVarName), str(aVarName), 1, self.ranges[aVarName][0], self.ranges[aVarName][1])) + aVarSampling = ROOT.RooRealVar(str(aVarName), str(aVarName), self.ranges[aVarName][0], self.ranges[aVarName][1]) + rooVarSet.append(aVarSampling) - print(rooVarSet) # pdf = ROOT.RooGaussian("pdf", "pdf", x, mean, width) # try: self.func = getattr(self.module, self.function_name) - print(self.func) + # print(self.func) # self.f = PyFunctionObject(self.func,len(self.ranges)) self.f = PyFunctionObject(self.func, 3) # print(self.pyFuncObj) # self.f = ROOT.TPyMultiGenFunction(self.func) - print(self.f) + # print(self.f) self.bindFunc = ROOT.RooFit.bindFunction("test", self.f, ROOT.RooArgList(*rooVarSet)) # print(self.f.DoEval([1,1,1])) @@ -284,33 +294,33 @@ def definePdf(self, wspace): a0 = ROOT.RooRealVar("a0","a0",0); a0.setConstant(); # RooChebychev will make a first order polynomial, set to a constant - bkg = ROOT.RooChebychev("a0_bkg","a0_bkg",x, ROOT.RooArgList(a0)) + bkg = ROOT.RooChebychev("a0_bkg","a0_bkg",aVarSampling, ROOT.RooArgList(a0)) # RooAbsReal *LSfcn = bindFunction(SH,m, RooArgList(deltaM,cw));//deltaM and cw are parameters self.pdf = ROOT.RooRealSumPdf("pdf","pdf",self.bindFunc,bkg,ROOT.RooFit.RooConst(1.)) #; //combine the constant term (bkg) print("pdf:", self.pdf) getattr(wspace, 'import')(self.pdf) - pdf = wspace.pdf("pdf") - x = wspace.var("x") - a = wspace.var("a") - b = wspace.var("b") + # pdf = wspace.pdf("pdf") + # x = wspace.var("x") + # a = wspace.var("a") + # b = wspace.var("b") - wspace.Print() - x.Print() - a.Print() - # mean = ROOT.RooRealVar("mean", "mean", 0, -1, 1) - # width = ROOT.RooRealVar("width", "width", 1., 0, 2) + # wspace.Print() + # x.Print() + # a.Print() + # # mean = ROOT.RooRealVar("mean", "mean", 0, -1, 1) + # # width = ROOT.RooRealVar("width", "width", 1., 0, 2) - print(ROOT.RooArgSet(*rooVarSet)) - # pdf = ROOT.RooGaussian("pdf", "pdf", x, mean, width) - data2 = pdf.generate(ROOT.RooArgSet(*rooVarSet),10000) - data2.Print() - xframe2 = x.frame(ROOT.RooFit.Title("Tailored Gaussian pdf")) ; - data2.plotOn(xframe2) - can = ROOT.TCanvas("can", "can", 600, 400) - xframe2.Draw() - can.SaveAs("test.pdf") + # print(ROOT.RooArgSet(*rooVarSet)) + # # pdf = ROOT.RooGaussian("pdf", "pdf", x, mean, width) + # data2 = pdf.generate(ROOT.RooArgSet(*rooVarSet),10000) + # data2.Print() + # xframe2 = x.frame(ROOT.RooFit.Title("Tailored Gaussian pdf")) ; + # data2.plotOn(xframe2) + # can = ROOT.TCanvas("can", "can", 600, 400) + # xframe2.Draw() + # can.SaveAs("test.pdf") # self.pdf.generate(100) # try: # self.results = getattr(self.module, self.function_name)(self.config_dict) diff --git a/tests/sampling/myModule.py b/tests/sampling/myModule.py index 57f75318..2e990499 100644 --- a/tests/sampling/myModule.py +++ b/tests/sampling/myModule.py @@ -2,11 +2,11 @@ A super simple module for unittest ''' -from math import sin +from math import sin, cos from morpho.utilities import morphologging logger = morphologging.getLogger(__name__) def myFunction(x, a, b): # logger.info("This is my function: {},{},{} -> {}".format(x,a,b,abs(a*sin(b*x)))) - return abs(a*sin(b*x)) + return abs(a*cos(b*x)) From 7f8d61de9713338669fdadd2c584edb238a68f27 Mon Sep 17 00:00:00 2001 From: Mathieu GUIGUE Date: Tue, 28 May 2019 12:56:40 +0200 Subject: [PATCH 05/12] My Module is now a tritium spectrum!! :) --- .../sampling/PyBindRooFitProcessor.py | 84 +------------------ tests/sampling/myModule.py | 9 +- tests/sampling/sampling_test.py | 9 +- 3 files changed, 11 insertions(+), 91 deletions(-) diff --git a/morpho/processors/sampling/PyBindRooFitProcessor.py b/morpho/processors/sampling/PyBindRooFitProcessor.py index 130d6b1c..2eed58c5 100644 --- a/morpho/processors/sampling/PyBindRooFitProcessor.py +++ b/morpho/processors/sampling/PyBindRooFitProcessor.py @@ -121,7 +121,6 @@ def __init__(self, pythonFunction, dimension=2): self.dimension = dimension def NDim(self): - logger.info('PYTHON NDim called: {}'.format(self.dimension)) return self.dimension # def __call__(self, args): @@ -219,7 +218,6 @@ def _defineDataset(self, wspace): for aVarName in self.ranges.keys(): rooVarSet.add(ROOT.RooRealVar(str(aVarName), str(aVarName), min( self._data[aVarName]), max(self._data[aVarName]))) - print(rooVarSet) # data = ROOT.RooDataSet(self.datasetName, self.datasetName, ROOT.RooArgSet(*rooVarSet)) # for x in self._data["x"]: # varX.setVal(x) @@ -239,8 +237,6 @@ def definePdf(self, wspace): rooVarSet = list() aVarSampling = 0 for aVarName in self.ranges.keys(): - print(aVarName) - if aVarName in self.fixedParameters.keys(): logger.debug("{} is fixed".format(aVarName)) rooVarSet.append(ROOT.RooRealVar(str(aVarName), str(aVarName), self.fixedParameters[aVarName])) @@ -248,49 +244,11 @@ def definePdf(self, wspace): aVarSampling = ROOT.RooRealVar(str(aVarName), str(aVarName), self.ranges[aVarName][0], self.ranges[aVarName][1]) rooVarSet.append(aVarSampling) - # pdf = ROOT.RooGaussian("pdf", "pdf", x, mean, width) - # try: self.func = getattr(self.module, self.function_name) - # print(self.func) - # self.f = PyFunctionObject(self.func,len(self.ranges)) - self.f = PyFunctionObject(self.func, 3) - # print(self.pyFuncObj) - # self.f = ROOT.TPyMultiGenFunction(self.func) - # print(self.f) + self.f = PyFunctionObject(self.func, len(rooVarSet)) self.bindFunc = ROOT.RooFit.bindFunction("test", self.f, ROOT.RooArgList(*rooVarSet)) - # print(self.f.DoEval([1,1,1])) - - - # // Construct parameter mean2 and sigma - # mean = ROOT.RooRealVar("mean","mean",10,0,200) ; - # sigma = ROOT.RooRealVar("sigma","sigma",3,0.1,10) ; - - # # // Construct interpreted function mean = sqrt(mean^2) - # def testfunc(x, mean, sigma): - # return mean*sin(x*sigma) #1/sqrt(2*3.14159)*exp(-0.5*(x-mean)*(x-mean)) - # g2 = ROOT.RooFormulaVar("mean","mean","testfunc(x,mean,sigma)",ROOT.RooArgList(x,mean,sigma)) ; - # gen = ROOT.RooGenericPdf("gen","@0",ROOT.RooArgList(g2)) - # // Construct a gaussian g2(x,sqrt(mean2),sigma) ; - # g2 = ROOT.RooGaussian("g2","h2",x,mean,sigma) ; - - # // G e n e r a t e t o y d a t a - # // --------------------------------- - - # // Construct a separate gaussian g1(x,10,3) to generate a toy Gaussian dataset with mean 10 and width 3 - # g1 = ROOT.RooGaussian("g1","g1",x,ROOT.RooFit.RooConst(10),ROOT.RooFit.RooConst(3)) ; - # data2 = gen.generate(ROOT.RooArgSet(x),1000) ; - # data2.Print() - # xframe2 = x.frame(ROOT.RooFit.Title("Tailored Gaussian pdf")) ; - # data2.plotOn(xframe2) - # can = ROOT.TCanvas("can", "can", 600, 400) - # xframe2.Draw() - # can.SaveAs("test.pdf") - - # except: - # logger.critical("I failed") - x = ROOT.RooRealVar("x", "x", 0, -5, 5) a0 = ROOT.RooRealVar("a0","a0",0); a0.setConstant(); # RooChebychev will make a first order polynomial, set to a constant @@ -301,46 +259,6 @@ def definePdf(self, wspace): print("pdf:", self.pdf) getattr(wspace, 'import')(self.pdf) - # pdf = wspace.pdf("pdf") - # x = wspace.var("x") - # a = wspace.var("a") - # b = wspace.var("b") - - # wspace.Print() - # x.Print() - # a.Print() - # # mean = ROOT.RooRealVar("mean", "mean", 0, -1, 1) - # # width = ROOT.RooRealVar("width", "width", 1., 0, 2) - - # print(ROOT.RooArgSet(*rooVarSet)) - # # pdf = ROOT.RooGaussian("pdf", "pdf", x, mean, width) - # data2 = pdf.generate(ROOT.RooArgSet(*rooVarSet),10000) - # data2.Print() - # xframe2 = x.frame(ROOT.RooFit.Title("Tailored Gaussian pdf")) ; - # data2.plotOn(xframe2) - # can = ROOT.TCanvas("can", "can", 600, 400) - # xframe2.Draw() - # can.SaveAs("test.pdf") - # self.pdf.generate(100) - # try: - # self.results = getattr(self.module, self.function_name)(self.config_dict) - # return True - # except Exception as err: - # logger.critical(err) - # return False - - # Save pdf: this will save all required variables and functions - # getattr(wspace, 'import')(self.pdf) - # print(wspace) - # TF1 SH = new TF1("SH",signal_shape,x1,x2,nPar); // Define TF1 SH->SetParameters(a,b,c...); // set parameters - my function has 2 - # 3.) Bind the TF1 to a RooFit RooAbsReal, create a constant term, then use RooRealSumPdf (root.cern.ch/root/html/RooRealSumPdf.html 17) to create a PDF out of the function and the constant: - # RooRealVar a0("a0","a0",0); - # a0.setConstant(kTRUE); - # //RooChebychev will make a first order polynomial, set to a constant - # RooChebychev bkg("bkg","bkg",m, RooArgList(a0)); //this is now a constant - # RooAbsReal *LSfcn = bindFunction(SH,m, RooArgList(deltaM,cw));//deltaM and cw are parameters - # RooAbsPdf *LSPdf = new RooRealSumPdf("LSPdf","LSPdf",*LSfcn,bkg,c1); //combine the constant term (bkg) and the term RooAbsReal (which is a function) into a PDF. - # exit() return wspace diff --git a/tests/sampling/myModule.py b/tests/sampling/myModule.py index 2e990499..6e6ddf51 100644 --- a/tests/sampling/myModule.py +++ b/tests/sampling/myModule.py @@ -2,11 +2,12 @@ A super simple module for unittest ''' -from math import sin, cos +from math import sin, cos, sqrt, pow from morpho.utilities import morphologging logger = morphologging.getLogger(__name__) -def myFunction(x, a, b): - # logger.info("This is my function: {},{},{} -> {}".format(x,a,b,abs(a*sin(b*x)))) - return abs(a*cos(b*x)) +def myFunction(x, a, b, c): + # logger.info("This is my function: {},{},{},{} -> {}".format(x,a,b,c,abs(a*sin(b*x)))) + # return abs(a*cos(b*x)+c) + return (x Date: Tue, 28 May 2019 20:37:08 +0200 Subject: [PATCH 06/12] for some very very strange reasons, the fitting works super well but the lsampling doesn't find the right values for the example... need more cleanup before releasing... --- morpho/processors/plots/TimeSeries.py | 8 +- .../sampling/PyBindRooFitProcessor.py | 25 +++-- .../sampling/RooFitInterfaceProcessor.py | 23 +++- tests/sampling/myModule.py | 4 +- tests/sampling/sampling_test.py | 103 ++++++++++++++++-- 5 files changed, 140 insertions(+), 23 deletions(-) diff --git a/morpho/processors/plots/TimeSeries.py b/morpho/processors/plots/TimeSeries.py index dd139677..e8d783cb 100644 --- a/morpho/processors/plots/TimeSeries.py +++ b/morpho/processors/plots/TimeSeries.py @@ -75,6 +75,9 @@ def InternalRun(self): pass # Histograms must still be in memory when the pdf is saved for iName, name in enumerate(self.namedata): + if len(self.data[name])==0: + logger.warning("No data for variable {}".format(name)) + continue self.rootcanvas.cd(iName+1) listGraph.append(ROOT.TGraph()) listGraphWarmup.append(ROOT.TGraph()) @@ -91,7 +94,10 @@ def InternalRun(self): iWarmup += 1 listGraph[iName].Draw("AP") listGraph[iName].SetMarkerStyle(7) - listGraphWarmup[iName].Draw("sameP") + if iSample==0: + listGraphWarmup[iName].Draw("AP") + else: + listGraphWarmup[iName].Draw("sameP") listGraphWarmup[iName].SetMarkerStyle(7) listGraphWarmup[iName].SetMarkerColor(2) listGraph[iName].SetTitle(";Iteration;{}".format(name)) diff --git a/morpho/processors/sampling/PyBindRooFitProcessor.py b/morpho/processors/sampling/PyBindRooFitProcessor.py index 2eed58c5..3523cbb9 100644 --- a/morpho/processors/sampling/PyBindRooFitProcessor.py +++ b/morpho/processors/sampling/PyBindRooFitProcessor.py @@ -159,6 +159,8 @@ class PyBindRooFitProcessor(RooFitInterfaceProcessor): varName (required): name(s) of the variable in the data nuisanceParams (required): parameters to be discarded at end of sampling interestParams (required): parameters to be saved in the results variable + paramRange (required): range of parameters (defined as <{'a': [a_min, a_max]}>) + initValues: initial value (defined as <{'a': a_init}>) iter (required): total number of iterations (warmup and sampling) warmup: number of warmup iterations (default=iter/2) chain: number of chains (default=1) @@ -181,6 +183,7 @@ class PyBindRooFitProcessor(RooFitInterfaceProcessor): def InternalConfigure(self, config_dict): super().InternalConfigure(config_dict) self.ranges = reader.read_param(config_dict, "paramRange", "required") + self.initParamValues = reader.read_param(config_dict, "initValues", dict()) self.module_name = reader.read_param( config_dict, "module_name", "required") self.function_name = reader.read_param( @@ -212,18 +215,18 @@ def InternalConfigure(self, config_dict): # return True - def _defineDataset(self, wspace): - rooVarSet = set() + # def _defineDataset(self, wspace): + # rooVarSet = set() - for aVarName in self.ranges.keys(): - rooVarSet.add(ROOT.RooRealVar(str(aVarName), str(aVarName), min( - self._data[aVarName]), max(self._data[aVarName]))) + # for aVarName in self.ranges.keys(): + # rooVarSet.add(ROOT.RooRealVar(str(aVarName), str(aVarName), min( + # self._data[aVarName]), max(self._data[aVarName]))) # data = ROOT.RooDataSet(self.datasetName, self.datasetName, ROOT.RooArgSet(*rooVarSet)) # for x in self._data["x"]: # varX.setVal(x) # data.add(ROOT.RooArgSet(varX)) # getattr(wspace, 'import')(data) - return wspace + # return wspace def definePdf(self, wspace): ''' @@ -237,12 +240,19 @@ def definePdf(self, wspace): rooVarSet = list() aVarSampling = 0 for aVarName in self.ranges.keys(): + logger.info(aVarName) if aVarName in self.fixedParameters.keys(): logger.debug("{} is fixed".format(aVarName)) rooVarSet.append(ROOT.RooRealVar(str(aVarName), str(aVarName), self.fixedParameters[aVarName])) + logger.info(aVarName) + elif aVarName in self.initParamValues.keys(): + aVarSampling = ROOT.RooRealVar(str(aVarName), str(aVarName), self.initParamValues[aVarName], self.ranges[aVarName][0], self.ranges[aVarName][1]) + rooVarSet.append(aVarSampling) + logger.info(aVarName) else: aVarSampling = ROOT.RooRealVar(str(aVarName), str(aVarName), self.ranges[aVarName][0], self.ranges[aVarName][1]) rooVarSet.append(aVarSampling) + logger.info(aVarName) self.func = getattr(self.module, self.function_name) self.f = PyFunctionObject(self.func, len(rooVarSet)) @@ -256,7 +266,8 @@ def definePdf(self, wspace): # RooAbsReal *LSfcn = bindFunction(SH,m, RooArgList(deltaM,cw));//deltaM and cw are parameters self.pdf = ROOT.RooRealSumPdf("pdf","pdf",self.bindFunc,bkg,ROOT.RooFit.RooConst(1.)) #; //combine the constant term (bkg) - print("pdf:", self.pdf) + logger.debug("pdf: {}".format(self.pdf)) + wspace.Print() getattr(wspace, 'import')(self.pdf) return wspace diff --git a/morpho/processors/sampling/RooFitInterfaceProcessor.py b/morpho/processors/sampling/RooFitInterfaceProcessor.py index 986e8786..d8d5c005 100644 --- a/morpho/processors/sampling/RooFitInterfaceProcessor.py +++ b/morpho/processors/sampling/RooFitInterfaceProcessor.py @@ -60,7 +60,7 @@ def _defineDataset(self, wspace): self.datasetName, self.datasetName, ROOT.RooArgSet(var)) for value in self._data[self.varName]: var.setVal(value) - data.add(ROOT.RooArgSet(var)) + data.add(ROOT.RooArgSet(var)) getattr(wspace, 'import')(data) return wspace @@ -84,6 +84,7 @@ def data(self, value): def _getArgSet(self, wspace, listNames): argSet = ROOT.RooArgSet() for name in listNames: + logger.debug("{} -> {}".format(name, wspace.var(name))) argSet.add(wspace.var(name)) return argSet @@ -91,11 +92,14 @@ def InternalConfigure(self, config_dict): self.varName = reader.read_param(config_dict, "varName", "required") self.mode = reader.read_param(config_dict, "mode", "generate") logger.debug("Mode {}".format(self.mode)) - self.iter = int(reader.read_param(config_dict, "iter", "required")) + if self.mode == "lsampling" or self.mode == "generate": + self.iter = int(reader.read_param(config_dict, "iter", "required")) if self.mode == "lsampling": self.binned = int(reader.read_param(config_dict, "binned", False)) self.nuisanceParametersNames = reader.read_param(config_dict, "nuisanceParams", "required") self.warmup = int(reader.read_param(config_dict, "warmup", self.iter/2.)) + elif self.mode == "fit": + self.binned = int(reader.read_param(config_dict, "binned", False)) self.numCPU = int(reader.read_param(config_dict, "n_jobs", 1)) self.options = reader.read_param(config_dict, "options", dict()) if self.mode not in ['generate', 'lsampling', 'fit']: @@ -132,10 +136,20 @@ def _Fit(self): wspace = self._FixParams(wspace) pdf = wspace.pdf("pdf") dataset = wspace.data(self.datasetName) + dataset.Print() paramOfInterest = self._getArgSet(wspace, self.paramOfInterestNames) result = pdf.fitTo(dataset, ROOT.RooFit.Save()) result.Print() + + can = ROOT.TCanvas("can","can",600,400) + var = wspace.var(self.varName) + frame = var.frame() + dataset.plotOn(frame) + pdf.plotOn(frame) + frame.Draw() + can.SaveAs("plots/results_fit.pdf") + self.result = {} for varName in self.paramOfInterestNames: self.result.update({str(varName): wspace.var(str(varName)).getVal()}) @@ -148,6 +162,7 @@ def _FixParams(self, wspace): logger.debug("No fixed parameters given") return wspace for varName, value in self.fixedParameters.items(): + logger.debug(varName, value) wspace.var(str(varName)).setVal(float(value)) wspace.var(str(varName)).setConstant() logger.debug("Value of {} set to {}".format(varName, wspace.var(str(varName)).getVal())) @@ -177,7 +192,7 @@ def _Generator(self): self.data[item].append( data.get(i).getRealValue(item)) self.data.update({"is_sample": [1]*(self.iter)}) - + return True def _LikelihoodSampling(self): @@ -186,8 +201,8 @@ def _LikelihoodSampling(self): ''' wspace = ROOT.RooWorkspace() wspace = self._defineDataset(wspace) - wspace = self._FixParams(wspace) wspace = self.definePdf(wspace) + wspace = self._FixParams(wspace) logger.debug("Workspace content:") wspace.Print() paramOfInterest = self._getArgSet(wspace, self.paramOfInterestNames) diff --git a/tests/sampling/myModule.py b/tests/sampling/myModule.py index 6e6ddf51..2bf669b3 100644 --- a/tests/sampling/myModule.py +++ b/tests/sampling/myModule.py @@ -9,5 +9,5 @@ def myFunction(x, a, b, c): # logger.info("This is my function: {},{},{},{} -> {}".format(x,a,b,c,abs(a*sin(b*x)))) - # return abs(a*cos(b*x)+c) - return (x Date: Wed, 5 Jun 2019 12:09:50 +0200 Subject: [PATCH 07/12] Missing dependency --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index b973853b..82c8c764 100644 --- a/setup.py +++ b/setup.py @@ -22,7 +22,7 @@ extras_require = { 'core': ['uproot>=2.8.13', 'colorlog', 'PyYAML>=3.13', 'pyparsing>=2.1.5', 'pystan==2.17.1.0', 'dnspython==1.12.0', - 'pbr==0.10.8', 'cycler==0.10.0', 'lz4', 'six', 'asteval'], + 'pbr==0.10.8', 'cycler==0.10.0', 'lz4', 'six', 'asteval', 'scipy'], 'doc': ['sphinx', 'sphinx_rtd_theme', 'sphinxcontrib-programoutput', 'six', 'colorlog'] } From ad59c3def4e853e526d1a43d381b39a2f17c195d Mon Sep 17 00:00:00 2001 From: Mathieu GUIGUE Date: Tue, 11 Jun 2019 09:37:17 +0200 Subject: [PATCH 08/12] Some cleaning up; need to install scipy at a deeper level --- .../sampling/PyBindRooFitProcessor.py | 122 ------------------ .../sampling/RooFitInterfaceProcessor.py | 1 - setup.py | 4 +- tests/sampling/myModule.py | 2 - tests/sampling/sampling_test.py | 106 +++++++-------- tests/test.sh | 4 + 6 files changed, 59 insertions(+), 180 deletions(-) diff --git a/morpho/processors/sampling/PyBindRooFitProcessor.py b/morpho/processors/sampling/PyBindRooFitProcessor.py index 3523cbb9..0024a659 100644 --- a/morpho/processors/sampling/PyBindRooFitProcessor.py +++ b/morpho/processors/sampling/PyBindRooFitProcessor.py @@ -9,110 +9,6 @@ Date: 06/26/18 ''' -try: - import ROOT -except ImportError: - pass -value = ROOT.gSystem.Load("libRooFit") -if value < 0: - print("Failed loading", value) - exit() - -logger = morphologging.getLogger(__name__) - -__all__ = [] -__all__.append(__name__) - - -c = 299792458 -m_kg = 9.10938291*1e-31 -me_keV = 510.998 -q_C = 1.60217657*1e-19 -Z0 = 119.917 * np.pi -Rc = 0.06e-2 -a = 1.006e-2/2 -B0 = 1 - - -def sinc(x): - if x == 0: - return 1 - return np.sin(x)/x - - -def cot(x): - return np.cos(x)/np.sin(x) - - -def beta_function(n, L0, E, H, theta): - - fc11 = 1.841 * c / (2 * np.pi * a) - gamma = 1 + E / me_keV - wc = q_C * H / (gamma * m_kg) - v0 = c * (1 - 1 / gamma ** 2)**0.5 - vz0 = v0 * np.cos(theta) - wa = v0 / L0 * np.sin(theta) - f = wc / (2 * np.pi) - Dw = 0.5 * wc * cot(theta)**2 - vp = c/(1-(2*np.pi*fc11/(wc+Dw))**2)**0.5 - k = (wc + Dw) / vp - zmax = L0 * cot(theta) - - return scs.jv(n, k * zmax) - - -class start_freq_func(ROOT.TPyMultiGenFunction): - def __init__(self): - print("f CREATED") - ROOT.TPyMultiGenFunction.__init__(self, self) - - def NDim(self): - print('PYTHON NDim Freq. called') - return 2 - - def DoEval(self, args): - #(E, theta) - E = args[0] - theta = args[1] - gamma = 1 + E / me_keV - wc = q_C * B0 / (gamma * m_kg) - #Dw = 0.5 * wc * cot(theta)**2 - return (wc) / (2e6 * np.pi) - - -class power_function(ROOT.TPyMultiGenFunction): - def __init__(self): - print("p CREATED") - ROOT.TPyMultiGenFunction.__init__(self, self) - - def NDim(self): - print('PYTHON NDim power called') - return 3 - - def DoEval(self, args): - #(E, theta, L0) - E = args[0] - theta = args[1] - L0 = args[2] - return beta_function(n=0, L0=L0, E=E, H=B0, theta=theta) - - -c1 = ROOT.TCanvas("c", "c", 800, 600) -L0 = ROOT.RooRealVar('L0', 'L0', 0.3, 0, 1) - -E = ROOT.RooRealVar('E', 'E', 17, 18) -theta = ROOT.RooRealVar('theta', 'theta', 89*np.pi/180, np.pi/2) - - -startfreq_obj = start_freq_func() -ff01 = ROOT.TPyMultiGenFunction(startfreq_obj) -# ff0 = ROOT.RooFit.bindFunction("ff0",ff01, ROOT.RooArgList(E, theta)) - -power_obj = power_function() -fp1 = ROOT.TPyMultiGenFunction(power_obj) -# fp = ROOT.RooFit.bindFunction("fp",fp1, ROOT.RooArgList(E, theta, L0)) - - class PyFunctionObject(ROOT.TPyMultiGenFunction): def __init__(self, pythonFunction, dimension=2): logger.info("Created PyFunctionObject") @@ -123,30 +19,12 @@ def __init__(self, pythonFunction, dimension=2): def NDim(self): return self.dimension - # def __call__(self, args): - # E = args[0] - # theta = args[1] - # L0 = args[2] - # return self.pythonFunction(E, theta, L0) - def DoEval(self, args): - # print(args) - # x = args[0] - # y = args[1]; - # tmp1 = y-x*x; - # # tmp2 = 1-x; - # return self.pythonFunction(args) - # E = args[0] - # theta = args[1] - # L0 = args[2] test_argv = list() for i in range(self.dimension): value = args[i] test_argv.append(value) - # print("argv",*argv) - # print("normal",E,theta,L0) return self.pythonFunction(*test_argv) - # return self.pythonFunction(E, theta, L0) class PyBindRooFitProcessor(RooFitInterfaceProcessor): diff --git a/morpho/processors/sampling/RooFitInterfaceProcessor.py b/morpho/processors/sampling/RooFitInterfaceProcessor.py index d8d5c005..e3fef2e5 100644 --- a/morpho/processors/sampling/RooFitInterfaceProcessor.py +++ b/morpho/processors/sampling/RooFitInterfaceProcessor.py @@ -162,7 +162,6 @@ def _FixParams(self, wspace): logger.debug("No fixed parameters given") return wspace for varName, value in self.fixedParameters.items(): - logger.debug(varName, value) wspace.var(str(varName)).setVal(float(value)) wspace.var(str(varName)).setConstant() logger.debug("Value of {} set to {}".format(varName, wspace.var(str(varName)).getVal())) diff --git a/setup.py b/setup.py index 82c8c764..869ee0c3 100644 --- a/setup.py +++ b/setup.py @@ -21,8 +21,8 @@ requirements = [] extras_require = { 'core': ['uproot>=2.8.13', 'colorlog', 'PyYAML>=3.13', 'pyparsing>=2.1.5', - 'pystan==2.17.1.0', 'dnspython==1.12.0', - 'pbr==0.10.8', 'cycler==0.10.0', 'lz4', 'six', 'asteval', 'scipy'], + 'pystan==2.19.0.0', 'dnspython==1.12.0', + 'pbr==0.10.8', 'cycler==0.10.0', 'lz4', 'six', 'asteval'], 'doc': ['sphinx', 'sphinx_rtd_theme', 'sphinxcontrib-programoutput', 'six', 'colorlog'] } diff --git a/tests/sampling/myModule.py b/tests/sampling/myModule.py index 2bf669b3..47e7f514 100644 --- a/tests/sampling/myModule.py +++ b/tests/sampling/myModule.py @@ -8,6 +8,4 @@ def myFunction(x, a, b, c): - # logger.info("This is my function: {},{},{},{} -> {}".format(x,a,b,c,abs(a*sin(b*x)))) return abs(cos(b*x)+c) - # return (x Date: Thu, 8 Aug 2019 12:51:46 -0700 Subject: [PATCH 09/12] Re-added importing of ROOT to get the PyBind test working again. --- morpho/processors/sampling/PyBindRooFitProcessor.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/morpho/processors/sampling/PyBindRooFitProcessor.py b/morpho/processors/sampling/PyBindRooFitProcessor.py index 0024a659..1d9e8996 100644 --- a/morpho/processors/sampling/PyBindRooFitProcessor.py +++ b/morpho/processors/sampling/PyBindRooFitProcessor.py @@ -9,6 +9,18 @@ Date: 06/26/18 ''' +try: + import ROOT +except ImportError: + pass +value = ROOT.gSystem.Load("libRooFit") +if value < 0: + print("Failed loading", value) + exit() + +logger = morphologging.getLogger(__name__) + + class PyFunctionObject(ROOT.TPyMultiGenFunction): def __init__(self, pythonFunction, dimension=2): logger.info("Created PyFunctionObject") From caf67cd216b5a2ea04c9263add875d9dc42a5f95 Mon Sep 17 00:00:00 2001 From: Mathieu GUIGUE Date: Mon, 26 Aug 2019 20:28:09 +0200 Subject: [PATCH 10/12] Reversing pystan version for testing --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 869ee0c3..b973853b 100644 --- a/setup.py +++ b/setup.py @@ -21,7 +21,7 @@ requirements = [] extras_require = { 'core': ['uproot>=2.8.13', 'colorlog', 'PyYAML>=3.13', 'pyparsing>=2.1.5', - 'pystan==2.19.0.0', 'dnspython==1.12.0', + 'pystan==2.17.1.0', 'dnspython==1.12.0', 'pbr==0.10.8', 'cycler==0.10.0', 'lz4', 'six', 'asteval'], 'doc': ['sphinx', 'sphinx_rtd_theme', 'sphinxcontrib-programoutput', 'six', 'colorlog'] } From 22d8991fb8ca1da3f01940c109f540ca58cf47da Mon Sep 17 00:00:00 2001 From: Mathieu GUIGUE Date: Mon, 26 Aug 2019 20:36:46 +0200 Subject: [PATCH 11/12] Removing pystan upgrade before tests --- tests/test.sh | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test.sh b/tests/test.sh index f199f152..53347177 100644 --- a/tests/test.sh +++ b/tests/test.sh @@ -2,7 +2,6 @@ # missing scipy (required only for testing) pip3 install scipy -pip3 install pystan --upgrade --no-cache cd IO && python IO_test.py && cd .. cd misc && python misc_test.py && cd .. From cb045727680bbea4119a820f832a86f2745aa08c Mon Sep 17 00:00:00 2001 From: Mathieu GUIGUE Date: Mon, 26 Aug 2019 21:03:47 +0200 Subject: [PATCH 12/12] Missed a merge conflict... --- morpho/processors/sampling/RooFitInterfaceProcessor.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/morpho/processors/sampling/RooFitInterfaceProcessor.py b/morpho/processors/sampling/RooFitInterfaceProcessor.py index 0559f17d..2eb24cd2 100644 --- a/morpho/processors/sampling/RooFitInterfaceProcessor.py +++ b/morpho/processors/sampling/RooFitInterfaceProcessor.py @@ -157,15 +157,6 @@ def _Fit(self): result = pdf.fitTo(dataset, ROOT.RooFit.Save()) result.Print() -<<<<<<< HEAD - can = ROOT.TCanvas("can","can",600,400) - var = wspace.var(self.varName) - frame = var.frame() - dataset.plotOn(frame) - pdf.plotOn(frame) - frame.Draw() - can.SaveAs("plots/results_fit.pdf") -======= if self.make_fit_plot: can = ROOT.TCanvas("can", "can", 600, 400) var = wspace.var(self.varName) @@ -174,7 +165,6 @@ def _Fit(self): pdf.plotOn(frame) frame.Draw() can.SaveAs("results_fit.pdf") ->>>>>>> develop self.result = {} for varName in self.paramOfInterestNames: