diff --git a/docs/examples/notebooks/custom_modifiers.ipynb b/docs/examples/notebooks/custom_modifiers.ipynb index 9ba6856ee1..e6c4c993c6 100644 --- a/docs/examples/notebooks/custom_modifiers.ipynb +++ b/docs/examples/notebooks/custom_modifiers.ipynb @@ -1 +1,222 @@ -{"cells":[{"cell_type":"code","execution_count":1,"metadata":{},"outputs":[],"source":["import pyhf\n","import matplotlib.pyplot as plt\n","import numpy as np\n","import scipy.stats"]},{"cell_type":"code","execution_count":5,"metadata":{},"outputs":[],"source":["class custom_modifier_add(object):\n"," op_code = 'addition'\n"," name = 'custom_modifier_add'\n"," def __init__(self,name = 'name'):\n"," self.required_parsets = {\n"," 'mean': [{\n"," 'paramset_type': pyhf.parameters.paramsets.unconstrained,\n"," 'n_parameters': 1,\n"," 'is_constrained': False,\n"," 'is_shared': True,\n"," 'inits': (0.0,),\n"," 'bounds': ((-5, 5),),\n"," 'fixed': False\n"," }]}\n"," def apply(self,pars):\n"," base = np.zeros((1,len(m.config.samples),1,sum(m.config.channel_nbins.values())))\n","\n"," bins = np.linspace(-5,5,20+1)\n"," mean = pars[self.config.par_slice('mean')][0]\n"," yields = 100*(scipy.stats.norm(loc = mean).cdf(bins[1:]) - scipy.stats.norm(loc = mean).cdf(bins[:-1]))\n"," base[0,self.config.samples.index('signal'),0,:] = yields\n"," return base\n","\n","class custom_modifier_mult(object):\n"," op_code = 'multiplication'\n"," name = 'custom_modifier_mult'\n"," def __init__(self,name = 'name'):\n"," self.required_parsets = {\n"," 'mean': [{\n"," 'paramset_type': pyhf.parameters.paramsets.unconstrained,\n"," 'n_parameters': 1,\n"," 'is_constrained': False,\n"," 'is_shared': True,\n"," 'inits': (0.0,),\n"," 'bounds': ((-5, 5),),\n"," 'fixed': False\n"," }]}\n"," def apply(self,pars):\n"," base = np.ones((1,len(m.config.samples),1,sum(m.config.channel_nbins.values())))\n","\n"," bins = np.linspace(-5,5,20+1)\n"," mean = pars[self.config.par_slice('mean')][0]\n"," yields = 100*(scipy.stats.norm(loc = mean).cdf(bins[1:]) - scipy.stats.norm(loc = mean).cdf(bins[:-1]))\n"," base[0,self.config.samples.index('signal'),0,:] = yields\n"," return base\n","\n","\n","\n","m = pyhf.Model(\n","{'channels': [{'name': 'singlechannel',\n"," 'samples': [{'name': 'signal',\n"," 'data': [0]*20,\n"," 'modifiers': [{'name': 'mu', 'type': 'normfactor', 'data': None}]},\n"," {'name': 'background',\n"," 'data': [300]*20,\n"," 'modifiers': []}]}]},\n"," custom_modifiers=[custom_modifier_add()],poi_name = 'mu'\n",")\n","bp = pyhf.tensorlib.astensor(m.config.suggested_init())\n","bp[m.config.poi_index] = 5\n","bp[m.config.par_slice('mean')] = [3.0]\n","d = m.make_pdf(bp).sample()"]},{"cell_type":"code","execution_count":6,"metadata":{},"outputs":[{"output_type":"stream","name":"stdout","text":["[2.96791731]\n[5.44414314]\n"]}],"source":["bestfit = pyhf.infer.mle.fit(d,m)\n","print(bestfit[m.config.par_slice('mean')])\n","print(bestfit[m.config.par_slice('mu')])\n"]},{"cell_type":"code","execution_count":7,"metadata":{},"outputs":[{"output_type":"execute_result","data":{"text/plain":[""]},"metadata":{},"execution_count":7},{"output_type":"display_data","data":{"text/plain":"
","image/svg+xml":"\n\n\n\n \n \n \n \n 2020-11-26T19:25:41.550139\n image/svg+xml\n \n \n Matplotlib v3.3.0, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n","image/png":"iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcUElEQVR4nO3df5BV9Znn8fdDN5GBkIDQiwwNNjFmlCQFmvbHROOyuhN+rAWaOAZjJZhxc8durUpqJj9wrErczVKl0cQddtOkcHElI0YdE5FYRKNGSVIZjeC0BjCO4EDoLn5rkEhkpX32j/ttvDT3dp/b55577r3n86rq6nO+53vuefrc28859zm/zN0REZHGNyLtAEREpDqU8EVEMkIJX0QkI5TwRUQyQglfRCQjmtMOAGDixIne1taWdhgiInVl48aN+929JWr/mkj4bW1tbNiwIe0wRETqipntKKe/SjoiIhmhhC8ikhFK+CIiGVETNfxi3n77bXp6enjrrbfSDiVRo0aNorW1lZEjR6Ydiog0uJpN+D09PYwdO5a2tjbMLO1wEuHuHDhwgJ6eHqZPn552OCLS4Gq2pPPWW28xYcKEhk32AGbGhAkTGv5bjIjUhppN+EBDJ/t+WfgbRaQ21HTCFxGRylHCL8PNN9/M7bffXnL6mjVr2LJlSxUjEhGJTgm/gpTwRZI1e/ZsZs+enXYYdathEv7q1atpa2tjxIgRtLW1sXr16oq87tKlS/nQhz7EhRdeyMsvvwzAnXfeyTnnnMPMmTP59Kc/zeHDh/n1r3/N2rVr+epXv8qsWbPYtm1b0X4iImlpiIS/evVqcrkcO3bswN3ZsWMHuVwudtLfuHEj9913H93d3axbt47nnnsOgE996lM899xzvPDCC5x55pmsXLmSj3/84yxYsIDbbruN7u5uTjvttKL9RETS0hAJ/6abbjph7/nw4cPcdNNNsV73l7/8JZdffjmjR4/mfe97HwsWLABg06ZNfOITn+CjH/0oq1evZvPmzUXnj9pPRKQaavbCq3L8/ve/L6s9rmuuuYY1a9Ywc+ZM7r77bp5++ulY/UREqqEh9vCnTZtWVntUF110EWvWrOFPf/oThw4d4ic/+QkAhw4dYvLkybz99tvHlY3Gjh3LoUOHjo2X6icikoaGSPhLly5l9OjRx7WNHj2apUuXxnrds88+m8985jPMnDmTefPmcc455wDwrW99i/POO48LLriAM84441j/RYsWcdttt3HWWWexbdu2kv1ERFLh7qn/fOxjH/OBtmzZckLbYO655x4/9dRT3cz81FNP9Xvuuaes+dNU7t8qkkUdHR0OOOBNTU3e0dGRdkipAzZ4Gbm2IWr4AFdffTVXX3112mGISAI6OztZvnz5sfG+vr5j411dXWmFVXcil3TMrMnM/tXMHgnj083sWTPbamb3m9l7QvtJYXxrmN6WUOwikoI0Ln5asWJFWe1SXDk1/C8BLxWM3wrc4e4fBF4Hrg3t1wKvh/Y7Qj8RyZjJrdMws8g/k1tLn2TR19dXVrsUF6mkY2atwH8BlgJ/Z/lbPF4MfDZ0WQXcDCwHFoZhgAeB/21mFupNIpIRu3t3curXH4ncf8etlx4bntw6jd29OyPN13/H2VOmTGVXTzKnYjeKqDX8/wl8DRgbxicAf3D3o2G8B5gShqcAOwHc/aiZHQz99xe+oJnlgBzEP31SRBrLwI3F/se6eLN73Qn9xsyaz8Q5ncDxGwwpbsiSjpldCux1942VXLC7r3D3dndvb2lpqeRLi8gg0rwB2e57l7D73iVlzzdxTidjZs1/t8FGHJfsJZooNfwLgAVmth24j3wp5x+BcWbW/w2hFegNw73AVIAw/f3AgbiBllsPjFMv7Lds2TLOPPNMxo8fzy233ALojpgiaZk4p5OTpn6Ek6Z+hFO/tlbJfhiGLOm4+43AjQBmNhv4irtfbWb/DFxBfiOwGHg4zLI2jP9LmP7zStTvy60HDiXK17+uri6eeOIJWltbj7WtWbOGSy+9lBkzZlQsFhFJXv+3muHc4iTOvLUkzpW2Xyd/AHcr+Rp9/60gVwITQvvfAeV/f6sB1113Ha+++irz5s3jjjvu4IYbbih6C2SRLOns7GT9+vWsX7+e5uZmOju1l11Pyrrwyt2fBp4Ow68C5xbp8xbw1xWILVXf//73efTRR3nqqad45JH8N4v+WyBfeumlXHHFFSlHKFJduvip/jXEvXREJHm6+Kn+NcytFUQkWXEvftr/WBdHdm4CYMe3FzBm5tyyD7ye8tlbyuovx9MefpkG3gJZpJ6UW4MvPDtuMEOd/XbCefT+Dm92r2P/YyoFVVPd7OGfMmVqRS+sOGXK1GHNt2jRIr74xS+ybNkyHnzwQU477bSKxSSSpOHU4AvPjotz8dObLzxasr3WT6/s30gCNDc3k8vl6vaYRd0k/DQumd6+fTuQf3LVNddcA8AFF1yg8/ClLg1Wg4+SwPoT87GkbyOil2X8nfLaYxrs1gzFvq2Uui1Dox2orpuELyLxVOIGZBPndHL0QD4xllVPtxHFk7slU1Wu1DeTuBvJWqMavkhGNDU1ldVeSWNmzi2rvZIGKycNpdHu0lnTCT8LN9jMwt8otSGXy5XVXkmp3gsnRjkpzY1kEmo24Y8aNYoDBw40dEJ0dw4cOMCoUaPSDkUyoKuri46OjmPjTU1NdHR0VK00kdq9cEqVjSKUk9LcSCahZmv4ra2t9PT0sG/fvrRDSdSoUaOOu1ePSJK6urqOnXRQ7/eFiWrMzLnFa/gRykn9G8P+A7VNTU06SycJI0eOZPr06WmHISID1NvFT7HOLqKxNpI1m/BFRCpl2GcXNZiareGLiEhlKeGLiGSESjoiGZNmHTrL5ZRaoD18EZGMGHIP38xGAb8ATgr9H3T3b5rZ3cB/BA6Grte4e7flb1Txj8B84HBofz6J4EVEoorz7aLez87pF6WkcwS42N3/aGYjgV+Z2U/DtK+6+4MD+s8DTg8/5wHLw28REUnRkCUdz/tjGB0Zfga7/HUh8IMw3zPAODObHD9UESlX4f3so/6Uuqe91L9IB23NrAnYCHwQ+J67P2tmHcBSM/sG8CSwxN2PAFOAwvuS9oS2XQNeMwfkAKZN0wdMJAmFd42MqpLPnZDaEumgrbv3ufssoBU418w+AtwInAGcA5wMfL2cBbv7Cndvd/f2lpaW8qIWEUlIud+K6ukbUVmnZbr7H8zsKWCuu98emo+Y2f8FvhLGe4HCx0m1hjYRkZpX7reievpGNOQevpm1mNm4MPxnwF8Bv+uvy4ezci4DNoVZ1gKft7zzgYPuvuuEF07Y7NmzmT17drUXKyJSs6Ls4U8GVoU6/gjgAXd/xMx+bmYtgAHdwHWh/zryp2RuJX9a5hcqHnUD699INcppYNWQxXWWxb9Z4hsy4bv7i8BZRdovLtHfgevjhyYiIpWkK21FRDJCCV9EJCOU8EVEMqIhE35nZyfr169n/fr1NDc309lZpWdniojUsIZL+J2dnceePwnQ19fH8uXL6yLpa0NVviyusyz+zVIZDZfwV6xYUVZ7rajnDVVasrjOsvg3S+U0XMLv6+srq72YNC7aqtcNVZqyss4KL/UvTPaFli9fXpeX+mdFrVwI2nBPvGpqaiqa3JuamlKIJrpKbKiyJivrrPBS/8Eu44/SR5I1uXUau3t3lpyevzHB8U6ZMpVdPb9PMqxjGi7h53K5ontBuVwuhWiiS3tDVY9Xbqa9zlJhI8DfKd4uqSt1H57d9y4Bij+EpZob6Ib7lHR1ddHR0XFsvKmpiY6ODrq6ulKMamilNki1vqFKUxbX2ZiZc8tqFynUcHv4kE/6W7ZsAepnj7V/g9T/7aSpqYlcLlfzG6o0ZXGdTZyTPzj7Zve6fIONYMzMucfaRQbTkAk/jv5T3gCam5vLTiBxSiP1uKFKWz2vs+F+VibO6eTogXzNN85zWiV7Gq6kE4dOeRORStv/WBdHdm7iyM5N7Pj2AvY/lt43UCX8Alk5zW8gXciTvGJPUepf5/X+FCUpbf9jXe+W3wD8Hd7sXpda0ldJp0BWTvMrVOpbDdDQtfB+1To7qdjZG7Vy5oYk580XHi3ZnsZxlyhPvBplZr8xsxfMbLOZ/bfQPt3MnjWzrWZ2v5m9J7SfFMa3hultCf8NRT399NNl/xOXOp2vkU/zy+q3GpGqKHYK7WDtCYtS0jkCXOzuM4FZwNzw6MJbgTvc/YPA68C1of+1wOuh/Y7Qry7EPc2vEqWR4Wyo+g3nar5G+FYTZ52lJW5d95TP3qIDtvWg1PURKV03MeRSPe+PYXRk+HHgYuDB0L6K/HNtARaGccL0S6zY5WU1KM45/PV6wDeL32qGo1gNfrCfwWrwtVbXleTU2nUTkWr44Xm2G4EPAt8DtgF/cPejoUsPMCUMTwF2Arj7UTM7CEwA9g94zRyQA5g2bfgHqIa6lLmYwkuZB5u/P2kXJvKo8/YrnH/gJdTlxh512UNdvh1luX19fce9Trnzx1l2peYtZ/6h1lkla/C1VteV5NTadROREr679wGzzGwc8BBwRtwFu/sKYAVAe3u7D/d1Sl3KPJjCf8Y4/8hx73FSbuxJxA0D9jhLfCCHWnacuJOad6j5h/q7EztwWmN1XUlWLV03UdZZOu7+BzN7CvhLYJyZNYe9/FagN3TrBaYCPWbWDLwfOFDBmGtTHd/jpJY+kNVSqqwCJL/3VcefFalvUc7SaQl79pjZnwF/BbwEPAVcEbotBh4Ow2vDOGH6z9192Hvw9SLtWl3aF3fsvnfJsW8Y9WCwskrS0v6sSHZF2cOfDKwKdfwRwAPu/oiZbQHuM7P/AfwrsDL0Xwn8k5ltBV4DFiUQd81Js1aX6t5qvYpZVunfwALs+PaCst7rWqvrSnYMmfDd/UXgrCLtrwLnFml/C/jrikSXkuGWNdIqjegg4DDEKKtUYgObxTJaltXKe6yiYSPQQcCyxSmrpFkOEolDt1ZoBBU4CJjmHshgZxclJVZZRRtYqVNK+A1gzMy5x5cYCtqrIU49O03DLqvoLBupU0r4DUAHjKurUhvYWqnrSnYo4VdYWv/EOmBcPTrLRuqVEr7Ek9F6ts6ykXqkhC/xxKxnp13/V7KWLNFRJoklzumNumukSHVpD7+BpLG3GqeencX6v0ialPAltmHXs+u8/q9ykNQblXQkPTX2NCCRRqf/LEmN7hopUl0q6UhqdD67SHUp4UtF1NsdRkWySCUdEZGMiPLEq6lm9pSZbTGzzWb2pdB+s5n1mll3+JlfMM+NZrbVzF42szlJ/gEiIhJNlJLOUeDv3f15MxsLbDSzx8O0O9z99sLOZjaD/FOuPgz8OfCEmX0oPAhdRERSEuWJV7uAXWH4kJm9BEwZZJaFwH3ufgT49/Cow3OBf6lAvNKAVLsXqY6yavhm1kb+cYfPhqYbzOxFM7vLzMaHtinAzoLZeiiygTCznJltMLMN+/btKz9yEREpS+SEb2bvBX4EfNnd3wCWA6cBs8h/A/hOOQt29xXu3u7u7S0tLeXMKiIiwxAp4ZvZSPLJfrW7/xjA3fe4e5+7vwPcybsPNO8FphbM3hraREQkRVHO0jFgJfCSu3+3oH1yQbfLgU1heC2wyMxOMrPpwOnAbyoXsoiIDEeUs3QuAD4H/NbMukPbPwBXmdkswIHtwN8CuPtmM3sA2EL+DJ/rdYaOiEj6opyl8yvAikw68aGe786zFFgaIy4REakwXWkrIpIRSvgiIhmhhC8ikhFK+CIiGaGELyKSEUr4IiIZoYQvIpIRSvgiIhmhhC8ikhFK+CIiGaGELyKSEUr4IiIZoYQvIpIRSvgiIhmhhC8ikhFRnng11cyeMrMtZrbZzL4U2k82s8fN7JXwe3xoNzNbZmZbwwPOz076jxARkaFF2cM/Cvy9u88AzgeuN7MZwBLgSXc/HXgyjAPMI/9Yw9OBHPmHnYuISMqGTPjuvsvdnw/Dh4CXgCnAQmBV6LYKuCwMLwR+4HnPAOMGPP9WRERSUFYN38zagLOAZ4FJ7r4rTNoNTArDU4CdBbP1hDYREUlR5IRvZu8FfgR82d3fKJzm7k7+YeaRmVnOzDaY2YZ9+/aVM6uIiAxDpIRvZiPJJ/vV7v7j0Lynv1QTfu8N7b3A1ILZW0Pbcdx9hbu3u3t7S0vLcOMXEZGIopylY8BK4CV3/27BpLXA4jC8GHi4oP3z4Wyd84GDBaUfERFJSXOEPhcAnwN+a2bdoe0fgFuAB8zsWmAHcGWYtg6YD2wFDgNfqGTAIiIyPEMmfHf/FWAlJl9SpL8D18eMS0REKkxX2oqIZIQSvohIRijhi4hkhBK+iEhGKOGLiGSEEr6ISEYo4YuIZIQSvohIRijhi4hkhBK+iEhGKOGLiGSEEr6ISEYo4YuIZIQSvohIRijhi4hkhBK+iEhGRHnE4V1mttfMNhW03WxmvWbWHX7mF0y70cy2mtnLZjYnqcBFRKQ8Ufbw7wbmFmm/w91nhZ91AGY2A1gEfDjM02VmTZUKVkREhm/IhO/uvwBei/h6C4H73P2Iu/87+efanhsjPhERqZA4NfwbzOzFUPIZH9qmADsL+vSEthOYWc7MNpjZhn379sUIQ0REohhuwl8OnAbMAnYB3yn3Bdx9hbu3u3t7S0vLMMMQEZGohpXw3X2Pu/e5+zvAnbxbtukFphZ0bQ1tIiKSsmElfDObXDB6OdB/Bs9aYJGZnWRm04HTgd/EC1FERCqheagOZvZDYDYw0cx6gG8Cs81sFuDAduBvAdx9s5k9AGwBjgLXu3tfIpGLiEhZhkz47n5VkeaVg/RfCiyNE5SIiFSerrQVEckIJXwRkYxQwhcRyQglfBGRjFDCFxHJCCV8EZGMUMIXEckIJXwRkYxQwhcRyQglfBGRjFDCFxHJCCV8EZGMUMIXEckIJXwRkYxQwhcRyYghE354SPleM9tU0HaymT1uZq+E3+NDu5nZMjPbGh5wfnaSwYuISHRR9vDvBuYOaFsCPOnupwNPhnGAeeQfa3g6kCP/sHMREakBQyZ8d/8F8NqA5oXAqjC8CrisoP0HnvcMMG7A829FRCQlw63hT3L3XWF4NzApDE8Bdhb06wltJzCznJltMLMN+/btG2YYIiISVeyDtu7u5B9mXu58K9y93d3bW1pa4oYhIiJDGG7C39Nfqgm/94b2XmBqQb/W0CYiIikbbsJfCywOw4uBhwvaPx/O1jkfOFhQ+hERkRQ1D9XBzH4IzAYmmlkP8E3gFuABM7sW2AFcGbqvA+YDW4HDwBcSiFlERIZhyITv7leVmHRJkb4OXB83KBERqTxdaSsikhFK+CIiGaGELyKSEUr4IiIZoYQvIpIRSvgiIhmhhC8ikhFK+CIiGaGELyKSEUr4IiIZoYQvIpIRSvgiIhmhhC8ikhFK+CIiGaGELyKSEUPeD38wZrYdOAT0AUfdvd3MTgbuB9qA7cCV7v56vDBFRCSuSuzh/yd3n+Xu7WF8CfCku58OPBnGRUQkZUmUdBYCq8LwKuCyBJYhIiJlipvwHfiZmW00s1xom1Tw4PLdwKSYyxARkQqIVcMHLnT3XjP7D8DjZva7wonu7mbmxWYMG4gcwLRp02KGISIiQ4m1h+/uveH3XuAh4Fxgj5lNBgi/95aYd4W7t7t7e0tLS5wwREQkgmEnfDMbY2Zj+4eBTwKbgLXA4tBtMfBw3CBFRCS+OCWdScBDZtb/Ove6+6Nm9hzwgJldC+wArowfpoiIxDXshO/urwIzi7QfAC6JE5SIiFSerrQVEckIJXwRkYxQwhcRyQglfBGRjFDCFxHJCCV8EZGMUMIXEckIJXwRkYxQwhcRyQglfBGRjFDCFxHJCCV8EZGMUMIXEckIJXwRkYxQwhcRyQglfBGRjEgs4ZvZXDN72cy2mtmSpJYjIiLRJJLwzawJ+B4wD5gBXGVmM5JYloiIRJPUHv65wFZ3f9Xd/x9wH7AwoWWJiEgE5u6Vf1GzK4C57v5fw/jngPPc/YaCPjkgF0b/Ani54oHARGB/Aq8bl+IqX63GVqtxQe3GVqtxQe3GViquU929JeqLDPsh5nG5+wpgRZLLMLMN7t6e5DKGQ3GVr1Zjq9W4oHZjq9W4oHZjq1RcSZV0eoGpBeOtoU1ERFKSVMJ/DjjdzKab2XuARcDahJYlIiIRJFLScfejZnYD8BjQBNzl7puTWNYQEi0ZxaC4ylersdVqXFC7sdVqXFC7sVUkrkQO2oqISO3RlbYiIhmhhC8ikhF1n/CHuoWDmZ1kZveH6c+aWVuV4ppqZk+Z2RYz22xmXyrSZ7aZHTSz7vDzjSrFtt3MfhuWuaHIdDOzZWGdvWhmZ1cprr8oWBfdZvaGmX15QJ+qrDMzu8vM9prZpoK2k83scTN7JfweX2LexaHPK2a2uEqx3WZmvwvv10NmNq7EvIO+9wnEdbOZ9Ra8X/NLzJvorVhKxHZ/QVzbzay7xLxJrrOieSKxz5q71+0P+QPC24APAO8BXgBmDOjTCXw/DC8C7q9SbJOBs8PwWODfisQ2G3gkhfW2HZg4yPT5wE8BA84Hnk3pvd1N/sKSqq8z4CLgbGBTQdu3gSVheAlwa5H5TgZeDb/Hh+HxVYjtk0BzGL61WGxR3vsE4roZ+EqE93rQ/+MkYhsw/TvAN1JYZ0XzRFKftXrfw49yC4eFwKow/CBwiZlZ0oG5+y53fz4MHwJeAqYkvdwKWQj8wPOeAcaZ2eQqx3AJsM3dd1R5uQC4+y+A1wY0F36WVgGXFZl1DvC4u7/m7q8DjwNzk47N3X/m7kfD6DPkr32pqhLrLIrEb8UyWGwhH1wJ/LCSy4xikDyRyGet3hP+FGBnwXgPJybVY33CP8RBYEJVogtCGeks4Nkik//SzF4ws5+a2YerFJIDPzOzjZa/xcVAUdZr0hZR+h8wjXUGMMndd4Xh3cCkIn1qYd39DflvaMUM9d4n4YZQarqrRGki7XX2CWCPu79SYnpV1tmAPJHIZ63eE37NM7P3Aj8CvuzubwyY/Dz5ksVM4H8Ba6oU1oXufjb5u5leb2YXVWm5kVj+Yr0FwD8XmZzWOjuO579T19w5zWZ2E3AUWF2iS7Xf++XAacAsYBf50kmtuYrB9+4TX2eD5YlKftbqPeFHuYXDsT5m1gy8HzhQjeDMbCT5N3G1u/944HR3f8Pd/xiG1wEjzWxi0nG5e2/4vRd4iPxX6kJp3xpjHvC8u+8ZOCGtdRbs6S9thd97i/RJbd2Z2TXApcDVIUmcIMJ7X1Huvsfd+9z9HeDOEstLc501A58C7i/VJ+l1ViJPJPJZq/eEH+UWDmuB/qPXVwA/L/XPUEmhLrgSeMndv1uizyn9xxPM7Fzy70eiGyMzG2NmY/uHyR/s2zSg21rg85Z3PnCw4OtlNZTc40pjnRUo/CwtBh4u0ucx4JNmNj6ULz4Z2hJlZnOBrwEL3P1wiT5R3vtKx1V47OfyEstL81Ys/xn4nbv3FJuY9DobJE8k81lL4shzNX/In1Hyb+SP8t8U2v47+Q8+wCjypYGtwG+AD1QprgvJfw17EegOP/OB64DrQp8bgM3kz0p4Bvh4FeL6QFjeC2HZ/eusMC4j/wCbbcBvgfYqvp9jyCfw9xe0VX2dkd/g7ALeJl8bvZb8sZ8ngVeAJ4CTQ9924P8UzPs34fO2FfhClWLbSr6e2/9Z6z8z7c+BdYO99wnH9U/hM/Qi+SQ2eWBcYfyE/+OkYwvtd/d/tgr6VnOdlcoTiXzWdGsFEZGMqPeSjoiIRKSELyKSEUr4IiIZoYQvIpIRSvgiIhmhhC8ikhFK+CIiGfH/AUsFfjbA7joZAAAAAElFTkSuQmCC\n"},"metadata":{"needs_background":"light"}}],"source":["plt.bar(np.arange(20),m.expected_actualdata(bestfit), alpha = 1.0, facecolor = None, edgecolor = 'k', label = 'fit')\n","plt.scatter(np.arange(20),d[:m.config.nmaindata], alpha = 1.0, marker = 'o', c='k', label = 'data', zorder=99)\n","plt.errorbar(np.arange(20),d[:m.config.nmaindata],yerr=np.sqrt(d[:m.config.nmaindata]), marker='o', c='k', linestyle = '')\n","plt.legend()\n"]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["s"]}],"nbformat":4,"nbformat_minor":2,"metadata":{"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.2-final"},"orig_nbformat":2,"kernelspec":{"name":"python3","display_name":"Python 3"}}} \ No newline at end of file +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import pyhf\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import scipy.stats" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "class custom_modifier_add:\n", + " op_code = 'addition'\n", + " name = 'custom_modifier_add'\n", + "\n", + " def __init__(self, name='name'):\n", + " self.required_parsets = {\n", + " 'mean': [\n", + " {\n", + " 'paramset_type': pyhf.parameters.paramsets.unconstrained,\n", + " 'n_parameters': 1,\n", + " 'is_constrained': False,\n", + " 'is_shared': True,\n", + " 'inits': (0.0,),\n", + " 'bounds': ((-5, 5),),\n", + " 'fixed': False,\n", + " }\n", + " ]\n", + " }\n", + "\n", + " def apply(self, pars):\n", + " base = np.zeros(\n", + " (1, len(m.config.samples), 1, sum(m.config.channel_nbins.values()))\n", + " )\n", + "\n", + " bins = np.linspace(-5, 5, 20 + 1)\n", + " mean = pars[self.config.par_slice('mean')][0]\n", + " yields = 100 * (\n", + " scipy.stats.norm(loc=mean).cdf(bins[1:])\n", + " - scipy.stats.norm(loc=mean).cdf(bins[:-1])\n", + " )\n", + " base[0, self.config.samples.index('signal'), 0, :] = yields\n", + " return base\n", + "\n", + "\n", + "class custom_modifier_mult:\n", + " op_code = 'multiplication'\n", + " name = 'custom_modifier_mult'\n", + "\n", + " def __init__(self, name='name'):\n", + " self.required_parsets = {\n", + " 'mean': [\n", + " {\n", + " 'paramset_type': pyhf.parameters.paramsets.unconstrained,\n", + " 'n_parameters': 1,\n", + " 'is_constrained': False,\n", + " 'is_shared': True,\n", + " 'inits': (0.0,),\n", + " 'bounds': ((-5, 5),),\n", + " 'fixed': False,\n", + " }\n", + " ]\n", + " }\n", + "\n", + " def apply(self, pars):\n", + " base = np.ones(\n", + " (1, len(m.config.samples), 1, sum(m.config.channel_nbins.values()))\n", + " )\n", + "\n", + " bins = np.linspace(-5, 5, 20 + 1)\n", + " mean = pars[self.config.par_slice('mean')][0]\n", + " yields = 100 * (\n", + " scipy.stats.norm(loc=mean).cdf(bins[1:])\n", + " - scipy.stats.norm(loc=mean).cdf(bins[:-1])\n", + " )\n", + " base[0, self.config.samples.index('signal'), 0, :] = yields\n", + " return base\n", + "\n", + "\n", + "m = pyhf.Model(\n", + " {\n", + " 'channels': [\n", + " {\n", + " 'name': 'singlechannel',\n", + " 'samples': [\n", + " {\n", + " 'name': 'signal',\n", + " 'data': [0] * 20,\n", + " 'modifiers': [\n", + " {'name': 'mu', 'type': 'normfactor', 'data': None}\n", + " ],\n", + " },\n", + " {'name': 'background', 'data': [300] * 20, 'modifiers': []},\n", + " ],\n", + " }\n", + " ]\n", + " },\n", + " custom_modifiers=[custom_modifier_add()],\n", + " poi_name='mu',\n", + ")\n", + "bp = pyhf.tensorlib.astensor(m.config.suggested_init())\n", + "bp[m.config.poi_index] = 5\n", + "bp[m.config.par_slice('mean')] = [3.0]\n", + "d = m.make_pdf(bp).sample()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "[2.96791731]\n[5.44414314]\n" + ] + } + ], + "source": [ + "bestfit = pyhf.infer.mle.fit(d, m)\n", + "print(bestfit[m.config.par_slice('mean')])\n", + "print(bestfit[m.config.par_slice('mu')])" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "" + ] + }, + "metadata": {}, + "execution_count": 7 + }, + { + "output_type": "display_data", + "data": { + "text/plain": "
", + "image/svg+xml": "\n\n\n\n \n \n \n \n 2020-11-26T19:25:41.550139\n image/svg+xml\n \n \n Matplotlib v3.3.0, https://matplotlib.org/\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcUElEQVR4nO3df5BV9Znn8fdDN5GBkIDQiwwNNjFmlCQFmvbHROOyuhN+rAWaOAZjJZhxc8durUpqJj9wrErczVKl0cQddtOkcHElI0YdE5FYRKNGSVIZjeC0BjCO4EDoLn5rkEhkpX32j/ttvDT3dp/b55577r3n86rq6nO+53vuefrc28859zm/zN0REZHGNyLtAEREpDqU8EVEMkIJX0QkI5TwRUQyQglfRCQjmtMOAGDixIne1taWdhgiInVl48aN+929JWr/mkj4bW1tbNiwIe0wRETqipntKKe/SjoiIhmhhC8ikhFK+CIiGVETNfxi3n77bXp6enjrrbfSDiVRo0aNorW1lZEjR6Ydiog0uJpN+D09PYwdO5a2tjbMLO1wEuHuHDhwgJ6eHqZPn552OCLS4Gq2pPPWW28xYcKEhk32AGbGhAkTGv5bjIjUhppN+EBDJ/t+WfgbRaQ21HTCFxGRylHCL8PNN9/M7bffXnL6mjVr2LJlSxUjEhGJTgm/gpTwRZI1e/ZsZs+enXYYdathEv7q1atpa2tjxIgRtLW1sXr16oq87tKlS/nQhz7EhRdeyMsvvwzAnXfeyTnnnMPMmTP59Kc/zeHDh/n1r3/N2rVr+epXv8qsWbPYtm1b0X4iImlpiIS/evVqcrkcO3bswN3ZsWMHuVwudtLfuHEj9913H93d3axbt47nnnsOgE996lM899xzvPDCC5x55pmsXLmSj3/84yxYsIDbbruN7u5uTjvttKL9RETS0hAJ/6abbjph7/nw4cPcdNNNsV73l7/8JZdffjmjR4/mfe97HwsWLABg06ZNfOITn+CjH/0oq1evZvPmzUXnj9pPRKQaavbCq3L8/ve/L6s9rmuuuYY1a9Ywc+ZM7r77bp5++ulY/UREqqEh9vCnTZtWVntUF110EWvWrOFPf/oThw4d4ic/+QkAhw4dYvLkybz99tvHlY3Gjh3LoUOHjo2X6icikoaGSPhLly5l9OjRx7WNHj2apUuXxnrds88+m8985jPMnDmTefPmcc455wDwrW99i/POO48LLriAM84441j/RYsWcdttt3HWWWexbdu2kv1ERFLh7qn/fOxjH/OBtmzZckLbYO655x4/9dRT3cz81FNP9Xvuuaes+dNU7t8qkkUdHR0OOOBNTU3e0dGRdkipAzZ4Gbm2IWr4AFdffTVXX3112mGISAI6OztZvnz5sfG+vr5j411dXWmFVXcil3TMrMnM/tXMHgnj083sWTPbamb3m9l7QvtJYXxrmN6WUOwikoI0Ln5asWJFWe1SXDk1/C8BLxWM3wrc4e4fBF4Hrg3t1wKvh/Y7Qj8RyZjJrdMws8g/k1tLn2TR19dXVrsUF6mkY2atwH8BlgJ/Z/lbPF4MfDZ0WQXcDCwHFoZhgAeB/21mFupNIpIRu3t3curXH4ncf8etlx4bntw6jd29OyPN13/H2VOmTGVXTzKnYjeKqDX8/wl8DRgbxicAf3D3o2G8B5gShqcAOwHc/aiZHQz99xe+oJnlgBzEP31SRBrLwI3F/se6eLN73Qn9xsyaz8Q5ncDxGwwpbsiSjpldCux1942VXLC7r3D3dndvb2lpqeRLi8gg0rwB2e57l7D73iVlzzdxTidjZs1/t8FGHJfsJZooNfwLgAVmth24j3wp5x+BcWbW/w2hFegNw73AVIAw/f3AgbiBllsPjFMv7Lds2TLOPPNMxo8fzy233ALojpgiaZk4p5OTpn6Ek6Z+hFO/tlbJfhiGLOm4+43AjQBmNhv4irtfbWb/DFxBfiOwGHg4zLI2jP9LmP7zStTvy60HDiXK17+uri6eeOIJWltbj7WtWbOGSy+9lBkzZlQsFhFJXv+3muHc4iTOvLUkzpW2Xyd/AHcr+Rp9/60gVwITQvvfAeV/f6sB1113Ha+++irz5s3jjjvu4IYbbih6C2SRLOns7GT9+vWsX7+e5uZmOju1l11Pyrrwyt2fBp4Ow68C5xbp8xbw1xWILVXf//73efTRR3nqqad45JH8N4v+WyBfeumlXHHFFSlHKFJduvip/jXEvXREJHm6+Kn+NcytFUQkWXEvftr/WBdHdm4CYMe3FzBm5tyyD7ye8tlbyuovx9MefpkG3gJZpJ6UW4MvPDtuMEOd/XbCefT+Dm92r2P/YyoFVVPd7OGfMmVqRS+sOGXK1GHNt2jRIr74xS+ybNkyHnzwQU477bSKxSSSpOHU4AvPjotz8dObLzxasr3WT6/s30gCNDc3k8vl6vaYRd0k/DQumd6+fTuQf3LVNddcA8AFF1yg8/ClLg1Wg4+SwPoT87GkbyOil2X8nfLaYxrs1gzFvq2Uui1Dox2orpuELyLxVOIGZBPndHL0QD4xllVPtxHFk7slU1Wu1DeTuBvJWqMavkhGNDU1ldVeSWNmzi2rvZIGKycNpdHu0lnTCT8LN9jMwt8otSGXy5XVXkmp3gsnRjkpzY1kEmo24Y8aNYoDBw40dEJ0dw4cOMCoUaPSDkUyoKuri46OjmPjTU1NdHR0VK00kdq9cEqVjSKUk9LcSCahZmv4ra2t9PT0sG/fvrRDSdSoUaOOu1ePSJK6urqOnXRQ7/eFiWrMzLnFa/gRykn9G8P+A7VNTU06SycJI0eOZPr06WmHISID1NvFT7HOLqKxNpI1m/BFRCpl2GcXNZiareGLiEhlKeGLiGSESjoiGZNmHTrL5ZRaoD18EZGMGHIP38xGAb8ATgr9H3T3b5rZ3cB/BA6Grte4e7flb1Txj8B84HBofz6J4EVEoorz7aLez87pF6WkcwS42N3/aGYjgV+Z2U/DtK+6+4MD+s8DTg8/5wHLw28REUnRkCUdz/tjGB0Zfga7/HUh8IMw3zPAODObHD9UESlX4f3so/6Uuqe91L9IB23NrAnYCHwQ+J67P2tmHcBSM/sG8CSwxN2PAFOAwvuS9oS2XQNeMwfkAKZN0wdMJAmFd42MqpLPnZDaEumgrbv3ufssoBU418w+AtwInAGcA5wMfL2cBbv7Cndvd/f2lpaW8qIWEUlIud+K6ukbUVmnZbr7H8zsKWCuu98emo+Y2f8FvhLGe4HCx0m1hjYRkZpX7reievpGNOQevpm1mNm4MPxnwF8Bv+uvy4ezci4DNoVZ1gKft7zzgYPuvuuEF07Y7NmzmT17drUXKyJSs6Ls4U8GVoU6/gjgAXd/xMx+bmYtgAHdwHWh/zryp2RuJX9a5hcqHnUD699INcppYNWQxXWWxb9Z4hsy4bv7i8BZRdovLtHfgevjhyYiIpWkK21FRDJCCV9EJCOU8EVEMqIhE35nZyfr169n/fr1NDc309lZpWdniojUsIZL+J2dnceePwnQ19fH8uXL6yLpa0NVviyusyz+zVIZDZfwV6xYUVZ7rajnDVVasrjOsvg3S+U0XMLv6+srq72YNC7aqtcNVZqyss4KL/UvTPaFli9fXpeX+mdFrVwI2nBPvGpqaiqa3JuamlKIJrpKbKiyJivrrPBS/8Eu44/SR5I1uXUau3t3lpyevzHB8U6ZMpVdPb9PMqxjGi7h53K5ontBuVwuhWiiS3tDVY9Xbqa9zlJhI8DfKd4uqSt1H57d9y4Bij+EpZob6Ib7lHR1ddHR0XFsvKmpiY6ODrq6ulKMamilNki1vqFKUxbX2ZiZc8tqFynUcHv4kE/6W7ZsAepnj7V/g9T/7aSpqYlcLlfzG6o0ZXGdTZyTPzj7Zve6fIONYMzMucfaRQbTkAk/jv5T3gCam5vLTiBxSiP1uKFKWz2vs+F+VibO6eTogXzNN85zWiV7Gq6kE4dOeRORStv/WBdHdm7iyM5N7Pj2AvY/lt43UCX8Alk5zW8gXciTvGJPUepf5/X+FCUpbf9jXe+W3wD8Hd7sXpda0ldJp0BWTvMrVOpbDdDQtfB+1To7qdjZG7Vy5oYk580XHi3ZnsZxlyhPvBplZr8xsxfMbLOZ/bfQPt3MnjWzrWZ2v5m9J7SfFMa3hultCf8NRT399NNl/xOXOp2vkU/zy+q3GpGqKHYK7WDtCYtS0jkCXOzuM4FZwNzw6MJbgTvc/YPA68C1of+1wOuh/Y7Qry7EPc2vEqWR4Wyo+g3nar5G+FYTZ52lJW5d95TP3qIDtvWg1PURKV03MeRSPe+PYXRk+HHgYuDB0L6K/HNtARaGccL0S6zY5WU1KM45/PV6wDeL32qGo1gNfrCfwWrwtVbXleTU2nUTkWr44Xm2G4EPAt8DtgF/cPejoUsPMCUMTwF2Arj7UTM7CEwA9g94zRyQA5g2bfgHqIa6lLmYwkuZB5u/P2kXJvKo8/YrnH/gJdTlxh512UNdvh1luX19fce9Trnzx1l2peYtZ/6h1lkla/C1VteV5NTadROREr679wGzzGwc8BBwRtwFu/sKYAVAe3u7D/d1Sl3KPJjCf8Y4/8hx73FSbuxJxA0D9jhLfCCHWnacuJOad6j5h/q7EztwWmN1XUlWLV03UdZZOu7+BzN7CvhLYJyZNYe9/FagN3TrBaYCPWbWDLwfOFDBmGtTHd/jpJY+kNVSqqwCJL/3VcefFalvUc7SaQl79pjZnwF/BbwEPAVcEbotBh4Ow2vDOGH6z9192Hvw9SLtWl3aF3fsvnfJsW8Y9WCwskrS0v6sSHZF2cOfDKwKdfwRwAPu/oiZbQHuM7P/AfwrsDL0Xwn8k5ltBV4DFiUQd81Js1aX6t5qvYpZVunfwALs+PaCst7rWqvrSnYMmfDd/UXgrCLtrwLnFml/C/jrikSXkuGWNdIqjegg4DDEKKtUYgObxTJaltXKe6yiYSPQQcCyxSmrpFkOEolDt1ZoBBU4CJjmHshgZxclJVZZRRtYqVNK+A1gzMy5x5cYCtqrIU49O03DLqvoLBupU0r4DUAHjKurUhvYWqnrSnYo4VdYWv/EOmBcPTrLRuqVEr7Ek9F6ts6ykXqkhC/xxKxnp13/V7KWLNFRJoklzumNumukSHVpD7+BpLG3GqeencX6v0ialPAltmHXs+u8/q9ykNQblXQkPTX2NCCRRqf/LEmN7hopUl0q6UhqdD67SHUp4UtF1NsdRkWySCUdEZGMiPLEq6lm9pSZbTGzzWb2pdB+s5n1mll3+JlfMM+NZrbVzF42szlJ/gEiIhJNlJLOUeDv3f15MxsLbDSzx8O0O9z99sLOZjaD/FOuPgz8OfCEmX0oPAhdRERSEuWJV7uAXWH4kJm9BEwZZJaFwH3ufgT49/Cow3OBf6lAvNKAVLsXqY6yavhm1kb+cYfPhqYbzOxFM7vLzMaHtinAzoLZeiiygTCznJltMLMN+/btKz9yEREpS+SEb2bvBX4EfNnd3wCWA6cBs8h/A/hOOQt29xXu3u7u7S0tLeXMKiIiwxAp4ZvZSPLJfrW7/xjA3fe4e5+7vwPcybsPNO8FphbM3hraREQkRVHO0jFgJfCSu3+3oH1yQbfLgU1heC2wyMxOMrPpwOnAbyoXsoiIDEeUs3QuAD4H/NbMukPbPwBXmdkswIHtwN8CuPtmM3sA2EL+DJ/rdYaOiEj6opyl8yvAikw68aGe786zFFgaIy4REakwXWkrIpIRSvgiIhmhhC8ikhFK+CIiGaGELyKSEUr4IiIZoYQvIpIRSvgiIhmhhC8ikhFK+CIiGaGELyKSEUr4IiIZoYQvIpIRSvgiIhmhhC8ikhFRnng11cyeMrMtZrbZzL4U2k82s8fN7JXwe3xoNzNbZmZbwwPOz076jxARkaFF2cM/Cvy9u88AzgeuN7MZwBLgSXc/HXgyjAPMI/9Yw9OBHPmHnYuISMqGTPjuvsvdnw/Dh4CXgCnAQmBV6LYKuCwMLwR+4HnPAOMGPP9WRERSUFYN38zagLOAZ4FJ7r4rTNoNTArDU4CdBbP1hDYREUlR5IRvZu8FfgR82d3fKJzm7k7+YeaRmVnOzDaY2YZ9+/aVM6uIiAxDpIRvZiPJJ/vV7v7j0Lynv1QTfu8N7b3A1ILZW0Pbcdx9hbu3u3t7S0vLcOMXEZGIopylY8BK4CV3/27BpLXA4jC8GHi4oP3z4Wyd84GDBaUfERFJSXOEPhcAnwN+a2bdoe0fgFuAB8zsWmAHcGWYtg6YD2wFDgNfqGTAIiIyPEMmfHf/FWAlJl9SpL8D18eMS0REKkxX2oqIZIQSvohIRijhi4hkhBK+iEhGKOGLiGSEEr6ISEYo4YuIZIQSvohIRijhi4hkhBK+iEhGKOGLiGSEEr6ISEYo4YuIZIQSvohIRijhi4hkhBK+iEhGRHnE4V1mttfMNhW03WxmvWbWHX7mF0y70cy2mtnLZjYnqcBFRKQ8Ufbw7wbmFmm/w91nhZ91AGY2A1gEfDjM02VmTZUKVkREhm/IhO/uvwBei/h6C4H73P2Iu/87+efanhsjPhERqZA4NfwbzOzFUPIZH9qmADsL+vSEthOYWc7MNpjZhn379sUIQ0REohhuwl8OnAbMAnYB3yn3Bdx9hbu3u3t7S0vLMMMQEZGohpXw3X2Pu/e5+zvAnbxbtukFphZ0bQ1tIiKSsmElfDObXDB6OdB/Bs9aYJGZnWRm04HTgd/EC1FERCqheagOZvZDYDYw0cx6gG8Cs81sFuDAduBvAdx9s5k9AGwBjgLXu3tfIpGLiEhZhkz47n5VkeaVg/RfCiyNE5SIiFSerrQVEckIJXwRkYxQwhcRyQglfBGRjFDCFxHJCCV8EZGMUMIXEckIJXwRkYxQwhcRyQglfBGRjFDCFxHJCCV8EZGMUMIXEckIJXwRkYxQwhcRyYghE354SPleM9tU0HaymT1uZq+E3+NDu5nZMjPbGh5wfnaSwYuISHRR9vDvBuYOaFsCPOnupwNPhnGAeeQfa3g6kCP/sHMREakBQyZ8d/8F8NqA5oXAqjC8CrisoP0HnvcMMG7A829FRCQlw63hT3L3XWF4NzApDE8Bdhb06wltJzCznJltMLMN+/btG2YYIiISVeyDtu7u5B9mXu58K9y93d3bW1pa4oYhIiJDGG7C39Nfqgm/94b2XmBqQb/W0CYiIikbbsJfCywOw4uBhwvaPx/O1jkfOFhQ+hERkRQ1D9XBzH4IzAYmmlkP8E3gFuABM7sW2AFcGbqvA+YDW4HDwBcSiFlERIZhyITv7leVmHRJkb4OXB83KBERqTxdaSsikhFK+CIiGaGELyKSEUr4IiIZoYQvIpIRSvgiIhmhhC8ikhFK+CIiGaGELyKSEUr4IiIZoYQvIpIRSvgiIhmhhC8ikhFK+CIiGaGELyKSEUPeD38wZrYdOAT0AUfdvd3MTgbuB9qA7cCV7v56vDBFRCSuSuzh/yd3n+Xu7WF8CfCku58OPBnGRUQkZUmUdBYCq8LwKuCyBJYhIiJlipvwHfiZmW00s1xom1Tw4PLdwKSYyxARkQqIVcMHLnT3XjP7D8DjZva7wonu7mbmxWYMG4gcwLRp02KGISIiQ4m1h+/uveH3XuAh4Fxgj5lNBgi/95aYd4W7t7t7e0tLS5wwREQkgmEnfDMbY2Zj+4eBTwKbgLXA4tBtMfBw3CBFRCS+OCWdScBDZtb/Ove6+6Nm9hzwgJldC+wArowfpoiIxDXshO/urwIzi7QfAC6JE5SIiFSerrQVEckIJXwRkYxQwhcRyQglfBGRjFDCFxHJCCV8EZGMUMIXEckIJXwRkYxQwhcRyQglfBGRjFDCFxHJCCV8EZGMUMIXEckIJXwRkYxQwhcRyQglfBGRjEgs4ZvZXDN72cy2mtmSpJYjIiLRJJLwzawJ+B4wD5gBXGVmM5JYloiIRJPUHv65wFZ3f9Xd/x9wH7AwoWWJiEgE5u6Vf1GzK4C57v5fw/jngPPc/YaCPjkgF0b/Ani54oHARGB/Aq8bl+IqX63GVqtxQe3GVqtxQe3GViquU929JeqLDPsh5nG5+wpgRZLLMLMN7t6e5DKGQ3GVr1Zjq9W4oHZjq9W4oHZjq1RcSZV0eoGpBeOtoU1ERFKSVMJ/DjjdzKab2XuARcDahJYlIiIRJFLScfejZnYD8BjQBNzl7puTWNYQEi0ZxaC4ylersdVqXFC7sdVqXFC7sVUkrkQO2oqISO3RlbYiIhmhhC8ikhF1n/CHuoWDmZ1kZveH6c+aWVuV4ppqZk+Z2RYz22xmXyrSZ7aZHTSz7vDzjSrFtt3MfhuWuaHIdDOzZWGdvWhmZ1cprr8oWBfdZvaGmX15QJ+qrDMzu8vM9prZpoK2k83scTN7JfweX2LexaHPK2a2uEqx3WZmvwvv10NmNq7EvIO+9wnEdbOZ9Ra8X/NLzJvorVhKxHZ/QVzbzay7xLxJrrOieSKxz5q71+0P+QPC24APAO8BXgBmDOjTCXw/DC8C7q9SbJOBs8PwWODfisQ2G3gkhfW2HZg4yPT5wE8BA84Hnk3pvd1N/sKSqq8z4CLgbGBTQdu3gSVheAlwa5H5TgZeDb/Hh+HxVYjtk0BzGL61WGxR3vsE4roZ+EqE93rQ/+MkYhsw/TvAN1JYZ0XzRFKftXrfw49yC4eFwKow/CBwiZlZ0oG5+y53fz4MHwJeAqYkvdwKWQj8wPOeAcaZ2eQqx3AJsM3dd1R5uQC4+y+A1wY0F36WVgGXFZl1DvC4u7/m7q8DjwNzk47N3X/m7kfD6DPkr32pqhLrLIrEb8UyWGwhH1wJ/LCSy4xikDyRyGet3hP+FGBnwXgPJybVY33CP8RBYEJVogtCGeks4Nkik//SzF4ws5+a2YerFJIDPzOzjZa/xcVAUdZr0hZR+h8wjXUGMMndd4Xh3cCkIn1qYd39DflvaMUM9d4n4YZQarqrRGki7XX2CWCPu79SYnpV1tmAPJHIZ63eE37NM7P3Aj8CvuzubwyY/Dz5ksVM4H8Ba6oU1oXufjb5u5leb2YXVWm5kVj+Yr0FwD8XmZzWOjuO579T19w5zWZ2E3AUWF2iS7Xf++XAacAsYBf50kmtuYrB9+4TX2eD5YlKftbqPeFHuYXDsT5m1gy8HzhQjeDMbCT5N3G1u/944HR3f8Pd/xiG1wEjzWxi0nG5e2/4vRd4iPxX6kJp3xpjHvC8u+8ZOCGtdRbs6S9thd97i/RJbd2Z2TXApcDVIUmcIMJ7X1Huvsfd+9z9HeDOEstLc501A58C7i/VJ+l1ViJPJPJZq/eEH+UWDmuB/qPXVwA/L/XPUEmhLrgSeMndv1uizyn9xxPM7Fzy70eiGyMzG2NmY/uHyR/s2zSg21rg85Z3PnCw4OtlNZTc40pjnRUo/CwtBh4u0ucx4JNmNj6ULz4Z2hJlZnOBrwEL3P1wiT5R3vtKx1V47OfyEstL81Ys/xn4nbv3FJuY9DobJE8k81lL4shzNX/In1Hyb+SP8t8U2v47+Q8+wCjypYGtwG+AD1QprgvJfw17EegOP/OB64DrQp8bgM3kz0p4Bvh4FeL6QFjeC2HZ/eusMC4j/wCbbcBvgfYqvp9jyCfw9xe0VX2dkd/g7ALeJl8bvZb8sZ8ngVeAJ4CTQ9924P8UzPs34fO2FfhClWLbSr6e2/9Z6z8z7c+BdYO99wnH9U/hM/Qi+SQ2eWBcYfyE/+OkYwvtd/d/tgr6VnOdlcoTiXzWdGsFEZGMqPeSjoiIRKSELyKSEUr4IiIZoYQvIpIRSvgiIhmhhC8ikhFK+CIiGfH/AUsFfjbA7joZAAAAAElFTkSuQmCC\n" + }, + "metadata": { + "needs_background": "light" + } + } + ], + "source": [ + "plt.bar(\n", + " np.arange(20),\n", + " m.expected_actualdata(bestfit),\n", + " alpha=1.0,\n", + " facecolor=None,\n", + " edgecolor='k',\n", + " label='fit',\n", + ")\n", + "plt.scatter(\n", + " np.arange(20),\n", + " d[: m.config.nmaindata],\n", + " alpha=1.0,\n", + " marker='o',\n", + " c='k',\n", + " label='data',\n", + " zorder=99,\n", + ")\n", + "plt.errorbar(\n", + " np.arange(20),\n", + " d[: m.config.nmaindata],\n", + " yerr=np.sqrt(d[: m.config.nmaindata]),\n", + " marker='o',\n", + " c='k',\n", + " linestyle='',\n", + ")\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s" + ] + } + ], + "nbformat": 4, + "nbformat_minor": 2, + "metadata": { + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.2-final" + }, + "orig_nbformat": 2, + "kernelspec": { + "name": "python3", + "display_name": "Python 3" + } + } +} diff --git a/src/pyhf/modifiers/histosys.py b/src/pyhf/modifiers/histosys.py index c266780f36..1aea66b6b6 100644 --- a/src/pyhf/modifiers/histosys.py +++ b/src/pyhf/modifiers/histosys.py @@ -29,8 +29,8 @@ class histosys_combined: def __init__( self, histosys_mods, pdfconfig, mega_mods, interpcode='code0', batch_size=None ): - self.name = histosys.name - self.op_code = histosys.op_code + self.name = histosys.name + self.op_code = histosys.op_code self.batch_size = batch_size self.interpcode = interpcode assert self.interpcode in ['code0', 'code2', 'code4p'] diff --git a/src/pyhf/modifiers/lumi.py b/src/pyhf/modifiers/lumi.py index 8d5cfbdcbc..7963f7b5a0 100644 --- a/src/pyhf/modifiers/lumi.py +++ b/src/pyhf/modifiers/lumi.py @@ -29,8 +29,8 @@ def required_parset(cls, sample_data, modifier_data): class lumi_combined: def __init__(self, lumi_mods, pdfconfig, mega_mods, batch_size=None): self.batch_size = batch_size - self.name = lumi.name - self.op_code = lumi.op_code + self.name = lumi.name + self.op_code = lumi.op_code keys = [f'{mtype}/{m}' for m, mtype in lumi_mods] lumi_mods = [m for m, _ in lumi_mods] diff --git a/src/pyhf/modifiers/normfactor.py b/src/pyhf/modifiers/normfactor.py index 8aaf1331d5..0a61b3203c 100644 --- a/src/pyhf/modifiers/normfactor.py +++ b/src/pyhf/modifiers/normfactor.py @@ -26,8 +26,8 @@ def required_parset(cls, sample_data, modifier_data): class normfactor_combined: def __init__(self, normfactor_mods, pdfconfig, mega_mods, batch_size=None): self.batch_size = batch_size - self.name = normfactor.name - self.op_code = normfactor.op_code + self.name = normfactor.name + self.op_code = normfactor.op_code keys = [f'{mtype}/{m}' for m, mtype in normfactor_mods] normfactor_mods = [m for m, _ in normfactor_mods] diff --git a/src/pyhf/modifiers/normsys.py b/src/pyhf/modifiers/normsys.py index 45b93979a3..11ac7136f3 100644 --- a/src/pyhf/modifiers/normsys.py +++ b/src/pyhf/modifiers/normsys.py @@ -29,8 +29,8 @@ class normsys_combined: def __init__( self, normsys_mods, pdfconfig, mega_mods, interpcode='code1', batch_size=None ): - self.name = normsys.name - self.op_code = normsys.op_code + self.name = normsys.name + self.op_code = normsys.op_code self.interpcode = interpcode assert self.interpcode in ['code1', 'code4'] diff --git a/src/pyhf/modifiers/shapefactor.py b/src/pyhf/modifiers/shapefactor.py index f73ac8e050..d2d3e1ed70 100644 --- a/src/pyhf/modifiers/shapefactor.py +++ b/src/pyhf/modifiers/shapefactor.py @@ -62,8 +62,8 @@ def __init__(self, shapefactor_mods, pdfconfig, mega_mods, batch_size=None): and at that point can be used to compute the effect of shapefactor. """ - self.name = shapefactor.name - self.op_code = shapefactor.op_code + self.name = shapefactor.name + self.op_code = shapefactor.op_code self.batch_size = batch_size keys = [f'{mtype}/{m}' for m, mtype in shapefactor_mods] shapefactor_mods = [m for m, _ in shapefactor_mods] diff --git a/src/pyhf/modifiers/shapesys.py b/src/pyhf/modifiers/shapesys.py index a199fab6dc..5660eb1501 100644 --- a/src/pyhf/modifiers/shapesys.py +++ b/src/pyhf/modifiers/shapesys.py @@ -37,8 +37,8 @@ def required_parset(cls, sample_data, modifier_data): class shapesys_combined: def __init__(self, shapesys_mods, pdfconfig, mega_mods, batch_size=None): - self.name = shapesys.name - self.op_code = shapesys.op_code + self.name = shapesys.name + self.op_code = shapesys.op_code self.batch_size = batch_size keys = [f'{mtype}/{m}' for m, mtype in shapesys_mods] diff --git a/src/pyhf/modifiers/staterror.py b/src/pyhf/modifiers/staterror.py index 56bc308b88..4d306242a5 100644 --- a/src/pyhf/modifiers/staterror.py +++ b/src/pyhf/modifiers/staterror.py @@ -26,8 +26,8 @@ def required_parset(cls, sample_data, modifier_data): class staterror_combined: def __init__(self, staterr_mods, pdfconfig, mega_mods, batch_size=None): - self.name = staterror.name - self.op_code = staterror.op_code + self.name = staterror.name + self.op_code = staterror.op_code self.batch_size = batch_size keys = [f'{mtype}/{m}' for m, mtype in staterr_mods] diff --git a/src/pyhf/pdf.py b/src/pyhf/pdf.py index 4adf87e3f0..aa87fc24f7 100644 --- a/src/pyhf/pdf.py +++ b/src/pyhf/pdf.py @@ -83,13 +83,14 @@ def _finalize_parameters(user_parameters, _paramsets_requirements, channel_nbins return _sets + class nominal_helper: def __init__(self, config): self.mega_samples = {} self.config = config - def append(self,c,s,defined_samp): - self.mega_samples.setdefault(s,{'name': f'mega_{s}', 'nom': []}) + def append(self, c, s, defined_samp): + self.mega_samples.setdefault(s, {'name': f'mega_{s}', 'nom': []}) nom = ( defined_samp['data'] if defined_samp @@ -119,26 +120,27 @@ def __init__(self, config): self.config = config self.required_parsets = {} - def collect(self, thismod,nom): + def collect(self, thismod, nom): maskval = True if thismod else False mask = [maskval] * len(nom) return {'mask': mask} - def append(self,key,c,s,thismod,defined_samp): - self.mega_mods.setdefault(key, {}).setdefault(s,{}).setdefault( - 'data',{'mask': []} + def append(self, key, c, s, thismod, defined_samp): + self.mega_mods.setdefault(key, {}).setdefault(s, {}).setdefault( + 'data', {'mask': []} ) nom = ( defined_samp['data'] if defined_samp else [0.0] * self.config.channel_nbins[c] ) - moddata = self.collect(thismod,nom) + moddata = self.collect(thismod, nom) self.mega_mods[key][s]['data']['mask'] += moddata['mask'] if thismod: - self.required_parsets.setdefault(thismod['name'],[]).append(modifiers.lumi.required_parset( - defined_samp['data'], thismod['data'] - )) + self.required_parsets.setdefault(thismod['name'], []).append( + modifiers.lumi.required_parset(defined_samp['data'], thismod['data']) + ) + class normfactor_helper: def __init__(self, config): @@ -146,26 +148,29 @@ def __init__(self, config): self.config = config self.required_parsets = {} - def collect(self, thismod,nom): + def collect(self, thismod, nom): maskval = True if thismod else False mask = [maskval] * len(nom) return {'mask': mask} - def append(self,key,c,s,thismod,defined_samp): - self.mega_mods.setdefault(key, {}).setdefault(s,{}).setdefault( - 'data',{'mask': []} + def append(self, key, c, s, thismod, defined_samp): + self.mega_mods.setdefault(key, {}).setdefault(s, {}).setdefault( + 'data', {'mask': []} ) nom = ( defined_samp['data'] if defined_samp else [0.0] * self.config.channel_nbins[c] ) - moddata = self.collect(thismod,nom) + moddata = self.collect(thismod, nom) self.mega_mods[key][s]['data']['mask'] += moddata['mask'] if thismod: - self.required_parsets.setdefault(thismod['name'],[]).append(modifiers.normfactor.required_parset( - defined_samp['data'], thismod['data'] - )) + self.required_parsets.setdefault(thismod['name'], []).append( + modifiers.normfactor.required_parset( + defined_samp['data'], thismod['data'] + ) + ) + class shapefactor_helper: def __init__(self, config): @@ -173,26 +178,29 @@ def __init__(self, config): self.config = config self.required_parsets = {} - def collect(self, thismod,nom): + def collect(self, thismod, nom): maskval = True if thismod else False mask = [maskval] * len(nom) return {'mask': mask} - def append(self,key,c,s,thismod,defined_samp): - self.mega_mods.setdefault(key, {}).setdefault(s,{}).setdefault( - 'data',{'mask': []} + def append(self, key, c, s, thismod, defined_samp): + self.mega_mods.setdefault(key, {}).setdefault(s, {}).setdefault( + 'data', {'mask': []} ) nom = ( defined_samp['data'] if defined_samp else [0.0] * self.config.channel_nbins[c] ) - moddata = self.collect(thismod,nom) + moddata = self.collect(thismod, nom) self.mega_mods[key][s]['data']['mask'] += moddata['mask'] if thismod: - self.required_parsets.setdefault(thismod['name'],[]).append(modifiers.shapefactor.required_parset( - defined_samp['data'], thismod['data'] - )) + self.required_parsets.setdefault(thismod['name'], []).append( + modifiers.shapefactor.required_parset( + defined_samp['data'], thismod['data'] + ) + ) + class shapesys_helper: def __init__(self, config): @@ -200,30 +208,33 @@ def __init__(self, config): self.config = config self.required_parsets = {} - def collect(self,thismod,nom): + def collect(self, thismod, nom): uncrt = thismod['data'] if thismod else [0.0] * len(nom) maskval = [(x > 0 and y > 0) for x, y in zip(uncrt, nom)] mask = [maskval] * len(nom) return {'mask': mask, 'nom_data': nom, 'uncrt': uncrt} - def append(self,key,c,s,thismod,defined_samp): - self.mega_mods.setdefault(key, {}).setdefault(s,{}).setdefault( - 'data',{'uncrt': [], 'nom_data': [], 'mask': []} + def append(self, key, c, s, thismod, defined_samp): + self.mega_mods.setdefault(key, {}).setdefault(s, {}).setdefault( + 'data', {'uncrt': [], 'nom_data': [], 'mask': []} ) nom = ( defined_samp['data'] if defined_samp else [0.0] * self.config.channel_nbins[c] ) - moddata = self.collect(thismod,nom) + moddata = self.collect(thismod, nom) self.mega_mods[key][s]['data']['mask'] += moddata['mask'] self.mega_mods[key][s]['data']['uncrt'] += moddata['uncrt'] self.mega_mods[key][s]['data']['nom_data'] += moddata['nom_data'] if thismod: - self.required_parsets.setdefault(thismod['name'],[]).append(modifiers.shapesys.required_parset( - defined_samp['data'], thismod['data'] - )) + self.required_parsets.setdefault(thismod['name'], []).append( + modifiers.shapesys.required_parset( + defined_samp['data'], thismod['data'] + ) + ) + class staterr_helper: def __init__(self, config): @@ -231,29 +242,32 @@ def __init__(self, config): self.config = config self.required_parsets = {} - def collect(self,thismod,nom): + def collect(self, thismod, nom): uncrt = thismod['data'] if thismod else [0.0] * len(nom) mask = [True if thismod else False] * len(nom) return {'mask': mask, 'nom_data': nom, 'uncrt': uncrt} - def append(self,key,c,s,thismod,defined_samp): - self.mega_mods.setdefault(key, {}).setdefault(s,{}).setdefault( - 'data',{'uncrt': [], 'nom_data': [], 'mask': []} + def append(self, key, c, s, thismod, defined_samp): + self.mega_mods.setdefault(key, {}).setdefault(s, {}).setdefault( + 'data', {'uncrt': [], 'nom_data': [], 'mask': []} ) nom = ( defined_samp['data'] if defined_samp else [0.0] * self.config.channel_nbins[c] ) - moddata = self.collect(thismod,nom) + moddata = self.collect(thismod, nom) self.mega_mods[key][s]['data']['mask'] += moddata['mask'] self.mega_mods[key][s]['data']['uncrt'] += moddata['uncrt'] self.mega_mods[key][s]['data']['nom_data'] += moddata['nom_data'] if thismod: - self.required_parsets.setdefault(thismod['name'],[]).append(modifiers.staterror.required_parset( - defined_samp['data'], thismod['data'] - )) + self.required_parsets.setdefault(thismod['name'], []).append( + modifiers.staterror.required_parset( + defined_samp['data'], thismod['data'] + ) + ) + class histosys_helper: def __init__(self, config): @@ -261,32 +275,35 @@ def __init__(self, config): self.config = config self.required_parsets = {} - def collect(self,thismod,nom): + def collect(self, thismod, nom): lo_data = thismod['data']['lo_data'] if thismod else nom hi_data = thismod['data']['hi_data'] if thismod else nom maskval = True if thismod else False mask = [maskval] * len(nom) return {'lo_data': lo_data, 'hi_data': hi_data, 'mask': mask, 'nom_data': nom} - def append(self,key,c,s,thismod,defined_samp): - self.mega_mods.setdefault(key, {}).setdefault(s,{}).setdefault( - 'data',{'hi_data': [], 'lo_data': [], 'nom_data': [], 'mask': []} + def append(self, key, c, s, thismod, defined_samp): + self.mega_mods.setdefault(key, {}).setdefault(s, {}).setdefault( + 'data', {'hi_data': [], 'lo_data': [], 'nom_data': [], 'mask': []} ) nom = ( defined_samp['data'] if defined_samp else [0.0] * self.config.channel_nbins[c] ) - moddata = self.collect(thismod,nom) + moddata = self.collect(thismod, nom) self.mega_mods[key][s]['data']['lo_data'] += moddata['lo_data'] self.mega_mods[key][s]['data']['hi_data'] += moddata['hi_data'] self.mega_mods[key][s]['data']['nom_data'] += moddata['nom_data'] self.mega_mods[key][s]['data']['mask'] += moddata['mask'] if thismod: - self.required_parsets.setdefault(thismod['name'],[]).append(modifiers.histosys.required_parset( - defined_samp['data'], thismod['data'] - )) + self.required_parsets.setdefault(thismod['name'], []).append( + modifiers.histosys.required_parset( + defined_samp['data'], thismod['data'] + ) + ) + class normsys_helper: def __init__(self, config): @@ -294,7 +311,7 @@ def __init__(self, config): self.config = config self.required_parsets = {} - def collect(self,thismod,nom): + def collect(self, thismod, nom): maskval = True if thismod else False lo_factor = thismod['data']['lo'] if thismod else 1.0 hi_factor = thismod['data']['hi'] if thismod else 1.0 @@ -304,9 +321,9 @@ def collect(self,thismod,nom): mask = [maskval] * len(nom) return {'lo': lo, 'hi': hi, 'mask': mask, 'nom_data': nom_data} - def append(self,key,c,s,thismod,defined_samp): - self.mega_mods.setdefault(key, {}).setdefault(s,{}).setdefault( - 'data',{'hi': [], 'lo': [], 'nom_data': [], 'mask': []} + def append(self, key, c, s, thismod, defined_samp): + self.mega_mods.setdefault(key, {}).setdefault(s, {}).setdefault( + 'data', {'hi': [], 'lo': [], 'nom_data': [], 'mask': []} ) nom = ( @@ -314,20 +331,21 @@ def append(self,key,c,s,thismod,defined_samp): if defined_samp else [0.0] * self.config.channel_nbins[c] ) - moddata = self.collect(thismod,nom) + moddata = self.collect(thismod, nom) self.mega_mods[key][s]['data']['nom_data'] += moddata['nom_data'] self.mega_mods[key][s]['data']['lo'] += moddata['lo'] self.mega_mods[key][s]['data']['hi'] += moddata['hi'] self.mega_mods[key][s]['data']['mask'] += moddata['mask'] if thismod: - self.required_parsets.setdefault(thismod['name'],[]).append(modifiers.normsys.required_parset( - defined_samp['data'], thismod['data'] - )) - + self.required_parsets.setdefault(thismod['name'], []).append( + modifiers.normsys.required_parset(defined_samp['data'], thismod['data']) + ) -def _nominal_and_modifiers_from_spec(poi_name,custom_modifiers,config, spec,batch_size): +def _nominal_and_modifiers_from_spec( + poi_name, custom_modifiers, config, spec, batch_size +): # the mega-channel will consist of mega-samples that subscribe to # mega-modifiers. i.e. while in normal histfactory, each sample might # be affected by some modifiers and some not, here we change it so that @@ -353,32 +371,32 @@ def _nominal_and_modifiers_from_spec(poi_name,custom_modifiers,config, spec,batc 'shapefactor': shapefactor_helper(config), 'lumi': lumi_helper(config), 'shapesys': shapesys_helper(config), - 'staterror': staterr_helper(config) + 'staterror': staterr_helper(config), } nominal = nominal_helper(config) - for c in config.channels: for s in config.samples: helper_data = helper.get(c, {}).get(s) - defined_samp,defined_mods = (None,None) if not helper_data else helper_data - nominal.append(c,s,defined_samp) + defined_samp, defined_mods = ( + (None, None) if not helper_data else helper_data + ) + nominal.append(c, s, defined_samp) for m, mtype in config.modifiers: key = f'{mtype}/{m}' # this is None if modifier doesn't affect channel/sample. thismod = defined_mods.get(key) if defined_mods else None - modifiers_helpers[mtype].append(key,c,s,thismod,defined_samp) + modifiers_helpers[mtype].append(key, c, s, thismod, defined_samp) nominal_rates = nominal.apply() - mega_mods = {k:v.mega_mods for k,v in modifiers_helpers.items()} - + mega_mods = {k: v.mega_mods for k, v in modifiers_helpers.items()} _parset_reqs = {} for v in list(modifiers_helpers.values()) + custom_modifiers: for pname, req_list in v.required_parsets.items(): - _parset_reqs.setdefault(pname,[]) + _parset_reqs.setdefault(pname, []) _parset_reqs[pname] += req_list user_parameters = spec.get('parameters', []) @@ -389,10 +407,8 @@ def _nominal_and_modifiers_from_spec(poi_name,custom_modifiers,config, spec,batc config.channel_nbins, ) - - poi_name = poi_name - config.set_parameters(poi_name,_required_paramsets) + config.set_parameters(poi_name, _required_paramsets) standard_modifiers = { k: c( @@ -434,8 +450,7 @@ def __init__(self, spec, **config_kwargs): self.auxdata_order = [] self.nmaindata = sum(self.channel_nbins.values()) - - def set_parameters(self,poi_name,_required_paramsets): + def set_parameters(self, poi_name, _required_paramsets): self._create_and_register_paramsets(_required_paramsets) if poi_name is not None: self.set_poi(poi_name) @@ -598,9 +613,9 @@ def logpdf(self, auxdata, pars): class _MainModel: """Factory class to create pdfs for the main measurement.""" - def __init__(self, config, modifiers, nominal_rates, batch_size = None): + def __init__(self, config, modifiers, nominal_rates, batch_size=None): self.config = config - + self._factor_mods = [] self._delta_mods = [] self.batch_size = batch_size @@ -611,7 +626,7 @@ def __init__(self, config, modifiers, nominal_rates, batch_size = None): self.modifiers_appliers = modifiers - for k,v in self.modifiers_appliers.items(): + for k, v in self.modifiers_appliers.items(): if v.op_code == 'addition': self._delta_mods.append(v.name) elif v.op_code == 'multiplication': @@ -699,9 +714,7 @@ def expected_data(self, pars, return_by_sample=False): pars = tensorlib.astensor(pars) deltas, factors = self._modifications(pars) - allsum = tensorlib.concatenate( - deltas + [self.nominal_rates] - ) + allsum = tensorlib.concatenate(deltas + [self.nominal_rates]) nom_plus_delta = tensorlib.sum(allsum, axis=0) nom_plus_delta = tensorlib.reshape( @@ -726,7 +739,7 @@ def expected_data(self, pars, return_by_sample=False): class Model: """The main pyhf model class.""" - def __init__(self, spec, custom_modifiers = None, batch_size=None, **config_kwargs): + def __init__(self, spec, custom_modifiers=None, batch_size=None, **config_kwargs): """ Construct a HistFactory Model. @@ -748,20 +761,20 @@ def __init__(self, spec, custom_modifiers = None, batch_size=None, **config_kwar log.info(f"Validating spec against schema: {self.schema:s}") utils.validate(self.spec, self.schema, version=self.version) # build up our representation of the specification - self.config = _ModelConfig(self.spec,**config_kwargs) - - - + self.config = _ModelConfig(self.spec, **config_kwargs) standard_modifiers, _nominal_rates = _nominal_and_modifiers_from_spec( - config_kwargs.pop('poi_name'), custom_modifiers, self.config, self.spec, self.batch_size + config_kwargs.pop('poi_name'), + custom_modifiers, + self.config, + self.spec, + self.batch_size, ) for custom in custom_modifiers: custom.config = self.config standard_modifiers[custom.name] = custom - self.main_model = _MainModel( self.config, modifiers=standard_modifiers,