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 new file mode 100644 index 00000000..1d9e8996 --- /dev/null +++ b/morpho/processors/sampling/PyBindRooFitProcessor.py @@ -0,0 +1,171 @@ +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 +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") + ROOT.TPyMultiGenFunction.__init__(self, self) + self.pythonFunction = pythonFunction + self.dimension = dimension + + def NDim(self): + return self.dimension + + def DoEval(self, args): + test_argv = list() + for i in range(self.dimension): + value = args[i] + test_argv.append(value) + return self.pythonFunction(*test_argv) + + +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 + 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) + 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.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( + 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]))) + # 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 = 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)) + self.bindFunc = ROOT.RooFit.bindFunction("test", self.f, ROOT.RooArgList(*rooVarSet)) + + + 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",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) + + logger.debug("pdf: {}".format(self.pdf)) + wspace.Print() + 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)) diff --git a/morpho/processors/sampling/RooFitInterfaceProcessor.py b/morpho/processors/sampling/RooFitInterfaceProcessor.py index 13c55c3f..2eb24cd2 100644 --- a/morpho/processors/sampling/RooFitInterfaceProcessor.py +++ b/morpho/processors/sampling/RooFitInterfaceProcessor.py @@ -66,7 +66,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) logger.info("Workspace after dataset:") wspace.Print() @@ -92,6 +92,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 @@ -150,6 +151,7 @@ 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()) @@ -198,6 +200,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() @@ -210,6 +213,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): @@ -218,8 +222,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/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/myModule.py b/tests/sampling/myModule.py new file mode 100644 index 00000000..47e7f514 --- /dev/null +++ b/tests/sampling/myModule.py @@ -0,0 +1,11 @@ +''' +A super simple module for unittest +''' + +from math import sin, cos, sqrt, pow +from morpho.utilities import morphologging +logger = morphologging.getLogger(__name__) + + +def myFunction(x, a, b, c): + return abs(cos(b*x)+c) diff --git a/tests/sampling/sampling_test.py b/tests/sampling/sampling_test.py index 63a972c7..683b62bc 100644 --- a/tests/sampling/sampling_test.py +++ b/tests/sampling/sampling_test.py @@ -116,6 +116,124 @@ 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 + from morpho.processors.plots import Histogram, APosterioriDistribution, TimeSeries + pybind_gene_config = { + "varName": "XY", + "paramRange": { + # "x": [17000, 20000], + "x": [-10, 10], + # "y": [-10, 50], + "a": [0,10], + "b": [0, 10], + "c": [1, 10000] + }, + "iter": 10000, + "fixedParams": {'a':1, 'b':1,'c':2}, + "interestParams": ['x'], + "module_name": "myModule", + "function_name": "myFunction", + "mode": "generate" + } + pybind_fit_config = { + "varName": "x", + "paramRange": { + "x": [-10, 10], + # "y": [-10, 50], + "a": [0, 5], + # "a": [18000, 19000], + "b": [0, 5], + "c": [0, 4] + }, + "initValues": { + # "y": [-10, 50], + "a": 1, + # "a": [18000, 19000], + "b": 1, + "c": 2 + }, + # "iter": 10000, + "fixedParams": {}, + "interestParams": ['a','b','c'], + "module_name": "myModule", + "function_name": "myFunction", + "binned": True, + "mode": "fit" + } + # pybind_lsampling_config = { + # "varName": "x", + # "paramRange": { + # "x": [-10, 10], + # # "y": [-10, 50], + # "a": [0,4], + # # "a": [18000, 19000], + # "b": [0.1, 4], + # "c": [1, 4] + # }, + # "initValues": { + # # "y": [-10, 50], + # # "a": 1, + # # "a": [18000, 19000], + # "b": 1 + # # "c": 2 + # }, + # "iter": 5000, + # "warmup": 1000, + # "fixedParams": {'a':1, 'c':2}, + # "interestParams": ['b'], + # "nuisanceParams": [], + # "module_name": "myModule", + # "function_name": "myFunction", + # "binned": False, + # "n_jobs": 3, + # "mode": "lsampling" + # } + histo_config = { + "variables": "x", + "n_bins_x": 300, + "output_path": "plots" + } + # aposteriori_config = { + # "n_bins_x": 100, + # "n_bins_y": 100, + # "variables": ['b', 'lp_prob'], + # "title": "aposteriori_distribution", + # "output_path": "plots" + # } + # timeSeries_config = { + # "variables": ['b', 'lp_prob'], + # "height": 600, + # "title": "timeseries", + # "output_path": "plots" + # } + sampler = PyBindRooFitProcessor("sampler") + fitter = PyBindRooFitProcessor("fitter") + # lsampler = PyBindRooFitProcessor("lsampler") + myhisto = Histogram("histo") + # timeSeriesPlotter = TimeSeries("timeSeries") + # aposterioriPlotter = APosterioriDistribution("posterioriDistrib") + self.assertTrue(sampler.Configure(pybind_gene_config)) + self.assertTrue(fitter.Configure(pybind_fit_config)) + # self.assertTrue(lsampler.Configure(pybind_lsampling_config)) + # self.assertTrue(timeSeriesPlotter.Configure(timeSeries_config)) + self.assertTrue(myhisto.Configure(histo_config)) + # self.assertTrue(aposterioriPlotter.Configure(aposteriori_config)) + + self.assertTrue(sampler.Run()) + myhisto.data = sampler.data + fitter.data = sampler.data + # lsampler.data = sampler.data + self.assertTrue(myhisto.Run()) + self.assertTrue(fitter.Run()) + # self.assertTrue(lsampler.Run()) + # aposterioriPlotter.data = lsampler.results + # timeSeriesPlotter.data = lsampler.results + # self.assertTrue(timeSeriesPlotter.Run()) + # self.assertTrue(aposterioriPlotter.Run()) + def test_GaussianSampler(self): logger.info("GaussianSampler test") from morpho.processors.sampling.GaussianSamplingProcessor import GaussianSamplingProcessor @@ -182,7 +300,7 @@ def test_linearModelGenerator(self): def test_GaussianRooFit(self): logger.info("GaussianRooFit test") from morpho.processors.sampling import GaussianRooFitProcessor - from morpho.processors.plots import TimeSeries, TimeSeries + from morpho.processors.plots import TimeSeries gaussGen_config = { "iter": 2000, diff --git a/tests/test.sh b/tests/test.sh index 00eb22b1..53347177 100644 --- a/tests/test.sh +++ b/tests/test.sh @@ -1,5 +1,8 @@ #!/bin/bash +# missing scipy (required only for testing) +pip3 install scipy + cd IO && python IO_test.py && cd .. cd misc && python misc_test.py && cd .. cd sampling && python sampling_test.py && cd .. \ No newline at end of file