diff --git a/.zenodo.json b/.zenodo.json index 99af30bbe6..71cabfaeb2 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -28,7 +28,6 @@ "fitting", "scipy", "numpy", - "pytorch", "jax", "auto-differentiation" ], diff --git a/CITATION.cff b/CITATION.cff index e09c96f692..a1dc9ce564 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -26,7 +26,6 @@ keywords: - fitting - scipy - numpy - - pytorch - jax - auto-differentiation license: "Apache-2.0" @@ -37,8 +36,8 @@ abstract: | of that statistical model for multi-bin histogram-based analysis and its interval estimation is based on the asymptotic formulas of "Asymptotic formulae for likelihood-based tests of new physics". pyhf supports modern - computational graph libraries such as PyTorch and JAX in order - to make use of features such as autodifferentiation and GPU acceleration. + computational graph libraries such as JAX in order to make use of features + such as automatic differentiation and GPU acceleration. references: - type: article authors: diff --git a/README.rst b/README.rst index 7074aca14e..023706d716 100644 --- a/README.rst +++ b/README.rst @@ -28,9 +28,8 @@ multi-bin histogram-based analysis and its interval estimation is based on the asymptotic formulas of “Asymptotic formulae for likelihood-based tests of new physics” [`arXiv:1007.1727 `__]. The aim is also -to support modern computational graph libraries such as PyTorch and -JAX in order to make use of features such as autodifferentiation -and GPU acceleration. +to support modern computational graph libraries such as JAX in order to +make use of features such as automatic differentiation and GPU acceleration. .. Comment: JupyterLite segment goes here in docs @@ -144,7 +143,6 @@ Implemented variations: Computational Backends: - ☑ NumPy - - ☑ PyTorch - ☑ JAX Optimizers: diff --git a/docker/gpu/install_backend.sh b/docker/gpu/install_backend.sh index 238f78acbe..2995e4b379 100644 --- a/docker/gpu/install_backend.sh +++ b/docker/gpu/install_backend.sh @@ -18,10 +18,7 @@ function get_JAXLIB_GPU_WHEEL { function install_backend() { # 1: the backend option name in setup.py local backend="${1}" - if [[ "${backend}" == "torch" ]]; then - # shellcheck disable=SC2102 - python3 -m pip install --no-cache-dir .[xmlio,torch] - elif [[ "${backend}" == "jax" ]]; then + if [[ "${backend}" == "jax" ]]; then python3 -m pip install --no-cache-dir .[xmlio] python3 -m pip install --no-cache-dir "$(get_JAXLIB_GPU_WHEEL)" python3 -m pip install --no-cache-dir jax diff --git a/docs/api.rst b/docs/api.rst index a7db15f1e7..6d63313a7b 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -65,7 +65,6 @@ The computational backends that :code:`pyhf` provides interfacing for the vector :nosignatures: numpy_backend.numpy_backend - pytorch_backend.pytorch_backend jax_backend.jax_backend Optimizers diff --git a/docs/conf.py b/docs/conf.py index ca52d89898..4844d323fb 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -146,7 +146,6 @@ def setup(app): # today_fmt = '%B %d, %Y' autodoc_mock_imports = [ - 'torch', 'jax', 'iminuit', ] @@ -195,13 +194,11 @@ def setup(app): 'examples/notebooks/StatError.ipynb', 'examples/notebooks/histogrammar.ipynb', 'examples/notebooks/histosys.ipynb', - 'examples/notebooks/histosys-pytorch.ipynb', 'examples/notebooks/importxml.ipynb', 'examples/notebooks/multichannel-coupled-normsys.ipynb', 'examples/notebooks/multichannel-normsys.ipynb', 'examples/notebooks/normsys.ipynb', 'examples/notebooks/pullplot.ipynb', - 'examples/notebooks/pytorch_tests_onoff.ipynb', ] # The reST default role (used for this markup: `text`) to use for all diff --git a/docs/examples/notebooks/histosys-pytorch.ipynb b/docs/examples/notebooks/histosys-pytorch.ipynb deleted file mode 100644 index 10a83d641a..0000000000 --- a/docs/examples/notebooks/histosys-pytorch.ipynb +++ /dev/null @@ -1,134 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# HistoSys with PyTorch" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Populating the interactive namespace from numpy and matplotlib\n" - ] - } - ], - "source": [ - "%pylab inline" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import pyhf\n", - "from pyhf import Model\n", - "from pyhf.simplemodels import uncorrelated_background" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[120.0, 180.0, 100.0, 225.0]\n", - "[[0, 10], [0, 10], [0, 10]]\n", - "['mu', 'uncorr_bkguncrt']\n" - ] - } - ], - "source": [ - "source = {\n", - " \"binning\": [2, -0.5, 1.5],\n", - " \"bindata\": {\n", - " \"data\": [120.0, 180.0],\n", - " \"bkg\": [100.0, 150.0],\n", - " \"bkgerr\": [10.0, 10.0],\n", - " \"sig\": [30.0, 95.0],\n", - " },\n", - "}\n", - "\n", - "pdf = uncorrelated_background(\n", - " source['bindata']['sig'], source['bindata']['bkg'], source['bindata']['bkgerr']\n", - ")\n", - "data = source['bindata']['data'] + pdf.config.auxdata\n", - "\n", - "init_pars = pdf.config.suggested_init()\n", - "par_bounds = pdf.config.suggested_bounds()\n", - "\n", - "print(data)\n", - "print(par_bounds)\n", - "print(pdf.config.par_order)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - "# NumPy\n", - " [-23.57960517]\n", - "\n", - "# PyTorch\n", - " tensor([-23.5796])\n" - ] - } - ], - "source": [ - "backends = [\n", - " pyhf.tensor.numpy_backend(),\n", - " pyhf.tensor.pytorch_backend(),\n", - "]\n", - "names = [\n", - " 'NumPy',\n", - " 'PyTorch',\n", - "]\n", - "\n", - "for backend, name in zip(backends, names):\n", - " print(f'\\n# {name}')\n", - " pyhf.set_backend(backend)\n", - " v = pdf.logpdf(init_pars, data)\n", - " print(type(v), v)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "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.6.6" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/docs/examples/notebooks/pytorch_tests_onoff.ipynb b/docs/examples/notebooks/pytorch_tests_onoff.ipynb deleted file mode 100644 index b1830d14ad..0000000000 --- a/docs/examples/notebooks/pytorch_tests_onoff.ipynb +++ /dev/null @@ -1,347 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# PyTorch Tests On/Off" - ] - }, - { - "cell_type": "code", - "execution_count": 33, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "import numpy as np\n", - "import torch\n", - "from torch.autograd import Variable\n", - "\n", - "\n", - "def _log_poisson_impl(n, lam):\n", - " normal = torch.distributions.Normal(lam, torch.sqrt(lam))\n", - " return normal.log_prob(n)\n", - "\n", - "\n", - "class shapesys_constraint:\n", - " def __init__(self, nom_data, mod_data):\n", - " self.auxdata = []\n", - " self.bkg_over_db_squared = []\n", - " for b, deltab in zip(nom_data, mod_data):\n", - " bkg_over_bsq = b * b / deltab / deltab # tau*b\n", - " self.bkg_over_db_squared.append(bkg_over_bsq)\n", - " self.auxdata.append(bkg_over_bsq)\n", - "\n", - " def alphas(self, pars):\n", - " return np.prod([pars, self.bkg_over_db_squared], axis=0)\n", - "\n", - " def logpdf(self, a, alpha):\n", - " return _log_poisson_impl(a, alpha)\n", - "\n", - " def expected_data(self, pars):\n", - " return self.alphas(pars)\n", - "\n", - "\n", - "class Model:\n", - " def __init__(self, spec):\n", - " self.samples = {}\n", - " self.samples['sig'] = spec['sig']\n", - " self.samples['bkg'] = spec['bkg']\n", - " self.constraint = shapesys_constraint(spec['bkg'], spec['bkgerr'])\n", - "\n", - " def expected_actualdata(self, pars):\n", - " return [\n", - " pars[0] * s + pars[1] * b\n", - " for s, b in zip(self.samples['sig'], self.samples['bkg'])\n", - " ]\n", - "\n", - " def expected_auxdata(self, pars):\n", - " modparslice = slice(1, 2)\n", - " return self.constraint.expected_data(pars[modparslice])\n", - "\n", - " def expected_data(self, pars, include_auxdata=True):\n", - " expected_actual = self.expected_actualdata(pars)\n", - "\n", - " if not include_auxdata:\n", - " return expected_actual\n", - " expected_constraints = self.expected_auxdata(pars)\n", - " return np.concatenate([expected_actual, expected_constraints])\n", - "\n", - " def logpdf(self, data, pars):\n", - " actual_data, auxdata = [data[0]], [data[1]]\n", - " lambdas_data = self.expected_actualdata(pars)\n", - "\n", - " sum = 0\n", - " for d, lam in zip(actual_data, lambdas_data):\n", - " sum = sum + _log_poisson_impl(d, lam)\n", - "\n", - " modauxslice = slice(0, 1)\n", - " thisauxdata = auxdata[slice(0, 1)]\n", - "\n", - " modparslice = slice(1, 2)\n", - " modalphas = pdf.constraint.alphas(pars[modparslice])\n", - "\n", - " for a, alpha in zip(thisauxdata, modalphas):\n", - " sum = sum + self.constraint.logpdf(a, alpha)\n", - " return sum\n", - "\n", - "\n", - "def loglambdav(pars, data, pdf):\n", - " return -2 * pdf.logpdf(data, pars)\n", - "\n", - "\n", - "def optimize(\n", - " objective,\n", - " init_pars,\n", - " par_bounds,\n", - " trf_to_unbounded=lambda x: x,\n", - " trf_from_unbounded=lambda x: x,\n", - " args=None,\n", - " constrain_index_value=None,\n", - "):\n", - " data = args[0]\n", - "\n", - " poi_index = None\n", - " constrained_var = None\n", - " if constrain_index_value is not None:\n", - " poi_index, poi_value = constrain_index_value\n", - "\n", - " upto = init_pars[:poi_index]\n", - " fromon = init_pars[poi_index + 1 :]\n", - " nuis_init = np.concatenate([upto, fromon])\n", - "\n", - " init_pars = nuis_init\n", - " constrained_var = Variable(torch.Tensor([poi_value]), requires_grad=True)\n", - "\n", - " init_pars = trf_to_unbounded(init_pars)\n", - "\n", - " # print('init from: %s' % init_pars)\n", - "\n", - " parameters = Variable(torch.Tensor(init_pars), requires_grad=True)\n", - " data = Variable(torch.Tensor(data))\n", - "\n", - " optimizer = torch.optim.Adam([parameters])\n", - "\n", - " def assemble_from_nuis(nuisance_pars):\n", - " tocat = []\n", - " if poi_index > 0:\n", - " tocat.append(nuisance_pars[:poi_index])\n", - " tocat.append(constrained_var)\n", - " if poi_index < len(parameters):\n", - " tocat.append(nuisance_pars[poi_index:])\n", - " return torch.cat(tocat)\n", - "\n", - " for i in range(1000):\n", - " objpars = None\n", - " if constrain_index_value is None:\n", - " objpars = parameters\n", - " else:\n", - " objpars = assemble_from_nuis(parameters)\n", - " loss = objective(objpars, data, args[1])\n", - " optimizer.zero_grad()\n", - " loss.backward()\n", - " optimizer.step()\n", - " # if i % 1000 == 0:\n", - " # print(trf_from_unbounded(parameters).data.numpy())\n", - "\n", - " result = parameters\n", - " if constrain_index_value is not None:\n", - " result = assemble_from_nuis(result)\n", - " result = trf_from_unbounded(result).data.numpy()\n", - " # print('result %s' % result)\n", - " return result\n", - "\n", - "\n", - "def unconstrained_bestfit(data, pdf, init_pars, par_bounds):\n", - " # print 'fit', data, 'expected for init pars', np.array(pdf.expected_actualdata(init_pars))\n", - " result = optimize(\n", - " loglambdav,\n", - " init_pars,\n", - " par_bounds,\n", - " args=(\n", - " data,\n", - " pdf,\n", - " ),\n", - " )\n", - " return result\n", - "\n", - "\n", - "def constrained_bestfit(constrained_mu, data, pdf, init_pars, par_bounds):\n", - " # print 'fit', data, 'for pars', np.array(pdf.expected_actualdata(init_pars))\n", - " result = optimize(\n", - " loglambdav,\n", - " init_pars,\n", - " par_bounds,\n", - " args=(\n", - " data,\n", - " pdf,\n", - " ),\n", - " constrain_index_value=(0, constrained_mu),\n", - " )\n", - " return result" - ] - }, - { - "cell_type": "code", - "execution_count": 34, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "spec = {\n", - " 'obs': [55.0],\n", - " 'sig': [10.0],\n", - " 'bkg': [50.0],\n", - " 'bkgerr': [7.0],\n", - "}\n", - "\n", - "pdf = Model(spec)\n", - "\n", - "data = np.concatenate([spec['obs'], pdf.constraint.auxdata])\n", - "init_pars = [1.0, 1.0] # (mu, gamma)\n", - "par_bounds = [[0, 10], [0, 10]]\n", - "\n", - "# constrained_bestfit(1.0, data, pdf, init_pars, par_bounds)\n", - "# unconstrained_bestfit(data, pdf, init_pars, par_bounds)" - ] - }, - { - "cell_type": "code", - "execution_count": 38, - "metadata": {}, - "outputs": [], - "source": [ - "### The Test Statistic\n", - "def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds):\n", - " bestfit_nuisance_asimov = constrained_bestfit(\n", - " asimov_mu, data, pdf, init_pars, par_bounds\n", - " )\n", - " return pdf.expected_data(bestfit_nuisance_asimov)\n", - "\n", - "\n", - "def qmu(mu, data, pdf, init_pars, par_bounds):\n", - " mubhathat = constrained_bestfit(mu, data, pdf, init_pars, par_bounds)\n", - " muhatbhat = unconstrained_bestfit(data, pdf, init_pars, par_bounds)\n", - "\n", - " v_mbhh = Variable(torch.Tensor(mubhathat))\n", - " v_mhbh = Variable(torch.Tensor(muhatbhat))\n", - " v_data = Variable(torch.Tensor(data))\n", - " qmu = loglambdav(v_mbhh, v_data, pdf) - loglambdav(v_mhbh, v_data, pdf)\n", - " qmu = qmu.data.numpy()[0]\n", - " if muhatbhat[0] > mu:\n", - " return 0.0\n", - " if -1e-6 < qmu < 0:\n", - " log.warning('WARNING: qmu negative: %s', qmu)\n", - " return 0.0\n", - " return qmu\n", - "\n", - "\n", - "def pvals_from_teststat(sqrtqmu_v, sqrtqmuA_v):\n", - " from scipy.stats import norm\n", - "\n", - " CLsb = 1 - norm.cdf(sqrtqmu_v)\n", - " CLb = norm.cdf(sqrtqmuA_v - sqrtqmu_v)\n", - " CLs = CLb / CLsb\n", - " return CLsb, CLb, CLs\n", - "\n", - "\n", - "def runOnePoint(muTest, data, pdf, init_pars, par_bounds):\n", - " asimov_mu = 0.0\n", - " asimov_data = generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds)\n", - "\n", - " qmu_v = qmu(muTest, data, pdf, init_pars, par_bounds)\n", - " qmuA_v = qmu(muTest, asimov_data, pdf, init_pars, par_bounds)\n", - "\n", - " sqrtqmu_v = np.sqrt(qmu_v)\n", - " sqrtqmuA_v = np.sqrt(qmuA_v)\n", - "\n", - " sigma = muTest / sqrtqmuA_v if sqrtqmuA_v > 0 else None\n", - "\n", - " CLsb, CLb, CLs = pvals_from_teststat(sqrtqmu_v, sqrtqmuA_v)\n", - "\n", - " CLs_exp = []\n", - " for nsigma in [-2, -1, 0, 1, 2]:\n", - " sqrtqmu_v_sigma = sqrtqmuA_v - nsigma\n", - " CLs_exp.append(pvals_from_teststat(sqrtqmu_v_sigma, sqrtqmuA_v)[-1])\n", - " return qmu_v, qmuA_v, sigma, CLsb, CLb, CLs, CLs_exp" - ] - }, - { - "cell_type": "code", - "execution_count": 44, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Populating the interactive namespace from numpy and matplotlib\n" - ] - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XlU1Nf5+PH3nRkY9kU2WUXAHRSUAYwKLkRB4opEE7ekidakSZMmX5ukbdI2vyZpkia1aWsSY02r0aiIiHFDoxijMe5GDcFdFNxFcUMR+fz+wCGoLDPDIDPDfZ3Tc4T5zOc+k3P6zOXe+3keoSgKkiRJkm1RNXcAkiRJkvnJ5C5JkmSDZHKXJEmyQTK5S5Ik2SCZ3CVJkmyQTO6SJEk2qMHkLoSYJYQ4K4TYV8frQgjxkRDikBBijxCiu/nDlCRJkoxhyMz9v0BKPa+nAu3u/G8y8HHjw5IkSZIao8HkrijKBqCknkuGAbOVKt8DHkIIf3MFKEmSJBlPY4Z7BAInavxcdOd3p+69UAgxmarZPc7Ozj06duxo9GAnT56s/rdGo8He3h57e3u0Wi1qtdro+0mSJFmTHTt2nFcUxaeh68yR3A2mKMoMYAZAbGyssn37dqPvsWHDBvLy8gDw8fHh6tWrlJWVAeDm5kZoaCht2rShTZs2tGrVCiGE+T6AJElSMxNCFBpynTmSezEQXOPnoDu/axIJCQnk5eWhKArnzp1j0KBBhIWFcezYMQoLCzl8+DB79uwBwMXFpTrZh4aG4uXlJZO9JEktgjmS+1LgOSHEfCAeKFUU5b4lGXPRL8HcvHmTCxcukJuby+DBg4mLiyMuLg5FUbhw4UJ1sj927Bj79lUd9HF2dq6e1YeGhuLj4yOTvSRJNqnB5C6E+BLoC3gLIYqAPwJ2AIqifAKsAAYDh4DrwJNNFaxeQEAAR48excHBgYqKClasWIFaraZ79+4IIfD29sbb25vY2FgURaGkpITCwsLqZJ+fnw+Ak5NTdbJv06YNfn5+MtlLkmQTGkzuiqI81sDrCvArs0VkgA4dOnD06FGcnZ3JyclhypQpfPXVV6jVarp163bXtUIIvLy88PLyonv37iiKwqVLl+5K9j/99BMADg4OdyX71q1bo1LJ57wkSbI+D3RD1VzCwsIAUKlU6HQ6Vq5cyciRI8nJyUGlUhEVFVXne4UQeHp64unpSXR0NAClpaXVyziFhYXs378fAK1WS0hISPUyjr+/v0z2kiRZBatM7t7e3qjVam7fvk1AQADLly9nypQpVFZWkp2djVqtpnPnzgbfz93dnW7dulXP+i9fvlyd6AsLCzl48CBQtd4fFxdHUlISGo1V/qeTJKmFsMoMJYTA19eXU6dOodFoGDBgAC+//DJbt25l/vz5ZGVloVKpMOUcPVQdqYyKiqr+C+Dq1asUFhZSUFDAxo0bKSgoYNiwYQQFBZnzY0mSJJmN1a4xtG3bFgBPT0/CwsI4fvw4//3vfxk7diz+/v5kZmZWz7gby8XFhS5dupCens7YsWMpLy9n1qxZrF69mlu3bpllDEmSJHOy2uQeHFx1tP7q1asIIRg1ahRvvPEG169fZ9y4cfj5+bFgwQIOHz5s1nEjIiJ49tln6d69O5s3b+bTTz/l+PHjZh1DkiSpsaw2ueuXRMrKyggNDaVdu3ZcuXKFP//5zzg4ODB+/Hi8vb2ZP38+R48eNevYWq2WRx55hPHjx3P79m0+//xzVq1aJWfxkiRZDKtN7i4uLri6ugJVG6zl5eU8++yz/Pvf/6agoABHR0fGjx+Pp6cnX375JYWFBj2xa5SwsDCeeeYZdDodW7Zs4eOPP+bYsWNmH0eSJMlYVpvcoWppRqVSUVJSgp+fH2FhYTg5OfHSSy8BVU+kTpgwAXd3d+bNm8eJEycauKPx7O3tGTx4MBMnTgTgf//7HytWrKC8vNzsY0mSJBnKqpN7YGAglZWVFBYWotPpKCkp4fe//z0rV65kxYoVQNUMf8KECbi4uDB37ty7qkqaU2hoKFOmTCE+Pp5t27bx8ccfc+TIkSYZS5IkqSFWn9wBKioqcHV1xcnJiaCgINq3b89LL71UvQbu6urKhAkTcHR0ZM6cOZw61TSlb+zt7UlJSeHJJ59EpVIxZ84cli1bxs2bN5tkPEmSpLpYdXIPCAgAqs69Hzt2jNjYWA4ePMjbb7/N/v37+fe//119rbu7OxMnTkSr1TJnzhzOnDnTZHGFhIQwZcoUevbsyc6dO5k+fbrZT+1IkiTVx6qTu52dHa1bt0ar1XL48GFiY2NRqVQ4OzszcOBA/vznP3P+/Pnq6z08PJgwYQIajYbZs2dz7ty5Jo1t4MCB/OIXv8De3p4vvviCnJwcbty40WRjSpIk6Vl1coeq2futW7c4e/YsiqLQpUsXdu3axXvvvceVK1d444037rq+VatWTJw4EZVKxezZs7lw4UKTxhcUFMQvf/lLevXqxQ8//MD06dM5cOBAk44pSZJk9ck9KCiI27dvA3D48GESEhIoLy/nxo0bPPPMM3z66afs3bv3rvd4eXkxYcIEKisr+d///kdJSX0tYhtPo9GQnJzM008/jaOjI19++SVLliyp7iAlSZJkblaf3PWbqlqtliNHjhAQEEBwcDBbtmzhjTfewN3dnRdffJGqysQ/8/HxYcKECVRUVDB79mwuXbrU5LEGBAQwadIk+vTpw549e5g+fToFBQVNPq4kSS2P1Sd3b29v7O3tcXFx4fDhw1RWVhIfH8/FixcpKSnhzTffZN26deTk5Nz3Xj8/P8aPH8/Nmzf53//+R2lpaZPHq9Fo6N+/P5MmTcLZ2ZkFCxaQlZXF9evXm3xsSZJaDqtP7iqVioCAAG7fvk1ZWRmnTp2iU6dOuLm5sWXLFqZMmULnzp35v//7v1qPJPr7+zNu3DjKysqYPXs2V65ceSBx+/v7M2nSJPr27Ut+fj7Tp0+v7hAlSZLUWFaf3KFqaUY/6z58+HB1E4+jR49y4cIF/v73v3P48GH+8Y9/1Pn+sWPHcvXqVWbPns3Vq1cfSNxqtZqkpCQmT56Mm5sbmZmZZGZmcu3atQcyviRJtssmkntQUBCKouDl5VV9nrxHjx5oNBq2bNnCwIEDeeSRR/jLX/5S5/n24OBgHn/8cUpLS5k9e/YDTbB+fn489dRT9O/fn/379zN9+nT27dt33z6BJEmSoWwiues3Vd3c3Dhx4gQ3btzA0dGRbt26sXfvXq5fv84HH3zAjRs3+P3vf1/nfdq0acNjjz3GxYsXmTNnzgM9zaJWq+nTpw+TJ0/G09OTrKwsFi5c+MD+ipAkybbYRHJ3dXXFzc0NAEVRqkv8xsfHU1FRwY4dO2jfvj3PP/88s2bNYteuXXXeq23btowZM4bz588zZ86cB/7Qka+vL7/4xS9ITk7m4MGDTJ8+nT179shZvCRJRrGJ5A5Vs/eSkhLs7e2rl2Z8fHwICwtj27Zt3L59m9dffx0vLy9eeOGFepNleHg4o0eP5syZM3zxxRcPvDaMSqWiV69eTJkyBS8vL7Kzs8nMzKSiouKBxiFJkvWyqeReWlpKcHAwhw8frk7eCQkJXLlyhZ9++gkPDw/eeustvv32WzIzM+u9X7t27cjIyODUqVPMnTu3WUr4ent78+STT5KcnMxPP/1EVlZW9QNbkiRJ9bGZ5K7vzOTp6cmlS5eqnzqNiIjAy8uL77//HoCnnnqKbt26MXXq1AbX1Dt27Eh6ejpFRUXMmzevWTot6WfxqampFBQUsGTJEiorKx94HJIkWRebSe7+/v4IIRBCAHDo0CGgqmJkXFwcxcXFFBUVoVarmTZtGsePH+eDDz5o8L6dO3dmxIgRHD9+nPnz5zdbK724uDgGDBjAvn37WLZsmVyDlySpXjaT3O3t7fH19aWkpARPT8+7Sux269YNrVbLli1bAOjbty8jR47knXfeobi4uMF7R0VFMWzYMI4cOcKaNWua7DM0pHfv3iQmJrJr1y5WrVolE7wkSXWymeQOVevuxcXFhIeHc+zYseoNSK1WS0xMDPn5+Vy+fBmA999/n4qKCl577TWD7t2tWzd0Oh3bt29v0lrwDenbty8JCQls3bqVdevWNVsckiRZNptL7jdu3MDX15dbt27d1TM1Li4ORVHYtm0bUNXc+qWXXmLOnDnVM/qG9OvXDwcHh2adNQshGDhwID169GDjxo1s2LChWeKQJMmy2VRy12+qCiFQqVTV6+5QtdHaoUMHduzYUb1u/rvf/Y7WrVvzwgsvGLRJ6ejoSP/+/Tl27Fiz1oERQpCWlkbXrl3Jy8ur3iyWJEnSs6nkrq8Qefbs2eojkTXFx8dTVlZWXd/d1dWVd955hy1btjBv3jyDxujevTutW7dm9erVzba5ClUJftiwYXTu3Jnc3Fx27NjRbLFIkmR5bCq56ytE6tfdz5w5c9fj+23atMHPz48tW7ZUL6tMmDCB2NhYXn31VYPqyahUKlJTU7l8+TIbN25sss9iCJVKxciRI2nXrh3Lli1jz549zRqPJEmWw6aSO1Stu58+fZrQ0FCAu2bvQgji4+M5e/Ysx44dA6oS5LRp0yguLubdd981aIyQkBCioqLYtGkTFy9eNPdHMIparebRRx+lbdu2LFmyRJYNliQJsNHkrl8/d3Jyum9pJioqCicnp7s2UXv16sWYMWN4//33KSwsNGic5ORkVCoVq1evNl/wJtJoNIwZM4agoCCysrI4ePBgc4ckSVIzs7nkrt9U1S/N1CxFAFWJMDY2lv3799/VO/Xdd99FCMErr7xi0Dhubm706dOHgoKC+75AmoO9vT2PP/44fn5+LFiwoLp4miRJLZPNJXdXV1dcXV2rk/v169c5derUXdfExsaiUqnYunVr9e9CQkKYOnUqCxYsMHgtvWfPnnh6erJq1SqLqPni4ODAuHHj8PLy4ssvv7zrKKgkSS2LQcldCJEihNgvhDgkhHi1ltdDhBB5QohdQog9QojB5g/VcEFBQdXJHbhvZu3q6kpkZCS7du26q+Ljb3/7WwIDAw0+GqnRaBg0aBDnz5+vPj/f3JycnBg/fjxubm7MnTuXkydPNndIkiQ1gwaTuxBCDfwbSAU6A48JITrfc9kfgIWKosQAY4Dp5g7UGIGBgVy8eBEhBK1bt6512SQ+Pp7y8vK7ars7Ozvz3nvvsXPnTv773/8aNFb79u2JiIhg/fr1FtMez8XFhfHjx+Pg4MAXX3zB2bNnmzskSZIeMENm7nHAIUVRjiiKUg7MB4bdc40CuN35tzvQrNNFfWemkydPEh4ezokTJ+6ryR4QEEBwcDBbt269a5b+2GOP0bNnT373u99VlyqojxCCQYMGcevWLdauXWveD9II7u7uTJw4EY1Gw+zZs7lw4UJzhyRJ0gNkSHIPBGou3hbd+V1NfwLGCSGKgBXA87XdSAgxWQixXQix/dy5cyaEa5iAgACEEBQVFREeHk5lZWX10cea4uPjuXjx4l2nS4QQ/OMf/+DMmTO8/fbbBo3n7e1NQkICu3btMqgQ2YPi6enJhAkTUBSF2bNnc+nSpeYOSZKkB8RcG6qPAf9VFCUIGAzMEULcd29FUWYoihKrKEqsj4+PmYa+n729PT4+PhQXFxMcHIydnd1dpQj0OnXqhJub2321ZXQ6HRMmTODvf/+7wSdhEhMTcXFxsbhqjd7e3owfP57y8nJmz55t0F8jkiRZP0OSezEQXOPnoDu/q+kpYCGAoiibAQfA2xwBmkpfIVKtVhMaGlprklapVOh0Oo4ePXpfpcd33nkHOzs7pk6datB4Wq2W5ORkioqK+OGHH8zyGcyldevWjBs3jmvXrjFnzhyL2RuQJKnpGJLctwHthBBthRD2VG2YLr3nmuPAAAAhRCeqknvTrbsYICgoiBs3blBSUkJ4eDgXL16861y7Xo8ePdBoNPfN3gMCAnjttdfIzs42uLRu165dCQoK4uuvv37gfVcbEhgYyOOPP86lS5eYM2dOg12oJEmybg0md0VRKoDngFzgJ6pOxfwohHhTCDH0zmUvA5OEED8AXwJPKM28NqHfVC0uLiYiIgK4/0gkVFV67NatG3v37uX69et3vfbSSy/Rpk0bXnzxRYOaUwshSElJ4dq1a3zzzTdm+BTm1aZNG8aMGcP58+eZO3euxX0BSZJkPgatuSuKskJRlPaKooQrivLWnd+9oSjK0jv/zlcUpZeiKN0URYlWFKXZn8n38fHBzs6OoqIiWrVqhYeHR53r5/Hx8VRUVNxXWdHR0ZG//e1v7N27l5kzZxo0bmBgIDExMWzZsoXz5883+nOYW3h4eHXj7+bqCytJUtOzuSdU9WpWiBRCEB4eztGjR2t9ktTHx4ewsDC2bdt23+vp6ekkJiby+uuvG3zaZMCAAdjZ2Vnc5qpehw4dGDFiBCdOnGDBggUG/VUiSZJ1sdnkDj9XiKyoqCA8PJzy8vI6H8lPSEjgypUr91VVFEIwbdo0Lly4wJtvvmnQuM7OzvTt25fDhw9z4MCBRn+OphAZGcnQoUM5fPgwixYtsojyCZIkmY9NJ/egoCAqKys5ffo0bdu2RQhR59JMREQEXl5etbbci4mJ4amnnuKf//wn+/fvN2hsnU6Hj48Pubm5Fjszjo6OZvDgwezfv5/s7GyDSi5IkmQdbDq519xUdXBwqLU7k54Qgri4OIqLiykqKrrv9b/85S84ODjw5z//2aCx1Wo1KSkpXLx4kc2bN5v+IZqYTqcjOTmZH3/8ka+++soil5EkSTKeTSd3Nze36gqRULWZeOrUqTrPeUdHR6PVamudvfv5+TFlyhQWLFhg8INNYWFhdOrUiW+//daiHx7q1asXSUlJ7N69m5UrV8oEL0k2wKaTO1TN3vUz8bqqROrZ29sTExNDfn5+rcn4N7/5DRqNhr/97W8Gjz9w4EAURWHNmjUmRP/gJCUl0bNnT7Zt28bXX38tE7wkWbkWkdwvXrzI9evX8ff3x9HRsd6Zd1xcHIqi1FrCNyAggCeeeILPP/+c06dPGzS+h4cHDz30EPv27TO4y1NzEELw8MMPExsby3fffceGDRuaOyRJkhrB5pN7zc5MKpWq1u5MNXl6etKhQwd27NhR6xnwqVOncuvWLf7+978bHEPv3r1xc3Nj5cqVFr1pKYRg8ODBREdHs379er777rvmDkmSJBPZfHL39/cHuGvd/dq1a/fVkqkpPj6esrIy9u7de99rERERZGRk8PHHHxt87t3Ozo6BAwdy5swZdu7cacKneHCEEAwZMoQuXbqwZs0ai2lCIkmScWw+uWu1Wnx9fe9K7kCtVSL12rRpg5+fH1u2bKl1hv/qq69y5coVpk83vCdJ586dCQ0NZd26dRZf10WlUjFixAjat2/PihUr7jv7L0mS5bP55A4/V4hUFAVXV1d8fX3rXXcXQpCQkMDZs2drrQMfHR1Namoq06ZNu68eTX33TElJ4caNG+Tl5Zn6UR4YtVpNRkYGgYGB5OTk1Fp0TZIky9VikntZWRkXL14Eqmbvx48fp7y8vM73REZG4uTkxPfff1/r66+++irnzp1j1qxZBsfh5+eHTqdj+/btBm/INieNRkNGRgZqtZqFCxfKOjSSZEVaTHIHqo9ERkRE1NmdSU+j0RAbG8uBAwdqnbX26dOHhx56iPfff9+opNe3b18cHBwstu7Mvdzd3RkxYgRnzpxh5cqVzR2OJEkGahHJ3dfXFzs7u+p195CQEDQaTYMPI8XGxqJSqdi6det9rwkheO211zh+/Djz5883OBZHR0cGDBhAYWEhP/74o3EfpJm0a9eO3r17s2vXLotrRCJJUu1aRHKvWSESqmbldXVnqsnV1ZXIyEh27dpVa+3ztLQ0IiMj+etf/2rUEceYmBhat27NmjVr6l0asiT9+vWjTZs2LF++nLNnzzZ3OJIkNaBFJHe4u0IkVK27X7hwocHjjPHx8ZSXl7Nr1677XhNC8Oqrr5Kfn89XX31lcCwqlYrU1FQuX77Mxo0bjfsgzUSlUpGeno69vT2ZmZlW86UkSS1Vi0rut2/frj7fru/OVN+RSKh6KjU4OJitW7fWOjsfPXo0bdu25Z133jFqDT0kJISoqCi+++676o1eS+fq6kp6ejoXLlxg2bJlVrFnIEktVYtK7vDzpqqXlxfu7u4GFQGLj4/n4sWLHDx48L7XNBoNU6dOZcuWLUa31ktOTkalUrF6dbM3rjJY27Zt6du3L3v37r2vc5UkSZajxSR3Nzc3XFxcqtfdhRCEhYXV2Z2ppk6dOuHm5lZrtUiAJ598Ej8/P9555x2jY0pMTKSgoMDgSpOWoE+fPkRERLBq1SpOnjzZ3OFIklSLFpPchRAEBQVVJ3eoWpq5efPmXb+rjUqlQqfTcfTo0VrLFjg4OPDiiy+yevVqo2ezCQkJtGrVilWrVllNNyQhBCNGjMDZ2ZnMzExu3LjR3CFJknSPFpPcoWpppqSkpPrxf313pobW3QF69OiBRqOp9VgkwDPPPIObmxt//etfjYpJo9EwaNAgzp8/X+e9LZGTkxOjRo3i8uXL5OTkyPV3SbIwLS65w89FxBwdHQkMDDRoScTR0ZGoqCj27t1b60zV3d2dX/3qV2RlZRndN7V9+/a0a9eO9evXc/XqVaPe25yCg4N5+OGHKSgosOhuU5LUErWo5B4QEABwVxu98PBwTp48aVCNGJ1Ox61bt9i9e3etr7/wwgtotVree+89o2MbNGgQFRUVrF271uj3Nqf4+Hg6derE119/zfHjx5s7HEmS7mhRyV2r1eLj43PXJqD+SOSRI0cafL+/vz9BQUFs37691mUIPz8/fvGLXzB79uxa+7DWx8vLi4SEBHbv3t3gHoAlEUIwdOhQPDw8WLRoUZ0tDCVJerBaVHKHn9vu6ZNzQEAADg4OBp9W0el0XLhwgaNHj9b6+tSpU6msrOTDDz80OrbExERcXFysro+pg4MDGRkZXL9+nezsbKuKXZJsVYtM7jUrRKpUKsLCwurtzlRT586dcXJyqrOJRWhoKI899hgzZszgwoULRsWm1WpJTk6muLjY6mq4+Pv7k5qayuHDh2WLPkmyAC0uuddsu6cXHh7OlStXDKqZotFoiImJYf/+/ZSWltZ6zSuvvMK1a9f417/+ZXR8Xbt2JSgoiK+//trqjhh2796drl27sn79eoOWuSRJajotLrnfWyESfl53N3RpJjY2FkVR6jzTHhkZyZAhQ/joo4+MPv0ihCA1NZVr164Z/cRrcxNCkJaWhre3N4sXL+bKlSvNHZIktVgtLrmrVCr8/f3vSu5ubm74+PgYnNw9PDxo3749O3furC5Edq/XXnuNkpISPvvsM6NjDAgIoHv37mzdupVz584Z/f7mZG9vz6OPPkp5eTlZWVkW3RBckmxZi0vuULXufurUqbueCA0PD6ewsNDgxhs6nY5r167x008/1fp6z549SUpK4oMPPjCpgmL//v2xs7MjNzfX6jYofXx8eOSRRygsLGTdunXNHY4ktUgtNrnfvn37rlZ34eHh3L59u97uTDWFh4fTqlWrOjdWoWr2XlxczBdffGF0jM7OzvTr14/Dhw+zf/9+o9/f3Lp27Ur37t3ZtGmT0Q91SZLUeC0yude2qdqmTRuDujPpCSGIjY3lxIkTdfZDHThwIDExMbz77rsm1Y2JjY3Fx8eH3NzcOpd/LFlqaiqtW7cmOzu7wbr5kiSZV4tM7vdWiASws7OjTZs2RlVnjI6ORqPR1Dl71zfzOHDgANnZ2UbHqVarSUlJ4dKlS3z33XdGv7+56RtsK4pCZmamVX5BSZK1apHJXQhBYGDgfU+ChoeHc/78+TqPON6roXozAOnp6bRr146//vWvJq2dh4WF0alTJ7799lurnP22atWKYcOGcfLkSdasWdPc4UhSi2FQchdCpAgh9gshDgkhXq3jmkeFEPlCiB+FEPPMG6b5BQYGcuHCheoKkVCV3MHwI5Hwc72Zuh46UqvV/Pa3v2XHjh18/fXXJsU6aNAghBCsXLnSpPc3t06dOpGQkMDWrVutpim4JFm7BpO7EEIN/BtIBToDjwkhOt9zTTvgNaCXoihdgBebIFazurdCJFSd8nB1dTUquevrzWzbtq3Omfn48eMJCAgwupmHnru7O0lJSRw4cICCggKT7tHckpOTCQoKYunSpUY/uStJkvEMmbnHAYcURTmiKEo5MB8Yds81k4B/K4pyEUBRlIYf9WxmtSV3IQTh4eEcOXLEqPPZsbGx9dab0Wq1vPTSS+Tl5dXZzakhCQkJ+Pr6smrVKqtsTq1Wqxk1ahRqtZqFCxcafORUkiTTGJLcA4ETNX4uuvO7mtoD7YUQm4QQ3wshUmq7kRBishBiuxBie3M/nKOvEHnvuntERAQ3btwwqjJjly5d6q03AzB58mQ8PT2Nbuahp1arSUtLo7S01Gprt7i7uzNy5EjOnj3LihUrmjscSbJp5tpQ1QDtgL7AY8BnQgiPey9SFGWGoiixiqLE+vj4mGlo0+k3VWsup4SFhSGEMGppxpB6M66urjz//PMsWbKE/Px8k+INCQkhOjqazZs3G1QHxxJFRETQp08fdu/eXWddfEmSGs+Q5F4MBNf4OejO72oqApYqinJLUZSjwAGqkr1FCwwM5Pr163edQnF0dCQgIMDohtUN1ZsBeP7553FycuLdd981OeaHH34YrVbLihUrrO7JVb2+ffsSGhrK8uXLa+1JK0lS4xmS3LcB7YQQbYUQ9sAYYOk91yyhataOEMKbqmUaiy8LWNu6O1SdmikuLr7rJE1DatabqeuBJW9vbyZNmsS8efMoLCw0KWYnJyeSk5MpLCy0urLAeiqVivT0dBwcHMjMzOTmzZvNHZIk2ZwGk7uiKBXAc0Au8BOwUFGUH4UQbwohht65LBe4IITIB/KAqYqiWPyRCD8/PzQazX1dk8LDw1EUxeiytfp6M/Utu7z88ssIIfjggw9MihkgJiaGoKAg1qxZY9QXkCVxcXEhPT2dkpISli1bZrV/hUiSpTJozV1RlBWKorRXFCVcUZS37vzuDUVRlt75t6IoykuKonRWFCVKUZT5TRkRv2G6AAAgAElEQVS0uahUKgICAu6buQcFBaHVao1emtHXm9m+fXud1wQHBzNu3DhmzpxpcsVHIQSPPPIIZWVlJp+dtwShoaH069ePffv21fvfTJIk47XIJ1RrCggIuK9CpLHdmfT09WaOHz9e71ryb3/7W27cuME//vEPk+P28/MjPj6enTt3Gt2v1ZL07t2bdu3akZube1dvW0mSGqfFJ/egoCBu3759XzIODw/n8uXLnD9/3qj76evNbN26tc5rOnbsyIgRI/jXv/7F5cuXTYobqjYmXV1dWbZsmdXWTRdCMHz4cFxcXMjMzLTaZSZJsjQtPrnXt6kKcOjQIaPu5+joSGRkZL31ZqCqHHBpaSmffvqpkRH/TKvVkpKSwpkzZ+r9MrF0Tk5OjBo1isuXL5OTkyPX3yXJDFp8cnd3d8fZ2fm+5O7h4YGXl5fR6+4AcXFx9dabgaqjk8nJyXz44YeN6pXaqVMnIiIiyMvLa9RfAc0tKCiIgQMHsn//fqusgClJlqbFJ3d9hcja1q2N7c6k5+/vT2BgYL31ZgBeffVVTp8+zf/+9z+j49bT91ytrKxk9erVJt/HEsTFxdGpUyfWrl3L8ePHmzscSbJqLT65Q+0VIqHqacqKigqTEo1Op6u33gxUtdLT6XS89957jap13qpVK3r37s2PP/5o0l8alkIIwdChQ/H09GTRokVcu3atuUOSJKslkzs/d2a697RGmzZtUKvVRq+7g2H1ZoQQvPbaaxw5coRFixYZPUZNvXr1wsvLixUrVlh1UwwHBwcyMjK4fv06ixcvttqNYklqbjK5U3UcEu7fVLW3tyckJMSk2bAh9WYAhg0bRseOHU1u5lFzvMGDB1NSUsLGjRtNvo8laN26NYMHD+bIkSN8++23zR2OJFklmdypmi16e3vXWgkyPDycc+fOmbRZaUi9GZVKxSuvvMIPP/zQ6GYcYWFhREZGsnHjRquvmR4TE0PXrl1Zv3690U8KS5Ikk3s1/abqvbPniIgIwLjuTHqG1JsBePzxxwkODja5mUdNAwcORKPRWHVhMahaskpLS8PHx4fFixdz5cqV5g5JkqyKTO536CtE3ruE4uvri4uLi8kblfp6Mz/99FOd19jb2/N///d/bNy4sdFLKq6urvTr148jR45YfUs7e3t7MjIyKC8vJysrS66/S5IRZHK/Q7+peu+RSFO7M+mFh4fj6elZ78YqwNNPP423t7fJzTxq0ul0+Pv7k5uba/UVF318fHjkkUcoLCwkLy+vucORJKshk/sdvr6+aDSaOtfdy8rKOHXqlNH3NbTejJOTE7/+9a9Zvnw5e/bsMXqcmlQqFWlpaVy9epV169Y16l6WoGvXrnTv3p2NGzdy8ODB5g5HkqyCTO53qNVq/P3960zuYHwpAr2YmBg0Gk2Ds/fnnnsOFxeXRjXz0AsMDCQ2NpZt27aZ9KVkaVJSUvDz8yM7O7ve00eSJFWRyb2GwMDA+ypEQtWs2pTuTHr6ejN79uypt9SAp6cnU6ZMYf78+WY5ITJgwACcnJxYvny5VW+uAtjZ2ZGRkcHt27dZtGhRvRvUkiTJ5H6XwMBAKioqau1PGh4eTlFRkcl1YHQ6XYP1ZgB+85vfoNFoeP/9900apyYHBwcGDhxIcXFxvccxrYWXlxdDhw6lqKjIquvYS9KDIJN7DXVtqsLP3ZnqKydQn4CAAIPqzQQEBDBx4kQ+//xzTp8+bdJYNUVFRREaGsratWtt4nH+Ll26EBcXx/fff1/vCSRJaulkcq/B3d0dJyenWtfdg4KCcHJyatRmpyH1ZqCqmcetW7eYNm2ayWPp6c+Ll5eXs2bNmkbfzxI8/PDDBAQEkJOTQ0lJSXOHI0kWSSb3GoQQBAUF1Zrc1Wq1QeUE6mNIvRmoenAqIyOD6dOnc+nSJZPGqsnb25uHHnqIH374gWPHjjX6fs1No9GQkZGBEIJFixZZdS0dSWoqMrnfIzAwkPPnz9e6tm5IOYH61Kw301A5g1dffZUrV64wffp0k8a6V2JiIh4eHixfvtwmNiM9PDwYPnw4p06dIjc3t7nDkSSLI5P7PfSdmWrr5+nh4UGHDh3YsWOHybPFHj16oChKgw2ho6OjSUlJYdq0aVy/ft2ksWqys7MjNTWV8+fPs3nz5kbfzxJ06NCBhx56iO3bt7N3797mDkeSLIpM7vfQJ/e6mk7rdDquX79Ofn6+Sff39PQ0qN4MVLXiO3fuHJ9//rlJY92rffv2dOzYkW+++cYsyz2WoH///gQHB/PVV18Z3e9WkmyZTO73cHBwwMvLq9Z1d6iqvNiqVasG183rY0i9GYA+ffrw0EMP8f777xvdDaouKSkpCCEaXYHSUqjVakaNGoWdnR2ZmZlm++8kSdZOJvda6DdVazuyKIRAp9NRVFRk8pOfhtab0TfzKCwsZP78+SaNdS93d3eSkpI4cOAABQUFZrlnc3Nzc2PkyJGcPXuWFStWNHc4kmQRZHKvRWBgINeuXavzVEx0dDR2dnZs3brVpPsbWm8GYPDgwURGRvL222+bbVaakJCAr68vq1atory83Cz3bG7h4eEkJiaye/dudu3a1dzhSFKzk8m9Fvp197qWZhwcHIiKimLfvn0mb3YaWm9GpVLx9ttvU1BQYJanVqFqKSMtLY3S0lI2bNhglntagqSkJEJDQ1mxYkWDX5qSZOtkcq+Fn58farW6zk1VgLi4OCoqKti9e7dJYxhabwZgyJAhjBo1ijfffJMDBw6YNN69QkJCiI6OZvPmzbWWW7BGKpWK9PR0HBwcyMzMtPpyx5LUGDK510JfIbK245B6fn5+hISEsH37dpObSBhabwbgo48+wsHBgcmTJ5utacXDDz+MVqu1+q5NNbm4uJCenk5JSQnLli2zmc8lScaSyb0OgYGBnDx5st7jinFxcVy8eNHkUsD6ejPbt29vMAn5+/vz/vvv88033zBr1iyTxruXk5MTycnJFBYWGvQFYy1CQ0Pp168f+/bts4mCaZJkCpnc6xAUFFRnhUi9jh074uLi0uhjkefPnzeoINlTTz1FUlISU6dONUtRMaha+w8KCmLNmjWUlZWZ5Z6WoHfv3kRERLBq1SqbqGcvScaSyb0ODW2qQtXyTY8ePTh06JDJBay6dOmCo6Njg0+sQtWa8owZMygrK+PXv/61SePdS19YrKyszKbK6AohGDFiBM7OzmRmZppcqlmSrJVM7nXw8PCos0JkTT169EClUpk8e9fXmykoKGiw3gxUPWX6+uuvk5mZydKlS00a816tW7cmPj6enTt31ruJbG2cnJwYNWoUpaWl5OTkyPV3qUWRyb0OQggCAwMbTO6urq506tSJ3bt3m3wO3diCZFOnTiUqKopnn33WoC8EQ/Tt2xdXV1eWLVtmtg1bSxAcHExycjIFBQVs2bKlucORpAdGJvd6BAYGcu7cuQaP1Ol0Om7cuGFy8Sp9vZkdO3YYVLHR3t6ezz77jJMnT/K73/3OpDHvpdVqSUlJ4cyZMyY/nGWpEhIS6NChA2vWrLGpv0wkqT4yuddD35mpodl7SEgIvr6+DXZZqk9sbKxB9Wb04uPjef7555k+fTrfffedSWPeq1OnTkRERJCXl2e2vwgsgRCCYcOG4ebmRmZmplmqbEqSpTMouQshUoQQ+4UQh4QQr9ZzXboQQhFCxJovxOYTEBAANJzc9fVmTp8+bfLMMCIiwqB6MzX95S9/ISgoiEmTJpnlgR0hBKmpqVRWVrJ69epG38+SODo6kpGRwbVr11iyZIlcf5dsXoPJXQihBv4NpAKdgceEEJ1ruc4VeAGwmYVNR0fHeitE1tS1a1e0Wu0DqTej5+rqyscff0x+fj7vvvuuSePeq1WrVvTu3Zsff/yRw4cPm+WeliIgIIBBgwZx8OBBNm3a1NzhSFKTMmTmHgccUhTliKIo5cB8YFgt1/0/4F3Aps6c6TdVG5rp2dvbEx0dTX5+PlevXjVpLEPrzdSUlpbGmDFjeOutt8zWMLpXr160atWKFStW2FwLu9jYWLp06cK6detsouWgJNXFkOQeCJyo8XPRnd9VE0J0B4IVRVle342EEJOFENuFENvPnTtndLDNITAwkKtXrxq0Bq3T6aisrGTnzp0mjWVMvZmapk2bhrOzM5MmTTLLSReNRkNaWholJSVs3Lix0fezJEIIhgwZQqtWrcjKyjL5i1iSLF2jN1SFECrgQ+Dlhq5VFGWGoiixiqLE+vj4NHboB6Khzkw1eXl5ER4e/sDqzej5+fnx4YcfsmnTJmbMmGHSuPcKCwsjMjKSjRs32lyHI61WS0ZGBjdu3GDx4sU2dfRTkvQMSe7FQHCNn4Pu/E7PFYgE1gshjgEJwFJb2VRt3bo1arXaoHV3qErOV65cMbkRhjH1ZmqaOHEiAwYM4JVXXjE41oYMHDgQe3t7MjMzbabuu56fnx+DBw/m6NGjNlX2WJL0DEnu24B2Qoi2Qgh7YAxQ/WikoiiliqJ4K4oSqihKKPA9MFRRlIafp7cC+gqRhibMdu3a4e7ubpZ6M8asCQsh+OSTTygvL+e5554zeeyaXF1dSU9P59y5c3z11Vc2d8IkOjqabt268c0339jc5rEkNZjcFUWpAJ4DcoGfgIWKovwohHhTCDG0qQO0BIGBgZw6dcqgP99VKhWxsbEcO3bM5Drp+nozxn5BRERE8Kc//YklS5awePFik8a+V3h4eHWFRVt7wlMIweDBg/Hx8WHx4sU2dbZfkgxac1cUZYWiKO0VRQlXFOWtO797Q1GU+4qbKIrS11Zm7XqBgYHcunXL4GTdvXt31Gr1A6s3U9NLL71EdHQ0zz33HJcuXTJp/Hv17t2bjh07snr1agoLC81yT0thb29PRkYGt27dIisrS66/SzZDPqFqAGM2VaGqYJX+1IupDxcZW29Gz87OjpkzZ3LmzBleeeUVk8a+lxCC4cOH06pVKzIzM7ly5YpZ7mspfHx8GDJkCMePH2fdunXNHY4kmYVM7gbw9PQ0qEJkTTqdjvLycpObYHh6etKuXTuD683U1KNHD1588UVmzJhhts1CrVbL6NGjKS8vZ+HChUbHZOmioqLo0aMHmzZtMlsrQ0lqTjK5G8DQCpE1BQYGEhAQ0Kh6Mzqdzqh6MzW9+eabhIaGMnnyZLPVMvfx8WHYsGEUFRWRm5trlntakpSUFFq3bk12djYXLlxo7nAkqVFkcjeQoRUiazKmy1JtTKk3o+fs7Mynn37K/v37efvtt00avzZdunShZ8+ebNu2zaZa80HVXkdGRgYqlYo5c+bIDVbJqsnkbiD9unt9TbPvFRkZadKpFz1T6s3UNHDgQMaNG8c777zDvn37TIqhNsnJyYSGhrJs2TKztfuzFK1atWLs2LGUlZUxZ84cWUFSsloyuRtIn9yNOQ+t0Wjo3r07+/fvp7S01KRxo6Ojja43U9OHH36Iu7s7kyZNMts6uUqlYtSoUTg6OrJgwQKb6r0KVQ+SPfbYY1y8eJG5c+eapeKmJD1oMrkbyNHRkc6dO7Nt2zajZnOxsVUP6hrSI7U2+pM3u3fvNmn27uPjw7Rp0/j+++/5+OOPTYqhNs7Ozjz66KNcvnyZxYsX29wDTqGhoWRkZHDq1CkWLFhgcwXUJNsnk7sRkpKSKC8vZ/PmzQa/x8PDg/bt27Nz506TE0RycjIODg5kZWWZ1Mpv7NixDBo0iNdee40TJ040/AYDBQUFkZqayqFDh1i/fr3Z7mspOnTowLBhwzh69Kg8Ay9ZHZncjeDr60tkZCRbtmzh2rVrBr9Pp9Nx/fp18vPzTRrX2dmZESNGcO7cOZOaaAgh+Pjjj6msrOTZZ5816yy7R48eREdHs2HDBps8QtitWzdSUlIoKCiwyRIMku2Syd1ISUlJVFRUGNXaLiwsDC8vr0bVmwkPD6dnz55s377dpKJkbdu25f/9v//HsmXLyMzMNDmOe+kf4ff392fx4sWUlJSY7d6WIj4+nqSkJHbv3k1ubq5M8JJVkMndSN7e3kRFRbF161aDa4HrT70UFRUZddrmXgMGDMDf35+lS5eadEzv17/+NT169OD55583axK2s7Pj0UcfRaVSsWDBApurIAlVX+pxcXFs2bKFb7/9trnDkaQGyeRugsTERG7fvm1Uq7bo6Gjs7OwaNXtXq9Wkp6dTUVFBdna20WvAGo2GmTNncuHCBaZOnWpyHLXx8PAgPT2ds2fPsmzZMpub3QohSElJoVu3buTl5ZncTlGSHhSZ3E3g5eVFt27d2L59u8F1VhwcHOjatSv79u1r1NlpLy8vUlNTOXbsmEl9QKOjo3n55ZeZNWsWeXl5JsdRG30Fyb1799pk8hNCMHToUDp06MDKlSvZs2dPc4ckSXWSyd1E+tm7MW3odDodFRUV7N69u1FjR0dH06VLF/Ly8gwuZlbTH//4R8LDw5k8ebLZz6j36dOHDh06sHr1ao4fP27We1sC/Rn/0NBQlixZYpObyJJtkMndRJ6enkRHR7Njxw6D17/9/PwICQlh27ZtjTpWJ4TgkUcewc3NjaysLKMfsnFycuLTTz/l0KFDvPnmmybHUVdsw4cPx8PDwyYrSELV8taYMWPw9/cnMzNTNtqWLJJM7o2QmJiIoihGbbDFxcVx6dIlDh061KixHRwcGDlyJKWlpSxfXm9f8loNGDCAJ554gvfff9/sNWIcHBwYPXo0N2/eZNGiRTZXQRKqqmSOHTsWDw8PvvzyS06dOtXcIUnSXWRybwQPDw9iYmLYuXOnweUFOnbsiIuLS6M2VvVCQkJISkpi7969JiXov/3tb7Rq1Yqnn37a7AnY19eXoUOHcvz4cZPO5lsDJycnxo8fj6OjI1988YXNNRKXrJtM7o3Up08fhBAG101Xq9X06NGDQ4cOmeU4Yp8+fQgJCWHFihVG38/Ly4uPPvqI7du389FHHzU6lntFRkaSkJDA1q1b2bt3r9nvbwnc3NwYP348QgjmzJljcg0hSTI3mdwbyd3dne7du7N7926D29r16NEDlUplltm7SqVi5MiRqFQqFi9ebPQMfPTo0QwePJg//OEPTbJ2nJycTJs2bVi6dKlJtXGsgZeXF+PGjePmzZvMmTPHqKeXJampyORuBr179zZq9u7q6kqnTp3YvXu3SbVi7uXu7s6QIUMoLi42+nijvjSBEIJnnnnG7OfT1Wr1XRUkzdU4xNK0bt2axx9/nNLSUr744gub/ZyS9ZDJ3Qzc3NyIjY1l9+7dBi+N6HQ6bty4Ybblis6dOxMTE8OmTZuMbg4SEhLC22+/zapVq/jyyy/NEk9NLi4uZGRkUFpaSnZ2ts094KQXEhLCo48+ytmzZ5k/f75ZvrglyVQyuZtJ7969UavVBs/eQ0JC8PX1bVQbvnulpKTg7e3N4sWLjX5Q6le/+hXx8fG88MILTbIxGBwcTEpKCgcOHDBbX1dL1K5dO0aMGEFhYaHNnhSSrINM7mbi4uJCbGwse/bsMaj/phACnU7H6dOnzVaG197envT0dMrKyli6dKlRXxpqtZrPPvuMS5cu8fLLL5slnnvFxsbStWtX1q9fz8GDB5tkDEsQGRlJWloaBw4cICcnx2b/UpEsm0zuZtSrVy80Gg3ffPONQdd37doVrVZrlo1VvdatW5OcnMz+/fuNbhASFRXFb3/7W2bPns2aNWvMFpOe/uErPz8/Fi9ezMWLF80+hqWIjY2lf//+7N27l5UrV8oELz1wMrmbkYuLCzqdjr1793Lu3LkGr7e3tyc6Opr8/HyDK0waIj4+noiICFavXs3Zs2eNeu/rr79O+/bt+eUvf9kkpz7s7OwYPXo0AAsWLLDpdenevXtXNxO3xWYmkmWTyd3MevXqhb29vcHryjqdjsrKSnbs2GG2GIQQDBs2DK1Wy6JFi4xKoA4ODsyYMYOjR4/ypz/9yWwx1eTp6cnIkSM5c+aMTVaQ1BNC8PDDDxMTE8OGDRv4/vvvmzskqQWRyd3MnJyciIuLY9++fQbNmr28vAgPD2fHjh1m3XxzcXFh+PDhJnVvSkpK4umnn+bDDz9k586dZouppnbt2tG3b1/27Nlj1mUpS6NfiurUqRO5ubmNLhonSYaSyb0JPPTQQ9jb2xu89q7T6bhy5Qr79+83axwREREkJCSwfft2o+/93nvv4evry9NPP91kzaETExNp3749ubm5Zu3tamn0D5qFhYWxdOlSkzppSZKxZHJvAo6OjiQkJJCfn8/p06cbvL5du3a4u7s3yQx2wIABtG7dmpycHKO6N3l6evLPf/6TXbt2MWnSpCY50ieEYMSIEbi7u7Nw4UKz7jtYGo1Gw+jRowkMDGTRokUcOXKkuUOSbJxM7k0kISEBrVZr0OxdpVIRGxvLsWPHjN4AbYhGo6nu3rRkyRKjSg2PGjWKP/7xj/z3v/9l7NixTbL5qa8geePGDTIzM236XLi9vT2PP/44Xl5ezJ8/n+Li4uYOSbJhMrk3EUdHR3r27ElBQYFB5WC7d++OWq1uktm7t7c3KSkpHD161KjG3gB/+tOfeO+991iwYAGjRo1qksfq/fz8qitINsURTEvi6OjIuHHjcHFxYe7cuQadqpIkU8jk3oTi4+NxcHAw6Bick5MTkZGR7Nmzx+jmG4aIiYmhc+fO5OXlGT1jnDp1Kv/6179YunQpQ4cObVSbwLpERUVVN6C21QqSeq6urowfPx61Ws2cOXMMLjgnScaQyb0JOTg48NBDD3HgwAGDEqpOp6O8vNzszTOgan17yJAhuLq6mtS96Ve/+hWzZs1i7dq1pKamNkmHpYEDBxISEsJXX31l9uUpS+Pp6cn48eO5desWc+bMsen9Bql5yOTexOLi4nB0dDRo9h4YGEhAQIBZ683UpO/edOnSJVasWGH0+5988knmzp3Lpk2bSE5ONvsTpvoKklqt1qYrSOr5+voyduxYrly5IitJSmZnUHIXQqQIIfYLIQ4JIV6t5fWXhBD5Qog9Qoi1Qog25g/VOmm1Wh566CEOHTpk0HG/uLg4zp8/b3RlR0OFhISQmJjInj172LNnj9HvHzNmDFlZWezevZt+/fqZfYbt6upKRkYGly5dYuHChTaf8IKCghg9ejTnzp1j3rx5lJeXN3dIko1oMLkLIdTAv4FUoDPwmBCi8z2X7QJiFUXpCiwC3jN3oNYsLi4OJycng2bvXbp0wcnJqUkf7ElMTCQ4OJjly5ebNPseNmwYX331FQcOHCApKYmTJ0+aNb6QkBCGDh1KYWEhn332mc1vOoaHh5Oenk5RUZFcg5fMxpCZexxwSFGUI4qilAPzgWE1L1AUJU9RFP0u2/dAkHnDtG729vb06tWLI0eOUFhYWO+1Go2GmJgY9u/f32Qt2/QP1QghyMrKMun44cCBA1m1ahVFRUUkJiY2+LmM1a1bNyZOnMjNmzeZOXOmzT/407lzZ0aNGsW5c+f45JNP2LdvX3OHJFk5Q5J7IFBzPaHozu/q8hSwsrYXhBCThRDbhRDbbX02di+dToezs7NBs/fY2FgAo6s6GsPDw6O6e5OpRa0SExP5+uuvuXDhAn369DF7Gd+QkBAmT56Mt7c3CxYsYP369TZbhwaqEvwvf/lLfHx8yMrKYunSpXKZRjKZWTdUhRDjgFjg/dpeVxRlhqIosYqixPr4+JhzaItnZ2dH7969OXbsWIO9Sj08PGjfvj07d+5sskf/oWoJKCYmho0bN5q8xh8fH09eXh5lZWUkJiby448/mjVGNzc3nnzySbp168Y333zDggULmuSoqKXw9PTkiSeeoE+fPuzatYsZM2YY9JSzJN3LkOReDATX+Dnozu/uIoRIBn4PDFUUxXb/39cIPXr0wNXV1aAZqE6n4/r16+Tn5zdpTCkpKXh5eZGdnW3y+fXo6Gg2bNiAEIKkpCSzFxvTaDQMGzasupPTzJkzDWqIYq3UajX9+/dnwoQJlJeXM3PmTL7//nub/qtFMj9Dkvs2oJ0Qoq0Qwh4YAyyteYEQIgb4lKrEbtsHlBtBP3svLCxscKYcFhaGl5dXk1dM1HdvunbtmtHdm2rq1KkT3377LS4uLvTv35/NmzebNU4hBPHx8YwfP55r167x2Wef2XQ3J4C2bdsyZcoUwsPDyc3N5csvv2ySGvuSbWowuSuKUgE8B+QCPwELFUX5UQjxphBi6J3L3gdcgEwhxG4hxNI6btfide/eHTc3twZn70IIYmNjKSoqMvtplHv5+/ub3L2ppvDwcDZs2ICPjw8PP/wweXl5ZoyyStu2bZk8eTKenp7MmzePjRs32vSM1snJiTFjxpCamsqRI0f45JNPZNExySAGrbkrirJCUZT2iqKEK4ry1p3fvaEoytI7/05WFMVPUZToO/8bWv8dWy6NRkOfPn04ceIEhw8frvfa6Oho7OzsHki984SEBMLDw03q3lRTSEgIGzZsoE2bNgwePJhVq1aZMcoqHh4e/OIXvyAyMpK1a9eSlZVl0xuPQgji4uKYNGkSDg4OzJkzh6+//tqmi6xJjSefUG0GMTExuLu7Nzh7d3BwoGvXruzbt69J6rnUJIRg+PDhaLVasrKyGlUB0t/fn2+++YZOnToxdOhQsrOzzRhpFTs7O0aOHElycjL5+fnMmjXLpnuyQlWBtcmTJ9O9e3c2bdrE559/bvOfWTKdTO7NQK1Wk5iYSHFxcYPrxjqdjoqKCnbt2tXkcem7N509e7bR1Rm9vb1Zt24dsbGxZGRkMG/ePDNF+TMhBL169eLxxx+ntLSUzz77zOaXLOzs7BgyZAgZGRlcuHCBTz75xOYLrUmmkcm9mXTr1g0PD48GZ+9+fn60adOGzZs3c+bMmSaPS9+9adu2bY3uDOXh4cHq1avp06cP48aNY+bMmWaK8m4RERFMmjQJFxcXvvjiCzZv3mzT6/Dw85l4Pz8/Fi9ezJIlS2z6iKhkPIw2lpQAABjRSURBVJncm4larSYpKYlTp05x4MCBeq9NTU1FpVLxn//854E8qVmze1Njqz+6uLiwYsUKUlJSmDRpEh999JGZorxbq1ateOqpp+jYsSOrV69myZIlTdJcxJJ4eHjwxBNPVNcKmjFjRpNvvkvWQyb3ZtS1a1datWpFXl5eg7P3SZMm4ePjw4IFC/j222+bdGZas3tTdnZ2o8dydHQkOzubESNG8MILL/DOO++YKdK7abVaMjIy6NevH3v27OHzzz9vshIOlkKlUtGvXz8mTpxIRUUF//nPf/juu+9s/i8XqWEyuTcjlUpFYmIiZ86caXBG7urqyhNPPEFUVBTr1q0jOzu7SWemNbs3zZs3z6j+q7XRarUsXLiQsWPH8rvf/Y7XX3+9SRKQEILExETGjBlDSUkJM2bMMHvdG0vUpk0bpkyZQvv27VmzZg1z586VNeJbONFc3/CxsbFKU9ZOsRaVlZVMnz4dtVrNlClTEELUe72iKGzcuJF169YREBDAmDFjcHV1bZLYFEVh27ZtrFmzBrVaTUpKCt26dWswxvrcvn2bKVOmMHPmTH7zm9/wwQcfNOp+9Tl//jzz58/n4sWLDBo0CJ1O12RjWQpFUdixYwe5ublotVqGDx9OREREc4clmZEQYoeiKLENXieTe/Pbu3cvixcvZtSoUXTp0sWg9xQUFLB48eLqBtOBgfXVcmuckpIScnJyOH78OO3atavu6GQqRVF48cUX+eijj/jlL3/J9OnTUama5o/IGzdukJ2dzYEDB4iOjiYtLQ2NRtMkY1mSs2fPkpWVxdmzZ+nZsycDBgxArVY3d1iSGcjkbkUqKyv55JNPUBSFZ555xuBEd+bMmepH0ocNG0ZkZGSTxagoClu2bGHt2rVoNBpSUlLo2rWryTNhRVH4/e9/zzvvvMP48eOZNWtWkyVdRVHIy8vj22+/JTAwkNGjRzfZXzuW5NatW6xevZrt27fj7+9Peno6Xl5ezR2W1EgyuVuZH3/8kUWLFjFy5EiioqIMft+1a9dYuHAhx48fp0+fPvTr169Jlx4uXLhATk4OJ06coEOHDjzyyCO4uLiYfL+33nqLP/zhD4waNYq5c+dib29vxmjvlp+fz5IlS9BqtTz66KMEBwc3/CYbUFBQQE5ODrdv3yYtLa1RX8pS85PJ3cooisInn3zC7du3efbZZ41aprh9+zbLly9n165ddOzYkREjRjRpkqysrGTLli2sW7cOOzs7UlNTiYyMNDlhTJs2jd/85jekpaWxaNEiHBwczBzxz86ePcv8+fMpLS0lLS2N7t27N9lYlqS0tJTs7GwKCwuJiooiLS0NrVbb3GFJJpDJ3Qr99NNPLFy4kOHDh9OtWzej3qtfNlm9ejW+vr6MGTMGDw+PJoq0yvnz58nJyaGoqIiOHTuSlpZm8ix+xowZTJkyhf79+5OTk4Ozs7OZo/1ZWVkZWVlZHD58mNjYWFJSUlrEenRlZSUbN25k/fr1eHh4kJ6e3qR7NVLTkMndCimKwowZM7h58ybPPfecSZuMhw4dYtGiRajVakaPHk1ISEgTRPqzyspKNm/eTF5eHvb29gwePJguXbqYNIv/4osvmDhxIj179mT58uW4u7s3QcRVKisrWbt2Ld999x0hISFkZGQ0annJmhw/fpzFixdz5coV+vXrR69eveQyjRWRyd1K7d+/n/nz5zN06FBiYmJMusf58+f58ssvuXTpEo888ojJ9zHGuXPnyMnJobi4mE6dOpGWlmbS7HvRokU89thjREdHs2rVqibfANy3bx85OTk4OTkxevRoAgICmnQ8S1FWVsayZcvIz88nLCyM4cOHt4hNZlsgk7uVUhSFzz77jLKyMp577jmTlwvKyspYtGgRR44cISEhgYcffrjJjhvqVVZW8t1337F+/Xq0Wm31LN5Yy5cvJz09HW9vb6ZOncrTTz/dpMs0p0+fZv78+Vy9epUhQ4YYvSRmrRRFYdeuXaxcuRKNRkOPHj3o0aMHnp6ezR2aVA+Z3K3YwYMHmTdvHkOGDGnUhl9lZSW5ubls3bqV8PBwRo0a1aSblXpnz55lyZIlnDp1ii5dujB48GCcnJyMusf333/Pq6/+//bOPbat677jnyOKT1GyZMt6xLJkOXH8iCVbDzpFOwzt2szRbKcZujYxumEoAhSrlzTu0K5J/um2P+oURdMGDVC0aVokXjA3SbMo8IJsKdzUVrolpB62/K5sPSyJFEXrwafE19kfEi9EWbYlmRQt8nyAA957dXXu757L+z2/+7u/c/gMf/jDH1i3bh1PP/00Tz75ZNqEJxgM8uabb9LX18eDDz7I5z73uZx54Tg6OsqJEye4dOkSUkq2bNmCzWbj3nvvTbtDoFg6StxXMVJKXnnlFfx+P0899dQdv+xrb2/nvffeo6SkhIMHD65IrnM8Huejjz7iww8/xGw2s2/fPrZv377kev74xz9y5MgRjh8/jtVq5Rvf+Abf+ta3qKysTLnNsViMDz74gI8//hiDwUB9fT179uwhV37M3ev10t7eTkdHB36/n+LiYpqbm2loaFhy56xIH0rcVzk9PT28/vrr7Nu3j+bm217H29LX18cbb7yBlJIvf/nLbN68OQVW3p6RkRFaW1txOp3s3LmTlpaWZQlFd3c3zz//PMeOHUOv1/O1r32N73znO2k5j6GhIex2O2fPniUWi7Fp0yZsNhvbtm3LCU82Fotx8eJFHA4HfX196HQ6HnjgAWw2Gxs2bFAvXzOMEvdVjpRSm9XwqaeeSsnozfHxcY4dO8bo6CgPP/zwis21EovFaGtr4+TJk5jNZvbv38+2bduWVdeVK1f44Q9/yK9//WtisRiPP/44zzzzTFpG5wYCATo7O3E4HExOTlJUVERTUxONjY05k1njdrtxOBycPn2acDhMRUUFNpuNnTt3pnUsheLmKHHPAq5evcrRo0dpaWlhz549Kalzenqat99+m8uXL9PU1ERLS8uK5Xi7XC5aW1txuVzU1dXR0tKC2WxeVl3Dw8P8+Mc/5mc/+xmBQIADBw7w3HPP8alPfSrFVs+EmC5fvozdbufq1avk5eVpnmxVVVVOeLLT09N0d3djt9txu90YjUZ2795Nc3MzpaWlmTYvp1DingVIKXn11Ve5fv063/zmN9Hr9SmpNx6Pc+LECT766CNqamr4yle+smIx1VgsxqlTpzh16hQWi4X9+/ezdevWZdc3NjbGSy+9xIsvvsjY2Bif/exnefbZZ3nooYfSIroejwe73c7p06eZnp7WPNm6urqUXZ+7GSkl165dw263c/78eeLxOLW1tdhsNrZu3ZoTYatMo8Q9S+jr6+PVV19l7969KfdKz5w5w7vvvkthYSEHDx6krKwspfXfCqfTyTvvvIPb7WbXrl3s3bt32V48gN/v5+WXX+ZHP/oRQ0NDNDU18dxzz/Hoo4+mRXDC4TBnzpzRPFmTyURDQwPNzc2sXbs25ce7G/H7/VrYyuv1UlhYqIWtVM58+lDinkW89tpr9Pf3s2PHDpqbm6murk6ZVzo4OMhvfvMbwuEwX/rSl7j//vtTUu9iiMVinDx5klOnTmG1Wtm/f/8dH396epqjR4/ygx/8gJ6eHrZt28Z3v/tdvvrVr6bFs5ZSMjAwwCeffMKFCxeSUgnvu+++nAjZxONx/vSnP2G327ly5Qp5eXls27YNm81GTU1NTrTBSqLEPYvw+Xy0tbVpoYCysjKamprYtWtXSnKxvV4vx44dw+l08oUvfIFPf/rTK3pDDg8P09raitvtZvfu3ezdu/eO8/FjsRhvvfUWR44c4fTp01RXV/Ptb3+bJ554Im0hKJ/PR3t7O+3t7fj9fkpKSrRUwjt5KllNjI2N4XA46OzsZGpqitLSUmw2G/X19SsyxiIXUOKehYTDYc6ePYvD4cDpdKLX66mrq6O5ufmO874jkQitra2cO3eO+vp6Dhw4sKI/ahGNRjl58iRtbW1YrVYeeeSRlPyCkJSS999/n+9///u0tbWxfv16Dh8+zKFDh9I2sVosFuPChQvY7XYGBgbIz8+nrq4Om82Wlvz8u5FIJMK5c+ew2+0MDw+j1+upr6/HZrNRXl6eafNWNUrcs5yhoSEcDgdnz54lGo2yYcMGmpubeeCBB5YdfpBScvLkST788EOqqqp47LHHVjzlb2hoiNbWVkZHR7n//vvZsmULmzdvpqSk5I6fJtra2jhy5AjvvfceRUVFHDp0iMOHD6dVbFwuF3a7ne7ubiKRCFVVVezZs4cdO3bkxEyUcON3dePGjdhsNrZv354Tv4qVarJX3A8fhq6u1Bu0SonF4wT8fnw+H5FIhLy8PKxWK4WFhcsW+UAwiMfjIS8vj7KyMowrnM8spWRichK/308sGgUgPz8fk9mM2WTCZDaju4OXpH6/n4GBAdyjo+QJQUVlJdUbN6Y1bBCPx/H7/Xh9PqKRCDqdDmthIYWFheTniMjHZtvAN9sGeTodhVYrFosFg8GQW7H53bvhJz9Z1r8uVtxVt7nK0eXlUVRURGFREVNTU/h8Prw+H16vF5PZTGFhIRazeUk3ToHFgr6ighG3G5fTicViwWQyYTKZyNfrSfctKISgpLiY4uJiopEIoakppkIhgoEAfp8PAIPBoIm90WQibwnnZ7Va2bFjB5tCIa4NDOB0OnEOD1NWXk51dTUFaYjJ5829TqEQXp+PyYkJJicnsVgsFBYWYjKZ0t62mUSXl8eaoiKKZtvA5/Mx6fUyOTkJQmDQ6zEYjRgNBgxGIwa9PrcEP8WsPs9dcVt8Ph+dnZ20t7fj9XqxWq00NjbS1NREUVHRouvx+/387ne/48qVK/j9fmBGGGtqaqipqWHTpk2Ulpau2A0Yj8cZHh7m6tWrXL16lWvXrhGPx9HpdFRXV1NbW8vmzZuprKxcUvrj4OAgL7zwAj//+c8JBoMcOHCAz3/+8zQ2NrJr164ltdlSGB8f114+hkIh1q9fz/bt26moqKC8vDwloai7HZ/Px+DgIMPDw1qZmpoCQKfTUV5ezj333KOV9evX53wuffaGZRSLJpGi5nA46OnpQQjB1q1baWpq4t577120cEgpGRsbo6+vj/7+fvr6+vDNetAWi0UT+pqaGsrKylZMkMLhMAMDA1y5coXe3l5GRkYAMJlMmtAvJV7v8Xj46U9/yssvv4zT6dS2b9myhYaGBhobG2loaKChoSGlk4klXj46HA6Gh4dJ3JMGg4Hy8nKtVFRUUFZWltXD/qWUTExMMDw8zNDQEE6nk+HhYcLhMDATnqusrKSyslIT/HXr1uWU4CtxVyQxPj5Oe3s7nZ2dBINBSkpKaGpqWtaMf1JKxsfH6e/v18R+cnISALPZnOTZl5eXr5jY+/1+ent7Nc/e6/UCUFxcrIl9bW3touaGdzqddHZ20tHRQWdnJ52dnfT29mp/r6qqShL8xsbGlExFEIlEcLvdjIyM4HK5GBkZYWRkhOnpaW2ftWvXat59QvSLioqy1suXUnL9+nXNs3c6nTidTiKRCDDTCc4X/LVr12ZteyhxVyxINBrlwoULOBwOBgYG0Ol02uCojRs3LvuGmJiYSPLsJyYmgBkver7Yr4SXlXjaSAh9b2+vJpAVFRWaV19dXb3oF8/j4+N0dXVpgt/R0cGlS5eIx+MArFu3Lsm7b2xs5L777rvj85VSMjk5mST2LpeL8fFxbR+TyZQk9uXl5ZSVlWVtNko8Hsfj8SSFc1wuF7FYDACj0ZgUzrnnnntYs2ZNVgi+EnfFbZk/419ZWRnNzc3U19ff8eCoyclJTez7+/sZGxsDZm666upqTeyXGh9fLnPj9b29vQwMDKQkXh8IBDhz5ozm3Xd0dHD27FktjGC1Wtm9e3eSl79jx46UjJadnp7G7XYnif7IyIjm0QohKC0tvUH0rVZrVojcfGKxGKOjo0mCPzIyonW+ZrM5SewTbbHa5gRS4q5YNOFwmO7ubhwOBy6XC4PBoA2OqqioSMkxvF6v5tX39/dz/fp1YOaRer7Yr0T+dyJen/Ds58bra2pqWLNmDQUFBVitVgoKCpKWbycG4XCY8+fPJ4V1urq6CAQC2jnX1dUlCX59fX1KRs4mnljmhnVcLpcWooKZ9yTzwzqlpaVZmXcfjUZxu90MDQ1pIR23281c3dPr9UnXeX6Z+zfzEjPP0oESd8WSkVJqA07OnTunDY4qKyvDYrFgsVgwm803LJtMpiV73z6fLylm7/F4gJkbLSH25eXlWgpmoujTlB4XCAS0eP21a9fw+XxJce65GAyGG276m3UERqMRIQSxWIyenp4kwe/o6NCeaBLzsdTW1lJSUnLLUlxcTElJCRaLZdFtEQqFbojju91uLYwBMx2b2WzWPucvzy+Jv602zzcSieByufB4PAQCgaTi9/sJBAIEg0EW0kYhxA3X3GKx3LRzSEdYLKXiLoR4GHgR0AG/lFI+P+/vRuA1oAm4Djwmpey7VZ1K3O9uQqEQXV1ddHd34/P5CAaD2uPtQiSEfr74L9QZJJbneoqBQCDJs3e73QseRwhxg+CbTCaMRuOC2+eXpQyWiUajSTf8/OX5YrAQOp1uwY4gUYLBIL29vVy8eJHOzk4GBwcZHx9nYjYH/lb3p16vv20HcLNitVqRUuLxeBgZGcHj8RAKhZiamiIUCt2wfCs78vPzF+wE5m+bv240Gu/aLJd4PE4oFFrwWi/UIURnB9vNx2g0Lnjtt27duuypKFIm7kIIHXAZeAgYBOzAQSnl+Tn7HALqpZT/IIR4HPhrKeVjt6pXifvqQkpJOBwmGAwSDAYJhUI3LM//DAaDN/3Sw8wX/2adQX5+PpFIhGg0qpXEeiQS0Uo4HNY+w+HwLY8HM52DwWDQOoO5nUJCgBJPInl5eQgh0Ol0CCGSts1fhplwzPT0NNPT00xNTWllrlAm2mehjjJhW15eHvn5+Vr9ifaPx+PEYrGk9giHw0nHCgaDWmcTiUSIxWJJJRqNEovFkFJqbV1QUIDFYkGv16PT6dDpdOTn599QDAZD0vr8NkhoiZQSKaV2nFuh1+u1unQ6nVZnYjlhz9zlhdbz8/OT1udum/s5f1ui3eeXpWyHmaeBREc4975IXIu51yUUCnHgwAEaGxtv2Ta3+A6nbITqHqBHSnl1tuJjwBeB83P2+SLwL7PLbwEvCSGEzFTMR5FyhBAYjUaMRiMlJSWL/r9IJHLTzmD+ssfjIRgMai8j04GUUhPguXHou4GEbYtFp9NpHWImiMfjt3yaWwyJTnohliofmY6FL4Xjx48vW9wXy2LEfQNwbc76IPDgzfaRUkaFEJPAOsAzdychxNeBr8+u+oUQl5ZjNFA6v+4cQJ1zbqDOOTco/d73vrfcc65ZzE4rmgQrpfwF8Is7rUcI4VjMY0k2oc45N1DnnBusxDkv5m3GELBxznrV7LYF9xFC5ANrmHmxqlAoFIoMsBhxtwNbhBC1QggD8Djw7rx93gX+fnb5b4ATKt6uUCgUmeO2YZnZGPqTwH8zkwr5KynlOSHEvwEOKeW7wCvAUSFEDzDGTAeQTu44tLMKUeecG6hzzg3Sfs4ZG8SkUCgUivRxd44gUCgUCsUdocRdoVAospBVJ+5CiIeFEJeEED1CiGcybU+6EUL8SgjhFkKczbQtK4UQYqMQ4vdCiPNCiHNCiKczbVO6EUKYhBCfCCFOz57zv2bappVACKETQnQKIY5n2paVQAjRJ4ToFkJ0CSHSOkR/VcXcFzMVQrYhhPhzwA+8JqXcmWl7VgIhRCVQKaXsEEIUAu3Ao1l+nQVQIKX0CyH0QBvwtJTy/zJsWloRQvwT0AwUSSn3Z9qedCOE6AOapZRpH7S12jx3bSoEKWUYSEyFkLVIKU8yk4GUM0gpnVLKjtllH3CBmVHQWYucwT+7qp8tq8fzWgZCiCpgH/DLTNuSjaw2cV9oKoSsvulzHSHEJqAB+DizlqSf2RBFF+AGPpBSZvs5/wT4Z+DOJqhZXUjgf4QQ7bPTsaSN1SbuihxCCGEFfgscllLeXbN8pQEpZUxKuZuZUeB7hBBZG4YTQuwH3FLK9kzbssL8mZSyEWgB/nE27JoWVpu4L2YqBEUWMBt3/i3wupTy7Uzbs5JIKSeA3wMPZ9qWNPIZ4JHZGPQx4C+EEP+eWZPSj5RyaPbTDfwnM6HmtLDaxH0xUyEoVjmzLxdfAS5IKV/ItD0rgRBivRCieHbZzEzSwMXMWpU+pJTPSimrpJSbmLmPT0gp/zbDZqUVIUTBbIIAQogC4C+BtGXBrSpxl1JGgcRUCBeAN6SU5zJrVXoRQvwH8L/AViHEoBDiiUzbtAJ8Bvg7Zry5rtnyV5k2Ks1UAr8XQpxhxon5QEqZE+mBOUQ50CaEOA18AvyXlPL9dB1sVaVCKhQKhWJxrCrPXaFQKBSLQ4m7QqFQZCFK3BUKhSILUeKuUCgUWYgSd4VCochClLgrFApFFqLEXaFQKLKQ/wdGr1gR2H62XQAAAABJRU5ErkJggg==\n", - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "%pylab inline\n", - "\n", - "testmus = np.linspace(0, 5, 11)\n", - "results = []\n", - "for mu in testmus:\n", - " results.append(runOnePoint(mu, data, pdf, init_pars, par_bounds))\n", - "\n", - "obs = [1.0 / r[-2:][0] for r in results]\n", - "exp = [[1.0 / r[-2:][1][i] for r in results] for i in range(5)]\n", - "\n", - "\n", - "def plot_results(testmus, cls_obs, cls_exp, test_size=0.05):\n", - " plt.plot(testmus, cls_obs, c='k')\n", - " for i, c in zip(range(5), ['grey', 'grey', 'grey', 'grey', 'grey']):\n", - " plt.plot(testmus, cls_exp[i], c=c)\n", - " plt.plot(testmus, [test_size] * len(testmus), c='r')\n", - " plt.ylim(0, 1)\n", - "\n", - "\n", - "plot_results(testmus, obs, exp)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 2", - "language": "python", - "name": "python2" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 2 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.14" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/docs/faq.rst b/docs/faq.rst index 48339ebb8e..cf3d146f62 100644 --- a/docs/faq.rst +++ b/docs/faq.rst @@ -169,7 +169,7 @@ invitation for community contributions and new developers. Troubleshooting --------------- -- :code:`import torch` or :code:`import pyhf` causes a :code:`Segmentation fault (core dumped)` +- :code:`import pyhf` causes a :code:`Segmentation fault (core dumped)` This is may be the result of a conflict with the NVIDIA drivers that you have installed on your machine. Try uninstalling and completely removing diff --git a/docs/installation.rst b/docs/installation.rst index f856f2be07..10aa45b5f4 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -27,13 +27,6 @@ Install latest stable release from `PyPI `__... python -m pip install pyhf -... with PyTorch backend -++++++++++++++++++++++++ - -.. code-block:: console - - python -m pip install 'pyhf[torch]' - ... with JAX backend ++++++++++++++++++++ @@ -67,13 +60,6 @@ Install latest development version from `GitHub =1.10.0", # c.f. PR #1657 - "numpy<2.0" # c.f. https://github.com/pytorch/pytorch/issues/157973 -] jax = [ "jax>=0.4.1", # c.f. PR #2079 "jaxlib>=0.4.1", # c.f. PR #2079 @@ -82,7 +77,7 @@ contrib = [ "matplotlib>=3.0.0", "requests>=2.22.0", ] -backends = ["pyhf[torch,jax,minuit]"] +backends = ["pyhf[jax,minuit]"] all = ["pyhf[backends,xmlio,contrib,shellcomplete]"] # Developer extras @@ -182,18 +177,12 @@ markers = [ "fail_jax", "fail_numpy", "fail_numpy_minuit", - "fail_pytorch", - "fail_pytorch64", "only_jax", "only_numpy", "only_numpy_minuit", - "only_pytorch", - "only_pytorch64", "skip_jax", "skip_numpy", "skip_numpy_minuit", - "skip_pytorch", - "skip_pytorch64", ] filterwarnings = [ "error", @@ -202,8 +191,6 @@ filterwarnings = [ 'ignore: Exception ignored in:pytest.PytestUnraisableExceptionWarning', #FIXME: Exception ignored in: <_io.FileIO [closed]> 'ignore:invalid value encountered in (true_)?divide:RuntimeWarning', #FIXME 'ignore:invalid value encountered in add:RuntimeWarning', #FIXME - "ignore:In future, it will be an error for 'np.bool_' scalars to be interpreted as an index:DeprecationWarning", #FIXME: tests/test_tensor.py::test_pdf_eval[pytorch] - 'ignore:Creating a tensor from a list of numpy.ndarrays is extremely slow. Please consider converting the list to a single numpy.ndarray with:UserWarning', #FIXME: tests/test_optim.py::test_minimize[no_grad-scipy-pytorch-no_stitch] 'ignore:divide by zero encountered in (true_)?divide:RuntimeWarning', #FIXME: pytest tests/test_tensor.py::test_pdf_calculations[numpy] 'ignore:[A-Z]+ is deprecated and will be removed in Pillow 10:DeprecationWarning', # keras "ignore:ml_dtypes.float8_e4m3b11 is deprecated.", #FIXME: Can remove when jaxlib>=0.4.12 @@ -213,8 +200,6 @@ filterwarnings = [ "ignore:'MultiCommand' is deprecated and will be removed in Click 9.0. Use 'Group' instead.:DeprecationWarning", # Click "ignore:Jupyter is migrating its paths to use standard platformdirs:DeprecationWarning", # papermill "ignore:datetime.datetime.utcnow\\(\\) is deprecated:DeprecationWarning", # papermill - "ignore:In future, it will be an error for 'np.bool' scalars to be interpreted as an index:DeprecationWarning", # PyTorch - "ignore:__array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments.:DeprecationWarning", # PyTorch interacting with NumPy ] [tool.coverage.run] @@ -243,7 +228,6 @@ module = [ 'jax.*', 'matplotlib.*', 'scipy.*', - 'torch.*', 'uproot.*', ] ignore_missing_imports = true @@ -273,7 +257,6 @@ module = [ 'pyhf.tensor.common.*', 'pyhf.tensor', 'pyhf.tensor.jax_backend.*', - 'pyhf.tensor.pytorch_backend.*', ] ignore_errors = true diff --git a/src/pyhf/cli/infer.py b/src/pyhf/cli/infer.py index 2212cb1e78..9722c52d16 100644 --- a/src/pyhf/cli/infer.py +++ b/src/pyhf/cli/infer.py @@ -36,7 +36,7 @@ def cli(): ) @click.option( "--backend", - type=click.Choice(["numpy", "pytorch", "jax", "np", "torch"]), + type=click.Choice(["numpy", "jax", "np"]), help="The tensor backend used for the calculation.", default="numpy", ) @@ -81,9 +81,7 @@ def fit( } """ # set the backend if not NumPy - if backend in ["pytorch", "torch"]: - set_backend("pytorch", precision="64b") - elif backend in ["jax"]: + if backend in ["jax"]: set_backend("jax") tensorlib, _ = get_backend() @@ -148,7 +146,7 @@ def fit( ) @click.option( '--backend', - type=click.Choice(['numpy', 'pytorch', 'jax', 'np', 'torch']), + type=click.Choice(['numpy', 'jax', 'np']), help='The tensor backend used for the calculation.', default='numpy', ) @@ -211,9 +209,7 @@ def cls( ) # set the backend if not NumPy - if backend in ['pytorch', 'torch']: - set_backend("pytorch", precision="64b") - elif backend in ['jax']: + if backend in ["jax"]: set_backend("jax") tensorlib, _ = get_backend() diff --git a/src/pyhf/infer/calculators.py b/src/pyhf/infer/calculators.py index f9922c6d89..12bb6d5bce 100644 --- a/src/pyhf/infer/calculators.py +++ b/src/pyhf/infer/calculators.py @@ -121,7 +121,7 @@ def cdf(self, value): >>> pyhf.set_backend("numpy") >>> bkg_dist = pyhf.infer.calculators.AsymptoticTestStatDistribution(0.0) >>> bkg_dist.cdf(0.0) - 0.5 + np.float64(0.5) Args: value (:obj:`float`): The test statistic value. @@ -620,7 +620,7 @@ def expected_value(self, nsigma): >>> samples = normal.sample((100,)) >>> dist = pyhf.infer.calculators.EmpiricalDistribution(samples) >>> dist.expected_value(nsigma=1) - 6.15094381... + np.float64(6.15094381...) >>> import pyhf >>> import numpy.random as random diff --git a/src/pyhf/optimize/common.py b/src/pyhf/optimize/common.py index 1f95805d6c..36692d1a80 100644 --- a/src/pyhf/optimize/common.py +++ b/src/pyhf/optimize/common.py @@ -42,11 +42,6 @@ def _get_tensor_shim(): return numpy_shim - if tensorlib.name == 'pytorch': - from pyhf.optimize.opt_pytorch import wrap_objective as pytorch_shim - - return pytorch_shim - if tensorlib.name == 'jax': from pyhf.optimize.opt_jax import wrap_objective as jax_shim diff --git a/src/pyhf/optimize/opt_pytorch.py b/src/pyhf/optimize/opt_pytorch.py deleted file mode 100644 index 55b204ff41..0000000000 --- a/src/pyhf/optimize/opt_pytorch.py +++ /dev/null @@ -1,42 +0,0 @@ -"""PyTorch Backend Function Shim.""" - -from pyhf import get_backend -import torch - - -def wrap_objective(objective, data, pdf, stitch_pars, do_grad=False, jit_pieces=None): - """ - Wrap the objective function for the minimization. - - Args: - objective (:obj:`func`): objective function - data (:obj:`list`): observed data - pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json - stitch_pars (:obj:`func`): callable that stitches parameters, see :func:`pyhf.optimize.common.shim`. - do_grad (:obj:`bool`): enable autodifferentiation mode. Default is off. - - Returns: - objective_and_grad (:obj:`func`): tensor backend wrapped objective,gradient pair - """ - - tensorlib, _ = get_backend() - - if do_grad: - - def func(pars): - pars = tensorlib.astensor(pars) - pars.requires_grad = True - constrained_pars = stitch_pars(pars) - constr_nll = objective(constrained_pars, data, pdf) - grad = torch.autograd.grad(constr_nll, pars)[0] - return constr_nll.detach().numpy()[0], grad - - else: - - def func(pars): - pars = tensorlib.astensor(pars) - constrained_pars = stitch_pars(pars) - constr_nll = objective(constrained_pars, data, pdf) - return constr_nll[0] - - return func diff --git a/src/pyhf/probability.py b/src/pyhf/probability.py index 0cc0330272..0a37d55cf8 100644 --- a/src/pyhf/probability.py +++ b/src/pyhf/probability.py @@ -171,10 +171,10 @@ def log_prob(self, value): >>> independent = pyhf.probability.Independent(poissons) >>> values = pyhf.tensorlib.astensor([8.0, 9.0]) >>> independent.log_prob(values) - -4.26248380... + np.float64(-4.26248380...) >>> broadcast_value = pyhf.tensorlib.astensor([11.0]) >>> independent.log_prob(broadcast_value) - -4.34774364... + np.float64(-4.34774364...) Args: value (:obj:`tensor` or :obj:`float`): The value at which to evaluate the distribution diff --git a/src/pyhf/tensor/__init__.py b/src/pyhf/tensor/__init__.py index 832e7ce0f8..0caa994448 100644 --- a/src/pyhf/tensor/__init__.py +++ b/src/pyhf/tensor/__init__.py @@ -7,7 +7,6 @@ class _BackendRetriever: "_array_types", "jax_backend", "numpy_backend", - "pytorch_backend", ] def __init__(self): @@ -39,21 +38,6 @@ def __getattr__(self, name): "There was a problem importing JAX. The jax backend cannot be used.", e, ) - elif name == 'pytorch_backend': - try: - from pyhf.tensor.pytorch_backend import pytorch_backend - - assert pytorch_backend - # for autocomplete and dir() calls - self.pytorch_backend = pytorch_backend - self._array_types.add(pytorch_backend.array_type) - self._array_subtypes.add(pytorch_backend.array_subtype) - return pytorch_backend - except ImportError as e: - raise exceptions.ImportBackendError( - "There was a problem importing PyTorch. The pytorch backend cannot be used.", - e, - ) @property def array_types(self): diff --git a/src/pyhf/tensor/manager.py b/src/pyhf/tensor/manager.py index b9bf33fba5..a4a6ec3541 100644 --- a/src/pyhf/tensor/manager.py +++ b/src/pyhf/tensor/manager.py @@ -65,11 +65,6 @@ def set_backend( Example: >>> import pyhf - >>> pyhf.set_backend(b"pytorch", precision="32b") - >>> pyhf.tensorlib.name - 'pytorch' - >>> pyhf.tensorlib.precision - '32b' >>> pyhf.set_backend(pyhf.tensor.numpy_backend()) >>> pyhf.tensorlib.name 'numpy' @@ -77,7 +72,7 @@ def set_backend( '64b' Args: - backend (:obj:`str` or :obj:`bytes` or `pyhf.tensor` backend): One of the supported pyhf backends: NumPy, PyTorch, and JAX + backend (:obj:`str` or :obj:`bytes` or `pyhf.tensor` backend): One of the supported pyhf backends: NumPy and JAX custom_optimizer (:obj:`str` or :obj:`bytes` or `pyhf.optimize` optimizer or :obj:`None`): Optional custom optimizer defined by the user precision (:obj:`str` or :obj:`bytes` or :obj:`None`): Floating point precision to use in the backend: ``64b`` or ``32b``. Default is backend dependent. default (:obj:`bool`): Set the backend as the default backend additionally @@ -112,7 +107,7 @@ def set_backend( )(**backend_kwargs) except TypeError: raise exceptions.InvalidBackend( - f"The backend provided is not supported: {backend:s}. Select from one of the supported backends: numpy, pytorch" + f"The backend provided is not supported: {backend:s}. Select from one of the supported backends: numpy, jax" ) else: new_backend = backend diff --git a/src/pyhf/tensor/numpy_backend.py b/src/pyhf/tensor/numpy_backend.py index e63251f193..e843330bb3 100644 --- a/src/pyhf/tensor/numpy_backend.py +++ b/src/pyhf/tensor/numpy_backend.py @@ -324,7 +324,7 @@ def percentile( >>> pyhf.set_backend("numpy") >>> a = pyhf.tensorlib.astensor([[10, 7, 4], [3, 2, 1]]) >>> pyhf.tensorlib.percentile(a, 50) - 3.5 + np.float64(3.5) >>> pyhf.tensorlib.percentile(a, 50, axis=1) array([7., 2.]) @@ -388,7 +388,7 @@ def simple_broadcast(self, *args: Sequence[Tensor[T]]) -> Sequence[Tensor[T]]: ... pyhf.tensorlib.astensor([1]), ... pyhf.tensorlib.astensor([2, 3, 4]), ... pyhf.tensorlib.astensor([5, 6, 7])) - [array([1., 1., 1.]), array([2., 3., 4.]), array([5., 6., 7.])] + (array([1., 1., 1.]), array([2., 3., 4.]), array([5., 6., 7.])) Args: args (Array of Tensors): Sequence of arrays @@ -471,7 +471,7 @@ def poisson(self, n: Tensor[T], lam: Tensor[T]) -> ArrayLike: >>> import pyhf >>> pyhf.set_backend("numpy") >>> pyhf.tensorlib.poisson(5., 6.) - 0.16062314... + np.float64(0.16062314...) >>> values = pyhf.tensorlib.astensor([5., 9.]) >>> rates = pyhf.tensorlib.astensor([6., 8.]) >>> pyhf.tensorlib.poisson(values, rates) @@ -514,7 +514,7 @@ def normal(self, x: Tensor[T], mu: Tensor[T], sigma: Tensor[T]) -> ArrayLike: >>> import pyhf >>> pyhf.set_backend("numpy") >>> pyhf.tensorlib.normal(0.5, 0., 1.) - 0.35206532... + np.float64(0.35206532...) >>> values = pyhf.tensorlib.astensor([0.5, 2.0]) >>> means = pyhf.tensorlib.astensor([0., 2.3]) >>> sigmas = pyhf.tensorlib.astensor([1., 0.8]) @@ -542,7 +542,7 @@ def normal_cdf( >>> import pyhf >>> pyhf.set_backend("numpy") >>> pyhf.tensorlib.normal_cdf(0.8) - 0.78814460... + np.float64(0.78814460...) >>> values = pyhf.tensorlib.astensor([0.8, 2.0]) >>> pyhf.tensorlib.normal_cdf(values) array([0.7881446 , 0.97724987]) diff --git a/src/pyhf/tensor/pytorch_backend.py b/src/pyhf/tensor/pytorch_backend.py deleted file mode 100644 index 0f9f253fcd..0000000000 --- a/src/pyhf/tensor/pytorch_backend.py +++ /dev/null @@ -1,628 +0,0 @@ -"""PyTorch Tensor Library Module.""" - -import torch -import torch.autograd -from torch.distributions.utils import broadcast_all -import logging -import math - -log = logging.getLogger(__name__) - - -class pytorch_backend: - """PyTorch backend for pyhf""" - - __slots__ = ['default_do_grad', 'dtypemap', 'name', 'precision'] - - #: The array type for pytorch - array_type = torch.Tensor - - #: The array content type for pytorch - array_subtype = torch.Tensor - - def __init__(self, **kwargs): - self.name = 'pytorch' - self.precision = kwargs.get('precision', '64b') - self.dtypemap = { - 'float': torch.float64 if self.precision == '64b' else torch.float32, - 'int': torch.int64 if self.precision == '64b' else torch.int32, - 'bool': torch.bool, - } - self.default_do_grad = True - - def _setup(self): - """ - Run any global setups for the pytorch lib. - """ - torch.set_default_dtype(self.dtypemap["float"]) - - def clip(self, tensor_in, min_value, max_value): - """ - Clips (limits) the tensor values to be within a specified min and max. - - Example: - - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> a = pyhf.tensorlib.astensor([-2, -1, 0, 1, 2]) - >>> pyhf.tensorlib.clip(a, -1, 1) - tensor([-1., -1., 0., 1., 1.]) - - Args: - tensor_in (:obj:`tensor`): The input tensor object - min_value (:obj:`scalar` or :obj:`tensor` or :obj:`None`): The minimum value to be clipped to - max_value (:obj:`scalar` or :obj:`tensor` or :obj:`None`): The maximum value to be clipped to - - Returns: - PyTorch tensor: A clipped `tensor` - """ - return torch.clamp(tensor_in, min_value, max_value) - - def erf(self, tensor_in): - """ - The error function of complex argument. - - Example: - - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> a = pyhf.tensorlib.astensor([-2., -1., 0., 1., 2.]) - >>> pyhf.tensorlib.erf(a) - tensor([-0.9953, -0.8427, 0.0000, 0.8427, 0.9953]) - - Args: - tensor_in (:obj:`tensor`): The input tensor object - - Returns: - PyTorch Tensor: The values of the error function at the given points. - """ - return torch.erf(tensor_in) - - def erfinv(self, tensor_in): - """ - The inverse of the error function of complex argument. - - Example: - - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> a = pyhf.tensorlib.astensor([-2., -1., 0., 1., 2.]) - >>> pyhf.tensorlib.erfinv(pyhf.tensorlib.erf(a)) - tensor([-2.0000, -1.0000, 0.0000, 1.0000, 2.0000]) - - Args: - tensor_in (:obj:`tensor`): The input tensor object - - Returns: - PyTorch Tensor: The values of the inverse of the error function at the given points. - """ - return torch.erfinv(tensor_in) - - def conditional(self, predicate, true_callable, false_callable): - """ - Runs a callable conditional on the boolean value of the evaluation of a predicate - - Example: - - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> tensorlib = pyhf.tensorlib - >>> a = tensorlib.astensor([4]) - >>> b = tensorlib.astensor([5]) - >>> tensorlib.conditional((a < b)[0], lambda: a + b, lambda: a - b) - tensor([9.]) - - Args: - predicate (:obj:`scalar`): The logical condition that determines which callable to evaluate - true_callable (:obj:`callable`): The callable that is evaluated when the :code:`predicate` evaluates to :code:`true` - false_callable (:obj:`callable`): The callable that is evaluated when the :code:`predicate` evaluates to :code:`false` - - Returns: - PyTorch Tensor: The output of the callable that was evaluated - """ - return true_callable() if predicate else false_callable() - - def tolist(self, tensor_in): - try: - return tensor_in.data.numpy().tolist() - except AttributeError: - if isinstance(tensor_in, list): - return tensor_in - raise - - def tile(self, tensor_in, repeats): - """ - Repeat tensor data along a specific dimension - - Example: - - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> a = pyhf.tensorlib.astensor([[1.0], [2.0]]) - >>> pyhf.tensorlib.tile(a, (1, 2)) - tensor([[1., 1.], - [2., 2.]]) - - Args: - tensor_in (:obj:`tensor`): The tensor to be repeated - repeats (:obj:`tensor`): The tuple of multipliers for each dimension - - Returns: - PyTorch tensor: The tensor with repeated axes - """ - return tensor_in.tile(repeats) - - def outer(self, tensor_in_1, tensor_in_2): - """ - Outer product of the input tensors. - - Example: - - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> a = pyhf.tensorlib.astensor([1.0, 2.0, 3.0]) - >>> b = pyhf.tensorlib.astensor([1.0, 2.0, 3.0, 4.0]) - >>> pyhf.tensorlib.outer(a, b) - tensor([[ 1., 2., 3., 4.], - [ 2., 4., 6., 8.], - [ 3., 6., 9., 12.]]) - - Args: - tensor_in_1 (:obj:`tensor`): 1-D input tensor. - tensor_in_2 (:obj:`tensor`): 1-D input tensor. - - Returns: - PyTorch tensor: The outer product. - """ - return torch.outer(tensor_in_1, tensor_in_2) - - def astensor(self, tensor_in, dtype='float'): - """ - Convert to a PyTorch Tensor. - - Example: - - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> tensor = pyhf.tensorlib.astensor([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) - >>> tensor - tensor([[1., 2., 3.], - [4., 5., 6.]]) - >>> type(tensor) - - - Args: - tensor_in (Number or Tensor): Tensor object - - Returns: - torch.Tensor: A multi-dimensional matrix containing elements of a single data type. - """ - try: - dtype = self.dtypemap[dtype] - except KeyError: - log.error( - 'Invalid dtype: dtype must be float, int, or bool.', exc_info=True - ) - raise - - return torch.as_tensor(tensor_in, dtype=dtype) - - def gather(self, tensor, indices): - return tensor[indices.type(torch.LongTensor)] - - def boolean_mask(self, tensor, mask): - return torch.masked_select(tensor, mask) - - def reshape(self, tensor, newshape): - return torch.reshape(tensor, newshape) - - def shape(self, tensor): - return tuple(map(int, tensor.shape)) - - def ravel(self, tensor): - """ - Return a flattened view of the tensor, not a copy. - - Example: - - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> tensor = pyhf.tensorlib.astensor([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) - >>> pyhf.tensorlib.ravel(tensor) - tensor([1., 2., 3., 4., 5., 6.]) - - Args: - tensor (Tensor): Tensor object - - Returns: - `torch.Tensor`: A flattened array. - """ - return torch.ravel(tensor) - - def sum(self, tensor_in, axis=None): - return ( - torch.sum(tensor_in) - if (axis is None or tensor_in.shape == torch.Size([])) - else torch.sum(tensor_in, axis) - ) - - def product(self, tensor_in, axis=None): - return torch.prod(tensor_in) if axis is None else torch.prod(tensor_in, axis) - - def abs(self, tensor): - return torch.abs(tensor) - - def ones(self, shape, dtype="float"): - try: - dtype = self.dtypemap[dtype] - except KeyError: - log.error( - f"Invalid dtype: dtype must be one of {list(self.dtypemap)}.", - exc_info=True, - ) - raise - - return torch.ones(shape, dtype=dtype) - - def zeros(self, shape, dtype="float"): - try: - dtype = self.dtypemap[dtype] - except KeyError: - log.error( - f"Invalid dtype: dtype must be one of {list(self.dtypemap)}.", - exc_info=True, - ) - raise - - return torch.zeros(shape, dtype=dtype) - - def power(self, tensor_in_1, tensor_in_2): - return torch.pow(tensor_in_1, tensor_in_2) - - def sqrt(self, tensor_in): - return torch.sqrt(tensor_in) - - def divide(self, tensor_in_1, tensor_in_2): - return torch.div(tensor_in_1, tensor_in_2) - - def log(self, tensor_in): - return torch.log(tensor_in) - - def exp(self, tensor_in): - return torch.exp(tensor_in) - - def percentile(self, tensor_in, q, axis=None, interpolation="linear"): - r""" - Compute the :math:`q`-th percentile of the tensor along the specified axis. - - Example: - - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> a = pyhf.tensorlib.astensor([[10, 7, 4], [3, 2, 1]]) - >>> pyhf.tensorlib.percentile(a, 50) - tensor(3.5000) - >>> pyhf.tensorlib.percentile(a, 50, axis=1) - tensor([7., 2.]) - - Args: - tensor_in (`tensor`): The tensor containing the data - q (:obj:`float` or `tensor`): The :math:`q`-th percentile to compute - axis (`number` or `tensor`): The dimensions along which to compute - interpolation (:obj:`str`): The interpolation method to use when the - desired percentile lies between two data points ``i < j``: - - - ``'linear'``: ``i + (j - i) * fraction``, where ``fraction`` is the - fractional part of the index surrounded by ``i`` and ``j``. - - - ``'lower'``: Not yet implemented in PyTorch. - - - ``'higher'``: Not yet implemented in PyTorch. - - - ``'midpoint'``: Not yet implemented in PyTorch. - - - ``'nearest'``: Not yet implemented in PyTorch. - - Returns: - PyTorch tensor: The value of the :math:`q`-th percentile of the tensor along the specified axis. - - .. versionadded:: 0.7.0 - """ - # Interpolation options not yet supported - # c.f. https://github.com/pytorch/pytorch/pull/49267 - # c.f. https://github.com/pytorch/pytorch/pull/59397 - return torch.quantile(tensor_in, q / 100, dim=axis) - - def stack(self, sequence, axis=0): - return torch.stack(sequence, dim=axis) - - def where(self, mask, tensor_in_1, tensor_in_2): - return torch.where(mask, tensor_in_1, tensor_in_2) - - def concatenate(self, sequence, axis=0): - """ - Join a sequence of arrays along an existing axis. - - Args: - sequence: sequence of tensors - axis: dimension along which to concatenate - - Returns: - output: the concatenated tensor - - """ - return torch.cat(sequence, dim=axis) - - def isfinite(self, tensor): - return torch.isfinite(tensor) - - def simple_broadcast(self, *args): - """ - Broadcast a sequence of 1 dimensional arrays. - - Example: - - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> pyhf.tensorlib.simple_broadcast( - ... pyhf.tensorlib.astensor([1]), - ... pyhf.tensorlib.astensor([2, 3, 4]), - ... pyhf.tensorlib.astensor([5, 6, 7])) - [tensor([1., 1., 1.]), tensor([2., 3., 4.]), tensor([5., 6., 7.])] - - Args: - args (Array of Tensors): Sequence of arrays - - Returns: - list of Tensors: The sequence broadcast together. - """ - - args = [arg.view(1) if not self.shape(arg) else arg for arg in args] - max_dim = max(map(len, args)) - try: - assert not [arg for arg in args if 1 < len(arg) < max_dim] - except AssertionError: - log.error( - 'ERROR: The arguments must be of compatible size: 1 or %i', - max_dim, - exc_info=True, - ) - raise - - broadcast = [arg if len(arg) > 1 else arg.expand(max_dim) for arg in args] - return broadcast - - def einsum(self, subscripts, *operands): - """ - This function provides a way of computing multilinear expressions (i.e. - sums of products) using the Einstein summation convention. - - Args: - subscripts: str, specifies the subscripts for summation - operands: list of array_like, these are the tensors for the operation - - Returns: - tensor: the calculation based on the Einstein summation convention - """ - return torch.einsum(subscripts, operands) - - def poisson_logpdf(self, n, lam): - # validate_args=True disallows continuous approximation - return torch.distributions.Poisson(lam, validate_args=False).log_prob(n) - - def poisson(self, n, lam): - r""" - The continuous approximation, using :math:`n! = \Gamma\left(n+1\right)`, - to the probability mass function of the Poisson distribution evaluated - at :code:`n` given the parameter :code:`lam`. - - .. note:: - - Though the p.m.f of the Poisson distribution is not defined for - :math:`\lambda = 0`, the limit as :math:`\lambda \to 0` is still - defined, which gives a degenerate p.m.f. of - - .. math:: - - \lim_{\lambda \to 0} \,\mathrm{Pois}(n | \lambda) = - \left\{\begin{array}{ll} - 1, & n = 0,\\ - 0, & n > 0 - \end{array}\right. - - Example: - - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> pyhf.tensorlib.poisson(5., 6.) - tensor(0.1606) - >>> values = pyhf.tensorlib.astensor([5., 9.]) - >>> rates = pyhf.tensorlib.astensor([6., 8.]) - >>> pyhf.tensorlib.poisson(values, rates) - tensor([0.1606, 0.1241]) - - Args: - n (:obj:`tensor` or :obj:`float`): The value at which to evaluate the approximation to the Poisson distribution p.m.f. - (the observed number of events) - lam (:obj:`tensor` or :obj:`float`): The mean of the Poisson distribution p.m.f. - (the expected number of events) - - Returns: - PyTorch FloatTensor: Value of the continuous approximation to Poisson(n|lam) - """ - # validate_args=True disallows continuous approximation - return torch.exp( - torch.distributions.Poisson(lam, validate_args=False).log_prob(n) - ) - - def normal_logpdf(self, x, mu, sigma): - x = self.astensor(x) - mu = self.astensor(mu) - sigma = self.astensor(sigma) - - normal = torch.distributions.Normal(mu, sigma) - return normal.log_prob(x) - - def normal(self, x, mu, sigma): - r""" - The probability density function of the Normal distribution evaluated - at :code:`x` given parameters of mean of :code:`mu` and standard deviation - of :code:`sigma`. - - Example: - - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> pyhf.tensorlib.normal(0.5, 0., 1.) - tensor(0.3521) - >>> values = pyhf.tensorlib.astensor([0.5, 2.0]) - >>> means = pyhf.tensorlib.astensor([0., 2.3]) - >>> sigmas = pyhf.tensorlib.astensor([1., 0.8]) - >>> pyhf.tensorlib.normal(values, means, sigmas) - tensor([0.3521, 0.4648]) - - Args: - x (:obj:`tensor` or :obj:`float`): The value at which to evaluate the Normal distribution p.d.f. - mu (:obj:`tensor` or :obj:`float`): The mean of the Normal distribution - sigma (:obj:`tensor` or :obj:`float`): The standard deviation of the Normal distribution - - Returns: - PyTorch FloatTensor: Value of Normal(x|mu, sigma) - """ - x = self.astensor(x) - mu = self.astensor(mu) - sigma = self.astensor(sigma) - - normal = torch.distributions.Normal(mu, sigma) - return self.exp(normal.log_prob(x)) - - def normal_cdf(self, x, mu=0.0, sigma=1.0): - """ - The cumulative distribution function for the Normal distribution - - Example: - - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> pyhf.tensorlib.normal_cdf(0.8) - tensor(0.7881) - >>> values = pyhf.tensorlib.astensor([0.8, 2.0]) - >>> pyhf.tensorlib.normal_cdf(values) - tensor([0.7881, 0.9772]) - - Args: - x (:obj:`tensor` or :obj:`float`): The observed value of the random variable to evaluate the CDF for - mu (:obj:`tensor` or :obj:`float`): The mean of the Normal distribution - sigma (:obj:`tensor` or :obj:`float`): The standard deviation of the Normal distribution - - Returns: - PyTorch FloatTensor: The CDF - """ - # the implementation of torch.Normal.cdf uses torch.erf: - # 0.5 * (1 + torch.erf((value - self.loc) * self.scale.reciprocal() / math.sqrt(2))) - # (see https://github.com/pytorch/pytorch/blob/3bbedb34b9b316729a27e793d94488b574e1577a/torch/distributions/normal.py#L78-L81) - # we get a more numerically stable variant for low p-values/high significances using erfc(x) := 1 - erf(x) - # since erf(-x) = -erf(x) we can replace - # 1 + erf(x) = 1 - erf(-x) = 1 - (1 - erfc(-x)) = erfc(-x) - mu, sigma = broadcast_all(mu, sigma) - return 0.5 * torch.erfc(-((x - mu) * sigma.reciprocal() / math.sqrt(2))) - - def poisson_dist(self, rate): - r""" - The Poisson distribution with rate parameter :code:`rate`. - - Example: - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> rates = pyhf.tensorlib.astensor([5, 8]) - >>> values = pyhf.tensorlib.astensor([4, 9]) - >>> poissons = pyhf.tensorlib.poisson_dist(rates) - >>> poissons.log_prob(values) - tensor([-1.7403, -2.0869]) - - Args: - rate (:obj:`tensor` or :obj:`float`): The mean of the Poisson distribution (the expected number of events) - - Returns: - PyTorch Poisson distribution: The Poisson distribution class - - """ - # validate_args=True disallows continuous approximation - return torch.distributions.Poisson(rate, validate_args=False) - - def normal_dist(self, mu, sigma): - r""" - The Normal distribution with mean :code:`mu` and standard deviation :code:`sigma`. - - Example: - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> means = pyhf.tensorlib.astensor([5, 8]) - >>> stds = pyhf.tensorlib.astensor([1, 0.5]) - >>> values = pyhf.tensorlib.astensor([4, 9]) - >>> normals = pyhf.tensorlib.normal_dist(means, stds) - >>> normals.log_prob(values) - tensor([-1.4189, -2.2258]) - - Args: - mu (:obj:`tensor` or :obj:`float`): The mean of the Normal distribution - sigma (:obj:`tensor` or :obj:`float`): The standard deviation of the Normal distribution - - Returns: - PyTorch Normal distribution: The Normal distribution class - - """ - return torch.distributions.Normal(mu, sigma) - - def to_numpy(self, tensor_in): - """ - Convert the PyTorch tensor to a :class:`numpy.ndarray`. - - Example: - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> tensor = pyhf.tensorlib.astensor([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) - >>> tensor - tensor([[1., 2., 3.], - [4., 5., 6.]]) - >>> numpy_ndarray = pyhf.tensorlib.to_numpy(tensor) - >>> numpy_ndarray - array([[1., 2., 3.], - [4., 5., 6.]]) - >>> type(numpy_ndarray) - - - Args: - tensor_in (:obj:`tensor`): The input tensor object. - - Returns: - :class:`numpy.ndarray`: The tensor converted to a NumPy ``ndarray``. - - """ - return tensor_in.numpy() - - def transpose(self, tensor_in): - """ - Transpose the tensor. - - Example: - >>> import pyhf - >>> pyhf.set_backend("pytorch") - >>> tensor = pyhf.tensorlib.astensor([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) - >>> tensor - tensor([[1., 2., 3.], - [4., 5., 6.]]) - >>> pyhf.tensorlib.transpose(tensor) - tensor([[1., 4.], - [2., 5.], - [3., 6.]]) - - Args: - tensor_in (:obj:`tensor`): The input tensor object. - - Returns: - PyTorch FloatTensor: The transpose of the input tensor. - - .. versionadded:: 0.7.0 - """ - return tensor_in.transpose(0, 1) diff --git a/tests/conftest.py b/tests/conftest.py index a0e19207eb..7b7919f209 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -16,7 +16,7 @@ def pytest_addoption(parser): action="append", type=str, default=[], - choices=["pytorch", "jax", "minuit"], + choices=["jax", "minuit"], help="list of backends to disable in tests", ) parser.addoption( @@ -78,15 +78,13 @@ def reset_backend(): scope='function', params=[ (("numpy_backend", dict()), ("scipy_optimizer", dict())), - (("pytorch_backend", dict()), ("scipy_optimizer", dict())), - (("pytorch_backend", dict(precision="64b")), ("scipy_optimizer", dict())), (("jax_backend", dict()), ("scipy_optimizer", dict())), ( ("numpy_backend", dict(poisson_from_normal=True)), ("minuit_optimizer", dict()), ), ], - ids=['numpy', 'pytorch', 'pytorch64', 'jax', 'numpy_minuit'], + ids=['numpy', 'jax', 'numpy_minuit'], ) def backend(request): # a better way to get the id? all the backends we have so far for testing diff --git a/tests/constraints.txt b/tests/constraints.txt index 4596667894..70297547e7 100644 --- a/tests/constraints.txt +++ b/tests/constraints.txt @@ -11,8 +11,6 @@ numpy==1.21.0 # constrained by jax v0.4.1 uproot==4.1.1 # minuit iminuit==2.7.0 # c.f. PR #1895 -# torch -torch==1.10.0 # jax # Use Google Cloud Storage buckets for long term wheel support # c.f. https://github.com/google/jax/discussions/7608#discussioncomment-1269342 diff --git a/tests/test_backend_consistency.py b/tests/test_backend_consistency.py index 4d65799904..b2bfe5c2ae 100644 --- a/tests/test_backend_consistency.py +++ b/tests/test_backend_consistency.py @@ -105,7 +105,6 @@ def test_hypotest_qmu_tilde( backends = [ pyhf.tensor.numpy_backend(precision='64b'), - pyhf.tensor.pytorch_backend(precision='64b'), pyhf.tensor.jax_backend(precision='64b'), ] @@ -128,10 +127,6 @@ def test_hypotest_qmu_tilde( numpy_ratio = np.divide(test_statistic, test_statistic[0]) numpy_ratio_delta_unity = np.absolute(np.subtract(numpy_ratio, 1)) - # compare tensor libraries to each other - tensors_ratio = np.divide(test_statistic[1], test_statistic[2]) - tensors_ratio_delta_unity = np.absolute(np.subtract(tensors_ratio, 1)) - try: assert (numpy_ratio_delta_unity < tolerance['numpy']).all() except AssertionError: @@ -139,10 +134,3 @@ def test_hypotest_qmu_tilde( f"Ratio to NumPy+SciPy exceeded tolerance of {tolerance['numpy']}: {numpy_ratio_delta_unity.tolist()}" ) assert False - try: - assert (tensors_ratio_delta_unity < tolerance['tensors']).all() - except AssertionError: - print( - f"Ratio between tensor backends exceeded tolerance of {tolerance['tensors']}: {tensors_ratio_delta_unity.tolist()}" - ) - assert False diff --git a/tests/test_init.py b/tests/test_init.py index c247d40ac9..495b3aeab1 100644 --- a/tests/test_init.py +++ b/tests/test_init.py @@ -7,12 +7,6 @@ "param", [ ["numpy", "numpy_backend", "numpy_backend", pytest.raises(ImportError)], - [ - "torch", - "pytorch_backend", - "pytorch_backend", - pytest.raises(pyhf.exceptions.ImportBackendError), - ], [ "jax", "jax_backend", @@ -20,7 +14,7 @@ pytest.raises(pyhf.exceptions.ImportBackendError), ], ], - ids=["numpy", "pytorch", "jax"], + ids=["numpy", "jax"], ) def test_missing_backends(isolate_modules, param): backend_name, module_name, import_name, expectation = param diff --git a/tests/test_interpolate.py b/tests/test_interpolate.py index cb5d605368..d544aee37f 100644 --- a/tests/test_interpolate.py +++ b/tests/test_interpolate.py @@ -118,13 +118,7 @@ def test_interpolator(backend, interpcode, random_histosets_alphasets_pair): def test_validate_implementation(backend, interpcode, random_histosets_alphasets_pair): histogramssets, alphasets = random_histosets_alphasets_pair - # single-float precision backends, calculate using single-floats - if pyhf.tensorlib.name in ['pytorch']: - abs_tolerance = 1e-6 - histogramssets = np.asarray(histogramssets, dtype=np.float32) - alphasets = np.asarray(alphasets, dtype=np.float32) - else: - abs_tolerance = 1e-12 + abs_tolerance = 1e-12 histogramssets = histogramssets.tolist() alphasets = pyhf.tensorlib.astensor(alphasets.tolist()) diff --git a/tests/test_optim.py b/tests/test_optim.py index fdb6707718..1caf26862d 100644 --- a/tests/test_optim.py +++ b/tests/test_optim.py @@ -12,8 +12,6 @@ # from https://docs.scipy.org/doc/scipy/tutorial/optimize.html#nelder-mead-simplex-algorithm-method-nelder-mead -@pytest.mark.skip_pytorch -@pytest.mark.skip_pytorch64 @pytest.mark.skip_numpy_minuit def test_scipy_minimize(backend, capsys): tensorlib, _ = backend @@ -34,10 +32,9 @@ def rosen(x): 'tensorlib', [ pyhf.tensor.numpy_backend, - pyhf.tensor.pytorch_backend, pyhf.tensor.jax_backend, ], - ids=['numpy', 'pytorch', 'jax'], + ids=['numpy', 'jax'], ) @pytest.mark.parametrize( 'optimizer', @@ -61,17 +58,13 @@ def test_minimize(tensorlib, optimizer, do_grad, do_stitch): 'do_grad-minuit-numpy': None, # no grad, scipy, 64b 'no_grad-scipy-numpy': [0.49998815367220306, 0.9999696999038924], - 'no_grad-scipy-pytorch': [0.49998815367220306, 0.9999696999038924], 'no_grad-scipy-jax': [0.4999880886490433, 0.9999696971774877], # do grad, scipy, 64b - 'do_grad-scipy-pytorch': [0.49998837853531425, 0.9999696648069287], 'do_grad-scipy-jax': [0.49998837853531414, 0.9999696648069285], # no grad, minuit, 64b - quite consistent 'no_grad-minuit-numpy': [0.5000493563629738, 1.0000043833598724], - 'no_grad-minuit-pytorch': [0.5000493563758468, 1.0000043833508256], 'no_grad-minuit-jax': [0.5000493563528641, 1.0000043833614634], # do grad, minuit, 64b - 'do_grad-minuit-pytorch': [0.500049321728735, 1.00000441739846], 'do_grad-minuit-jax': [0.500049321731032, 1.0000044174002167], }[identifier] @@ -101,7 +94,7 @@ def test_optimizer_mixin_extra_kwargs(optimizer): @pytest.mark.parametrize( 'backend,backend_new', - itertools.permutations([('numpy', False), ('pytorch', True), ('jax', True)], 2), + itertools.permutations([('numpy', False), ('jax', True)], 2), ids=lambda pair: f'{pair[0]}', ) def test_minimize_do_grad_autoconfig(mocker, backend, backend_new): diff --git a/tests/test_public_api.py b/tests/test_public_api.py index 3fad05ba57..d788a3379a 100644 --- a/tests/test_public_api.py +++ b/tests/test_public_api.py @@ -17,7 +17,7 @@ def model_setup(backend): return model, data, init_pars -@pytest.mark.parametrize("backend_name", ["numpy", "pytorch", "PyTorch"]) +@pytest.mark.parametrize("backend_name", ["numpy", "jax"]) def test_set_backend_by_string(backend_name): pyhf.set_backend(backend_name) assert isinstance( @@ -43,7 +43,7 @@ def test_set_precision_by_string(precision_level): assert pyhf.tensorlib.precision == precision_level.lower() -@pytest.mark.parametrize("backend_name", [b"numpy", b"pytorch"]) +@pytest.mark.parametrize("backend_name", [b"numpy", b"jax"]) def test_set_backend_by_bytestring(backend_name): pyhf.set_backend(backend_name) assert isinstance( @@ -98,7 +98,7 @@ def test_supported_precision(precision_level): def test_custom_backend_name_supported(): class custom_backend: def __init__(self, **kwargs): - self.name = "pytorch" + self.name = "jax" self.precision = '64b' def _setup(self): @@ -143,14 +143,14 @@ def __init__(self, **kwargs): assert pyhf.optimizer.name == optimizer.name -@pytest.mark.parametrize("backend_name", ["numpy", "pytorch", "PyTorch"]) +@pytest.mark.parametrize("backend_name", ["numpy", "jax"]) def test_backend_no_custom_attributes(backend_name): pyhf.set_backend(backend_name) with pytest.raises(AttributeError): pyhf.tensorlib.nonslotted = True -@pytest.mark.parametrize("backend_name", ["numpy", "pytorch", "PyTorch"]) +@pytest.mark.parametrize("backend_name", ["numpy", "jax"]) def test_backend_slotted_attributes(backend_name): pyhf.set_backend(backend_name) for attr in ["name", "precision", "dtypemap", "default_do_grad"]: diff --git a/tests/test_scripts.py b/tests/test_scripts.py index 1d9f6075d1..7eef61c5cd 100644 --- a/tests/test_scripts.py +++ b/tests/test_scripts.py @@ -192,7 +192,7 @@ def test_import_usingMounts_badDelimitedPaths(datadir, tmp_path, script_runner): assert 'is not a valid colon-separated option' in ret.stderr -@pytest.mark.parametrize("backend", ["numpy", "pytorch", "jax"]) +@pytest.mark.parametrize("backend", ["numpy", "jax"]) def test_fit_backend_option(tmp_path, script_runner, backend): temp = tmp_path.joinpath("parsed_output.json") command = f"pyhf xml2json validation/xmlimport_input/config/example.xml --basedir validation/xmlimport_input/ --output-file {temp}" @@ -207,7 +207,7 @@ def test_fit_backend_option(tmp_path, script_runner, backend): assert "mle_parameters" in ret_json -@pytest.mark.parametrize("backend", ["numpy", "pytorch", "jax"]) +@pytest.mark.parametrize("backend", ["numpy", "jax"]) def test_cls_backend_option(tmp_path, script_runner, backend): temp = tmp_path.joinpath("parsed_output.json") command = f'pyhf xml2json validation/xmlimport_input/config/example.xml --basedir validation/xmlimport_input/ --output-file {temp}' diff --git a/tests/test_simplemodels.py b/tests/test_simplemodels.py index 29a2be33e6..e7201bf460 100644 --- a/tests/test_simplemodels.py +++ b/tests/test_simplemodels.py @@ -38,8 +38,6 @@ def test_uncorrelated_background(backend): # See https://github.com/scikit-hep/pyhf/issues/1654 -@pytest.mark.fail_pytorch -@pytest.mark.fail_pytorch64 @pytest.mark.fail_jax def test_correlated_background_default_backend(default_backend): model = pyhf.simplemodels.correlated_background( @@ -56,8 +54,6 @@ def test_correlated_background_default_backend(default_backend): # See https://github.com/scikit-hep/pyhf/issues/1654 -@pytest.mark.fail_pytorch -@pytest.mark.fail_pytorch64 @pytest.mark.fail_jax def test_uncorrelated_background_default_backend(default_backend): model = pyhf.simplemodels.uncorrelated_background( diff --git a/tests/test_tensor.py b/tests/test_tensor.py index 2b0698e22a..448cedadc2 100644 --- a/tests/test_tensor.py +++ b/tests/test_tensor.py @@ -251,8 +251,6 @@ def test_shape(backend): ) -@pytest.mark.fail_pytorch -@pytest.mark.fail_pytorch64 def test_pdf_calculations(backend): tb = pyhf.tensorlib # FIXME @@ -295,50 +293,6 @@ def test_pdf_calculations(backend): ) == pytest.approx([0.4151074974205947, 0.3515379040027489, 0.2767383316137298]) -# validate_args in torch.distributions raises ValueError not nan -@pytest.mark.only_pytorch -@pytest.mark.only_pytorch64 -def test_pdf_calculations_pytorch(backend): - tb = pyhf.tensorlib - - values = tb.astensor([0, 0, 1, 1]) - mus = tb.astensor([0, 1, 0, 1]) - sigmas = tb.astensor([0, 0, 0, 0]) - for x, mu, sigma in zip(values, mus, sigmas): - with pytest.raises(ValueError): - _ = tb.normal_logpdf(x, mu, sigma) - assert tb.tolist( - tb.normal_logpdf( - tb.astensor([0, 0, 1, 1]), - tb.astensor([0, 1, 0, 1]), - tb.astensor([1, 1, 1, 1]), - ) - ) == pytest.approx( - [ - -0.91893853, - -1.41893853, - -1.41893853, - -0.91893853, - ], - ) - - # Allow poisson(lambda=0) under limit Poisson(n = 0 | lambda -> 0) = 1 - assert tb.tolist( - tb.poisson(tb.astensor([0, 0, 1, 1]), tb.astensor([0, 1, 0, 1])) - ) == pytest.approx([1.0, 0.3678794503211975, 0.0, 0.3678794503211975]) - with pytest.warns(RuntimeWarning, match="divide by zero encountered in log"): - assert tb.tolist( - tb.poisson_logpdf(tb.astensor([0, 0, 1, 1]), tb.astensor([0, 1, 0, 1])) - ) == pytest.approx( - np.log([1.0, 0.3678794503211975, 0.0, 0.3678794503211975]).tolist() - ) - - # Ensure continuous approximation is valid - assert tb.tolist( - tb.poisson(n=tb.astensor([0.5, 1.1, 1.5]), lam=tb.astensor(1.0)) - ) == pytest.approx([0.4151074974205947, 0.3515379040027489, 0.2767383316137298]) - - def test_boolean_mask(backend): tb = pyhf.tensorlib assert tb.tolist( @@ -365,11 +319,6 @@ def test_percentile(backend): assert tb.tolist(tb.percentile(a, 50, axis=1)) == [7.0, 2.0] -# FIXME: PyTorch doesn't yet support interpolation schemes other than "linear" -# c.f. https://github.com/pytorch/pytorch/pull/59397 -# c.f. https://github.com/scikit-hep/pyhf/issues/1693 -@pytest.mark.fail_pytorch -@pytest.mark.fail_pytorch64 def test_percentile_interpolation(backend): tb = pyhf.tensorlib a = tb.astensor([[10, 7, 4], [3, 2, 1]]) @@ -532,7 +481,7 @@ def test_tensor_precision(backend): @pytest.mark.parametrize( 'tensorlib', - ['numpy_backend', 'jax_backend', 'pytorch_backend'], + ['numpy_backend', 'jax_backend'], ) @pytest.mark.parametrize('precision', ['64b', '32b']) def test_set_tensor_precision(tensorlib, precision): @@ -574,7 +523,7 @@ def test_trigger_tensorlib_changed_precision(mocker): @pytest.mark.parametrize( 'tensorlib', - ['numpy_backend', 'jax_backend', 'pytorch_backend'], + ['numpy_backend', 'jax_backend'], ) @pytest.mark.parametrize('precision', ['64b', '32b']) def test_tensorlib_setup(tensorlib, precision, mocker): diff --git a/tests/test_validation.py b/tests/test_validation.py index 6984f4e001..0c5b1bcc18 100644 --- a/tests/test_validation.py +++ b/tests/test_validation.py @@ -996,7 +996,6 @@ def test_shapesys_nuisparfilter_validation(): [ pyhf.tensor.numpy_backend, pyhf.tensor.jax_backend, - pyhf.tensor.pytorch_backend, ], ) @pytest.mark.parametrize('optimizer', ['scipy', 'minuit']) @@ -1020,7 +1019,6 @@ def test_optimizer_stitching(backend, optimizer): 'backend', [ pyhf.tensor.jax_backend, - pyhf.tensor.pytorch_backend, ], ) @pytest.mark.parametrize('optimizer,rtol', [('scipy', 1e-6), ('minuit', 1e-3)])