diff --git a/.github/workflows/build-wheels.yml b/.github/workflows/build-wheels.yml index 3fdb276..47065ab 100644 --- a/.github/workflows/build-wheels.yml +++ b/.github/workflows/build-wheels.yml @@ -85,23 +85,23 @@ jobs: name: cibw-sdist path: dist/*.tar.gz - upload_pypi: - name: Publish to PyPI - needs: [build_wheels, build_sdist] - runs-on: ubuntu-latest - environment: - name: pypi - permissions: - id-token: write - # if: github.event_name == 'release' && github.event.action == 'published' - steps: - - uses: actions/download-artifact@v4 - with: - pattern: cibw-* - path: dist - merge-multiple: true +# upload_pypi: +# name: Publish to PyPI +# needs: [build_wheels, build_sdist] +# runs-on: ubuntu-latest +# environment: +# name: pypi +# permissions: +# id-token: write +# # if: github.event_name == 'release' && github.event.action == 'published' +# steps: +# - uses: actions/download-artifact@v4 +# with: +# pattern: cibw-* +# path: dist +# merge-multiple: true - - uses: pypa/gh-action-pypi-publish@release/v1 +# - uses: pypa/gh-action-pypi-publish@release/v1 # upload_test_pypi: # name: Publish to TestPyPI diff --git a/docs/examples/explicit-model.ipynb b/docs/examples/explicit-model.ipynb index 7b5d4fd..fbe3d19 100644 --- a/docs/examples/explicit-model.ipynb +++ b/docs/examples/explicit-model.ipynb @@ -29,7 +29,7 @@ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", - "from odrpack import odr" + "from odrpack import odr_fit" ] }, { @@ -45,19 +45,19 @@ "metadata": {}, "outputs": [], "source": [ - "x = np.array([0.100, 0.300, 0.400, 0.500, 0.600, 0.700, 0.800])\n", - "y = np.array([0.059, 0.243, 0.364, 0.486, 0.583, 0.721, 0.824])" + "xdata = np.array([0.100, 0.300, 0.400, 0.500, 0.600, 0.700, 0.800])\n", + "ydata = np.array([0.059, 0.243, 0.364, 0.486, 0.583, 0.721, 0.824])" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray:\n", " b1, b2 = beta\n", - " return (b1*x**2 + x*(1-x))/(b1*x**2 + 2*x*(1-x) + b2*(1-x)**2)" + " return (b1*x**2 + x*(1 - x))/(b1*x**2 + 2*x*(1 - x) + b2*(1 - x)**2)" ] }, { @@ -95,15 +95,15 @@ "sigma_x = 0.01\n", "sigma_y = 0.05\n", "\n", - "wd = 1/sigma_x**2\n", - "we = 1/sigma_y**2" + "weight_x = 1/sigma_x**2\n", + "weight_y = 1/sigma_y**2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We can now launch the regression! If you want to see a brief computation report, set `iprint=1001`." + "We can now launch the regression! If you want to see a brief computation report, set `report='short'`." ] }, { @@ -112,7 +112,9 @@ "metadata": {}, "outputs": [], "source": [ - "sol = odr(f, beta0, y, x, lower=lower, upper=upper, we=we, wd=wd)" + "sol = odr_fit(f, xdata, ydata, beta0,\n", + " bounds=(lower, upper),\n", + " weight_x=weight_x, weight_y=weight_y)" ] }, { @@ -199,7 +201,7 @@ "_, ax = plt.subplots()\n", "\n", "# Plot experimental data\n", - "ax.plot(x, y, 'o')\n", + "ax.plot(xdata, ydata, 'o')\n", "\n", "# Plot fitted model\n", "xm = np.linspace(0.0, 1.0, 100)\n", diff --git a/docs/examples/implicit-model.ipynb b/docs/examples/implicit-model.ipynb index 0d9e883..281234e 100644 --- a/docs/examples/implicit-model.ipynb +++ b/docs/examples/implicit-model.ipynb @@ -37,7 +37,7 @@ "import numpy as np\n", "from matplotlib.patches import Ellipse\n", "\n", - "from odrpack import odr" + "from odrpack import odr_fit" ] }, { @@ -76,10 +76,10 @@ " [-3.44, -4.86]]\n", "\n", "# We need shape (2, n)\n", - "X = np.array(X).T\n", + "Xdata = np.array(X).T\n", "\n", - "# Y is not used, but is (currently) required by odr\n", - "Y = np.full(X.shape[-1], np.nan)" + "# Ydata is not used, but is required\n", + "Ydata = np.zeros(Xdata.shape[-1])" ] }, { @@ -128,23 +128,26 @@ "metadata": {}, "outputs": [], "source": [ - "wd = 1.0" + "weight_x = 1.0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We can now launch the regression! As the problem is implicit, we set `job=1`. If you want to see a brief computation report, set `iprint=1001`." + "We can now launch the regression! As the problem is implicit, we set `task='implicit-ODR'`. If you want to see a brief computation report, set `report=short`." ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "sol = odr(f, beta0, Y, X, lower=lower, upper=upper, wd=wd, job=1)" + "sol = odr_fit(f, Xdata, Ydata, beta0,\n", + " bounds=(lower, upper),\n", + " weight_x=weight_x,\n", + " task='implicit-ODR')" ] }, { @@ -182,7 +185,7 @@ { "data": { "text/plain": [ - "array([-0.99938056, -2.93104944, 3.86422426, 3.15663765, -0.90359046])" + "array([-0.99938103, -2.93104802, 3.86422581, 3.15663797, -0.90359144])" ] }, "execution_count": 8, @@ -218,7 +221,7 @@ }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGwCAYAAABRgJRuAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVm5JREFUeJzt3Xl8VPW9//HXzGQPyUAgIQmEJSBoAEHAIIgKipouqL2ttlqtWmuVq+1Ve1Ww7aXxttJWq96qF5dbl5/YulYp1sYFd0GiQlBAEJA1CwkBJiGQbeb8/vgmmYQkEGAyZ5b38/E4zfecOQkfUsy88z3fxWFZloWIiIhImHPaXYCIiIhIICjUiIiISERQqBEREZGIoFAjIiIiEUGhRkRERCKCQo2IiIhEBIUaERERiQgxdhcQTD6fj7KyMlJSUnA4HHaXIyIiIj1gWRa1tbVkZ2fjdHbfHxNVoaasrIycnBy7yxAREZFjsGPHDgYPHtzt61EValJSUgDzTUlNTbW5GhEREemJmpoacnJy2t7HuxNVoab1kVNqaqpCjYiISJg50tARDRQWERGRiKBQIyIiIhFBoUZEREQigkKNiIiIRASFGhEREYkICjUiIiISERRqREREJCIo1IiIiEhEiKrF90SikmVBcy00esDXAN4G8DWatq+x5bx9uxGsZnDEgMPlP5ztz9u1XfEQkwKxKf6Prni7/9YiEoUUakTCieWD+io4WAYHS83Hhmpo3Nty7DMfm/b5z5v2mc8LJmcsxPTpHHbi+kFCBsRnQOJA8zEhAxIGmo8xfUCbzYrIMVKoEQklDdVQ8xXs/9qElgOl/vByoBTqy8HXZHeVR+Zr8geto+FKbBd6siB5GPQZZj4mDzUf4/op+IhIlxRqRIKtqQZqN0LNRqj9yrRbPx5tCOiOwwVxfSG2n/kY1w9iU01ocMaZx0POOHDGtztv13a4wPKaw9fsb1te82iq/bm3HppqzSOu9h9b2837obmuZ3V7D0LdNnN0JyalXdBpOfoMh9TR0GeEHn2JRDGFGpHe4muCmvWwt6TlWAWedVC/69i+Xnw6JGZD0iBIHNTyMdv0bMT1g9i+5mNc39B7jOPzgrcOGvaYv399JTRU+tttH1vaDbsBq+uv1VwL+74wx6EcTkjOhdQTwX1iy8cx4M4zoU5EIppCjUggNO2Hfas7Bph9a8wA3J5KGgIpJ7QcIyEpxx9gErPCuwfC6QJnqgkWfYYd+X5fswk3dVth/1bzsfXYvxUObDcDmg9l+WD/JnOUvdrxtaQcE3D6joW+EyBtEqSOMkFIRCKCQo3I0bIsqNsClR9CVctRs55uexbaSxhoHpO0hZdR5mOfERCT2Oulhw1njAl0SYMg/fTOr1s+OFjhDzq1m6Bmg/n/oWY9eA90/pwDO8xRXuS/FtMH+p1iAk7rkTLKhDARCTsKNSJH4ms2vTDtQ0x9xRE+yWHCS78J5k2z3wToO97M+JHj53BCUrY50qd1fM3ymUHVNevN4z7PWv/R5Ol4b/N+qPrAHK1ikv1Bp99EGDDV9JyF0uM8EemSw7KsHvx6GRlqampwu914PB5SU/V8Xbph+cwjpPIi2PUO7F5++IGujpiWN8GJLSFmAvQdZ94cJXRYlplFtu8L2POZ/ziw/cifmzAQ0qdD+hmQMd0EVKd+JxQJlp6+fyvUiIBZ+6XiTSgrgorXzYDV7sSmwoBp5rFI+nTonw8xScGrVQKrvgr2rIS97YLO4WZfgXlsNWCq+f8/4wzoP0X/BkR6kUJNFxRqpI2vGaqLTW9MWRHs+ZRux8QkDjJvXOnTzeEeqzEXka6h2gSd6mKo+gh2f2Sm4nfHEWMeVw2cCVnnmdAbzgO7RUKMQk0XFGqinLcRKt6CHS/AzsXdrwkTkwwDz4HsAsg8D/rkajxFtPN5wbMGKj9oGVf1gXmU1Z2YZMiYYQJO5nlmfJX+DYkcM4WaLijURCFvg3mstL0lyBw6ULRV35Mhq8AEmQGngysuuHVKeLEsM+uq6kN/0Kn5svv7k3L8ASfzHIjvH7RSRSKBQk0XFGqihLceyt8wQab0H10/NojpA9nfgKxvQNb5ZhaNyPGorzQ9geVvQMUbcLC8mxsdkDYZBn0LBl9kArV6cUQOS6GmCwo1EcyyzCylzX+BHS92E2RSYPAFMORiE2RcCcGvU6KDZZkp5K0Bp/I9E7a7kjzMhJvBF5oxW5pVJdKJQk0XFGoi0MEK2PI0fP14ywJ4h4hNhUEXtgSZcxVkxB7eevOIqvwNc+xb3fV9cWkwaLYJOVnnaUaVSAuFmi4o1EQIXzOUvWaCTOmrZlPF9mJSIOc7JshknqtZKBJ6DpSaR6M7XoFdb5tNQg/lSjBjcAZfBDkXmX29RKKUQk0XFGrCXO1m2PQobPl/Xa/om3EW5P4YhnxXC99J+GjcB2X/gp2vmLDevL/zPc5YM/5r2GWmJ0c9OBJlFGq6oFAThizLTJ9df5+ZvXToWjKJ2ZB7FeRebZayFwln3gbTc7NzsTm6Cu8xyab3Zuhl5pGqMzboZYoEm0JNFxRqwoi3EbY/b8LM3pUdX3PGwqALTK9M1nkaWCmRyfLB7hVm4Pu2v3U9myq+P+RcbHpw0k/XjuMSsRRquqBQEwYa9sCmR+CrBzsvbpaYDaNuhBE/gYR0e+oTsYPPC1Xvw9a/wvYXoWlf53uScky4yf0xpI4KeokivUmhpgsKNSFs/1b48o/w9ZPgPdjxtbRJMPpmM/BXi+JJtPM2QPnrJuCU/qPzfy8AGWea8J/zXY2/kYigUNMFhZoQVLcN1vwOvn7ikBkgDrNux4m3mLU7tDiZSGdN+83Ym21/NUHn0JmAsW4Y9kMTcNJOsadGkQBQqOmCQk0IqdsOa+8y07J9Tf7rMX1M9/non0PKCPvqEwk39VVmzabN/9f1lg39JsLIn8DQSyGub9DLEzkeCjVdUKgJAQd2wtoFsPmxjmEmNhVG3wQn3qT1OESOR9vq2v8H254D74GOr7sSzaPckdfDgNPUCyphQaGmCwo1NjpYbsLMpkfA1+i/HpMCo/8DTrwZ4tPsq08kEjXVwLZnYdP/wZ5POr+edqrpFR1yicarSUhTqOmCQo0Nmuvgy3tg3R87/sYYkwyjfg4n/UI7FosEw97VZm+0LU93nj2VkAknXA8jr4PETFvKEzkchZouKNQEkeUzK/+u/mXHqdmuJBj9MzjxPyFhgH31iUSr5oOw/TnY8D+wt6Tja85YGPID03vTf7It5Yl0RaGmCwo1QVL1EXz6846L5jli4IQ5MPZXkJBhX20iYliW2WRzw59h59/NLyLtDZhmwk3Ov2nVYrFdT9+/tRSrBM6BUii5HbY+0/H64Athwh+1IJhIKHE4IOMMc9Rth43/C5seg8Y95vXdy8yRlAMn/aeZFq41byTEqadGjp+vyWxnsOZOM4amVd/xMOk+GDjTvtpEpOeaD5hfSjb8GTxrOr4WP8DMUBx1g6aES9D19P077DYKeeihhxg2bBgJCQlMmTKF4uJiu0uKbntWwev5poemNdDEpcGpC6HgMwUakXASkwQjr4Vvfg7nvA3Z3/a/1rAbPv8VvDIESubCwS422xSxWViFmueee45bbrmF+fPns3LlSsaPH8/5559PZWWl3aVFH289lMyD10/1DzZ0OOGEG2D2RjOTwumytUQROUYOh/mFZMYS+MZqs2Bf62aZzbWw7g+weBh88u+wf4utpYq0F1aPn6ZMmcKpp57Kgw8+CIDP5yMnJ4ef/exnzJ07t9P9DQ0NNDQ0tJ3X1NSQk5Ojx0/Hq/JDWHEN1H7lv9Z3HEz5C/Q/1b66RKT31G6GL+82W5q0X2vK4TKhJ28u9B1jX30S0SLu8VNjYyOfffYZs2bNarvmdDqZNWsWy5cv7/JzFixYgNvtbjtycnKCVW5kaqqFT26Et87wBxpnLIy7E87/VIFGJJKljID8h+GCLWbgcEyyuW55YesieG0cfPRDqNlob50S1cIm1OzevRuv18vAgQM7XB84cCAVFV0/2503bx4ej6ft2LFjRzBKjUwVS+GfY2DjQ/5r/U+DglUw7tdajVQkWiRlwyl3w4XbYVyhGUMHgIV367Msf+5iFr/4G5avXYfXFzYPAiRCRPSU7vj4eOLj4+0uI7z5muDz/zLP0Gn5AeVKgvF3wagbNW5GJFrFp8G4/zKrgm9cSNGHb1K4/VLKm9LN659uISvxS+bPzqNg4kn21ipRI2xCzYABA3C5XOzatavD9V27dpGZqWW9e8X+rfDRpVD9sf/awHNgymPQZ7htZYlICIlJpsj7Q+ZsPgmLjj0zFQfjmPP8ZhaWLqHgvGu0JYr0urB5/BQXF8ekSZNYunRp2zWfz8fSpUuZOnWqjZVFqO0vwL8m+AONI8Z0OZ/9hgKNiLTx+iwKl6xriTMdd/y2cAIWhcUZeF/Jhc9/YzbZFOklYRNqAG655RYee+wxnnrqKb788kvmzJlDXV0dV199td2lRY7mA7Dip/DhJdDkMdf65MK5H5nBgY6w+icjIr2seMseyj313b5u4aS8KZ1iz1BYUwj/GAEbF4KvOYhVSrQIm8dPAN///vepqqriv/7rv6ioqGDChAkUFRV1Gjwsx2jfGvjo++BZ57829Adw6sMQ57avLhEJWZW13QeaDvc1t2xg27DbrG/z1YNwyp8gu6AXq5NoE1br1BwvbZNwGDtehuVX+FcFdiXB5Acg92qzEJeISBeWb67m0sc+PuJ9f7s8h6meBbDt2Y4vZBXAKfdojRs5rIhbp0Z6iWXBmt/CB//mDzR9T4aCT2HEjxVoROSw8oenkeVOoLufFA4gy51Aft44OP1vcO4y6D/Ff0N5EfxrvOm9qa8KRskSwRRqolnzATO76fNf+68NvRTO+xjcmoIpIkfmcjqYPzsPOHSYsP98/uw8XM6Ws/SpcN5ymPZXSBpirlleM85myUhYdzd4GxA5Fgo10erATnjrTNj+XMsFh1l7ZtozEJNoa2kiEl4Kxmax8PKJZLoTOlzPdCew8PKJFIzN6vgJDgcMuxS+vR7G/w5i+pjrTTVQchu8ehLs/EeQqpdIojE10Wj3Cnj/IqhvWYk5po8JM4MvsLUsEQlvXp9F8ZY9VNbWk5GSQP7wNH8PzeEcrDA9xpv/Au3Xuhl0AUz+MyQP7bWaJTz09P1boSbabHsOll8Jvpbu3eRhcNY/zIaUIiJ22rsaVt4Mu97xX3MlmZWLT7zF7DUnUUkDhaWzjQ+bMTStgSbjLDj/EwUaEQkN/cbD2Uth2t8goWWleO8BKJlrFgOtfN/W8iT0KdREA8uCtQvgkzm0de3m/hhmvgEJA2wtTUSkA4cDhv3AjLcZdaN/wU/POnjrLPj4as2Skm4p1EQ6y4KS22H1Hf5rJ90GU/5PO2uLSOiKc5u1ss4vhrTJ/utfPwmvjoZNj4Hls608CU0KNZHM54Xia+HLu/3XJvweTvmD1p8RkfCQNsksMzH5IYhtWdm8cS8U/9TM4Kz5yt76JKQo1EQqbwN89IOW2QQADsh/BPJut7UsEZGj5nTBqH83j6SG/dB/veojs3Df+vvML3ES9RRqIpG3Ad7/Dux40Zw7Y81KniN/am9dIiLHIzETpi0yg4n7jDDXvPWw8hZYepZ6bUShJuL4mkwPTfm/zLkrEc78Bwz9vr11iYgESubZ8M3VMOrn/mvqtREUaiKL5TMzA3a+Ys5dSTDzde2CKyKRJyYZJv8PzHpPvTbSRqEmUliWmbK99Rlz7oyDsxZDxhn21iUi0psyzlSvjbRRqIkElgWr/hM2PWrOHS6Y/gJkzrK3LhGRYDhsr80MqNtma3kSPAo1keCLQlh/b8uJA6Yu0j5OIhJ9uuy1+RBeGw/bnrevLgkahZpwt/4+WFPoP5/ymFmNU0QkGrXvtUkeZq41eeCj78PH10BzXadP8foslm+uZnFJKcs3V+P1Rc2WiBEnxu4C5Djs/Aes/IX/fOL9MOIa28oREQkZGWfCN0rgk+th27Pm2tePm56b0/8GaRMBKFpTTuGSdZR76ts+NcudwPzZeRSMzbKhcDke6qkJV3tXw7LLaNvLaex8OPE/bC1JRCSkxLlh2l/htCdNDw5A7Vfwxmnw5b0UfVHGnEUrOwQagApPPXMWraRoTXnwa5bjolATjg5WwHuz/d2oQ38A4+bbW5OISChyOCD3SihYabZcAPA14V15K4Uvvk9XD5parxUuWadHUWFGoSbceOvh/YvgwA5z3j8fpjyuvZxERA4ndRScuwxO+k8AiuvGUN6Q0u3tFlDuqad4y54gFSiBoFATTiwLPv4xVK8w50mD4cxXICbR1rJERMKCKw5OuRtmvk6lY1iPPqWytv7IN0nI0EDhcLL2d7Dtb6btSoKzlkCiBrKJiByVrPPImJ4Lmzcc8daMlIQgFCSBop6acFFWBJ//2n8+7RnoN8G2ckREwln+6BFkuRNwdDmqBhyYWVD5w9OCW5gcF4WacHCwHJb/yH8+/i7Iuci2ckREwp3L6WD+7DzAwaEjEh34AIv5s/NwOTVeMZwo1IQ6ywfLroCGKnOe/S3Im2tvTSIiEaBgbBYLL59IprvjI6bM2GoWDr2Lgvrfg6/JpurkWGhMTahb9wfYtdS0E7PhtCc000lEJEAKxmZxbl4mxVv2ULmvhoyyheTX/AmXwwdfLYe9K2H68xq/GCYclmVFzST8mpoa3G43Ho+H1NRUu8s5sqpl8NaZYHkBB5zzNgycYXdVIiKRy7Jg82Pw6c/A12iuJWSamaYDpthaWjTr6fu3Hj+Fqsa98NGlLYEGGPtrBRoRkd7mcMDIn8KsD8yyGQD1FWa37+0v2FqaHJlCTSiyLFhxLRzYbs7TzzChRkREgmNAPhR8ZvaQArPw6YeXwNoF5me0hCSFmlC0/QXY8ZJpx/Uz07edGv4kIhJUCRkw800YfqX/2uo74OOrwdtoX13SLYWaUNOwBz77mf88/1FIzrGvHhGRaOaKMxM0xt/lv7blKXjnPGiotq8u6ZJCTahZdSvUV5r24Isg57u2liMiEvUcDhgzz8yCcrVM/658D96YCjUb7a1NOlCoCSW73oGvHzft2FSY/KCmb4uIhIohF8M570LCQHNeuxHeOA12vWdnVdKOQk2oaD4IK37qP5/we0gaZF89IiLS2YApcP4KcI8154174J1zYetf7a1LAIWa0LHmv2H/JtNOPx1GXmdvPSIi0rXkoXDeR5BVYM59TbDsh7D+fw77aV6fxfLN1SwuKWX55mq8Ps2iCrSwmVLzu9/9jn/+85+UlJQQFxfHvn377C4pcPZ9AV/ebdrOODM42KG8KSISsmJT4awl8OmNsOkRc23lTVC/C8b/rtPQgaI15RQuWUe5p77tWpY7gfmz8ygYq9WKAyVs3jkbGxu5+OKLmTNnjt2lBN7KX4DVbNpj7gB3nr31iIjIkTlj4NSFMPa//NfWLYAVPwFfc9ulojXlzFm0skOgAajw1DNn0UqK1pQHq+KIFzahprCwkJtvvplx48bZXUpglb8BFW+advJwbVYpIhJOHA44udBM7Gjd7/vrx81Cfd4GvD6LwiXr6OpBU+u1wiXr9CgqQMIm1ByLhoYGampqOhwhxeeFVbf5z8ffBa54++oREZFjM+oGOP1ZcMaa850vw3uzKd5U2qmHpj0LKPfUU7xlT3DqjHARHWoWLFiA2+1uO3JyQmwRu62LYN9q006bDEMvsbceERE5dkMvgbNeBVeSOa94k8qPf9OjT62s7T74SM/ZGmrmzp2Lw+E47LF+/fpj/vrz5s3D4/G0HTt27Ahg9cep+SB8/iv/+Sl/1OBgEZFwl3UenP0mxLoByKgv7tGnZaQk9GZVUcPW2U+/+MUvuOqqqw57T25u7jF//fj4eOLjQ/Rxzld/hgM7TTv7WzBwpr31iIhIYKRPg1nvwjvnk2+tJSu2ioqmAVh0XkzVAWS6E8gfnhb0MiORraEmPT2d9PR0O0uwR0O12ekVTO/MhN/bW4+IiARWvwkw6wNcS89mfvajzNl2Bw6sDsGmtTV/dh4up1aPD4Swed6xfft2SkpK2L59O16vl5KSEkpKSti/f7/dpR299fdDk8e0c6+GvmNtLUdERHpB6iiY9S4FmdtZOPQuMmN3d3g5053Awssnap2aAHJYlhUW88iuuuoqnnrqqU7X33nnHWbMmNGjr1FTU4Pb7cbj8ZCamhrgCnuoaT8sHgKNe8ERAxd8rV24RUQiWc1GWDoD74EKiuvGUBk3gYzTfkP+qOHqoemhnr5/h02oCYSQCDXr74eVN5v28Cth6pP21CEiIsHTEmw4WGbO+46Hc5ZCfH9bywoXPX3/DpvHTxHB1wTr7/Wfn3SrfbWIiEjwpJ4A57wDiS2PmvathrdnmTGWEjAKNcG07Vk40DKtPPvb0HeMvfWIiEjwpI6Cc971B5u9JfD2udCghfcCRaEmWCwL1v3Rf553W/f3iohIZEodZXpsEjLN+d5V8O43obnO3roihEJNsJS9Bp41pj1gKqRPt7ceERGxR+rojsGmegV88D0zREGOi0JNsHQYS3Nbp23pRUQkirhPhJmvt608THkRfHw1WD576wpzCjXBsP9r2PW2afcZCYMvsLceERGxX7+T4awl4GrZImHrM7DyF2a4ghwThZpg+Lrd+jojfqw9nkRExMg4A05/Dhwuc77hflj3B1tLCmd6d+1tPi98/YRpO5ww/Ef21iMiIqFl8AWQ/5j/fPU82PyXLm/1+iyWb65mcUkpyzdX4/WpV6c9W/d+igq73vZP484qgKRB9tYjIiKhZ8TV0FAFJbeb8+KfQlx/yLmo7ZaiNeUULllHuae+7VqWO4H5s/O01UIL9dT0tq8f97dzf2xfHSIiEtpOuhVOvMW0LR989AOo/BAwgWbOopUdAg1AhaeeOYtWUrSmPNjVhiSFmt7UsAd2vGza8QNg0Gx76xERkdDlcMApd8OwK8y5rwE++A7emq8pXLKOrh40tV4rXLJOj6JQqOld2/5m/lECDLscXHH21iMiIqHN4YTT/gKZ55nzht0UL7mpUw9NexZQ7qmneItWJlao6U3bnvO3c6+2rw4REQkfzliY/pxZpA+o3Ofp0adV1nYffKKFQk1vqd8Nuz8y7dTRZj0CERGRnojrC2e9CnH9yIjd26NPyUhJ6N2awoBCTW8p+6d/ZchBF9pbi4iIhJ+UkXDGS+T32UBWbBUOul5t2IGZBZU/PC249YUghZreUvoPf1srCIuIyLEYOBNX/gPMz34UcHQKNq0b7syfnYfLqe13FGp6g7ceyl837fh06H+avfWIiEj4GvlTCvLzWTj0LjJjqzu8lOlOYOHlE7VOTQstvtcbKt72byM/6NvgdNlbj4iIhLdT7qGgZjbnpl5Dcd0YKmPHkXH6H8gfOUg9NO0o1PSG9o+eBunRk4iIHCdnDJz+LK43pjLV8QXwBex2wQlPHfFTo4kePwWaZUHpEtN2JUDWufbWIyIikSHODWe+DDF9zPnWp2Hz/9lbU4hRqAm02q/gYJlpZ8yEmGR76xERkciROhqmtAsyn/4M9qyyr54Qo1ATaFUf+dsZZ9pXh4iIRKah34cTbjBtXwN8eDE09myBvkinUBNoVR/62+mn21eHiIhErol/grRTTXv/ZljxYzP8Icop1ARaa0+NMw76n2pvLSIiEplc8TD9eYjta853/B02/I+tJYUChZpAqq80Y2oA0iaZgcIiIiK9oc8wmPr//OerboWq5baVEwoUagKpapm/nT7dvjpERCQ6DJ4NJ91m2lYzfPSDqB5fo1ATSBpPIyIiwTb+t/5fpA9sh5U32VqOnRRqAml3u56aAdPsq0NERKKHMxamPg0xKeb86ydhxyt2VmQbhZpAsXyw73PT7pMLCen21iMiItGjzzCY1G6gcPFPzTjPKKNQEygHdvj3e3KPsbcWERGJPrlXweALTbuhygSbKJvmrVATKJ4v/e3Uk+yrQ0REopPDAfmPQnzLk4Kdi82jqCiiUBMoNe1CjVuhRkREbJCQAVMe859/9h+wf6tt5QSbQk2gdOipybOvDhERiW6DLzSPogCaa+Hjq8y4zyigUBMoNev8bfeJ9tUhIiIy8X5IGmLale/BV/9raznBolATCJbl76lJHASxqfbWIyIi0S3ODVOf8p+vvgMOlNlXT5Ao1ARCQxU07jFtjacREZFQMHAGjLjGtJtrYeXNtpYTDAo1gVC3zd9OOcG+OkRERNqb8AeIH2Da25+HsiJ76+llYRFqtm7dyjXXXMPw4cNJTExkxIgRzJ8/n8bGRrtLMw5W+NsJmfbVISIi0l58fzjlHv/5J/8OzQfsq6eXhUWoWb9+PT6fj0ceeYS1a9dy33338fDDD3PHHXfYXZpRv8vfTlSoERGREDL8R5Axw7TrtsDa39laTm+KsbuAnigoKKCgoKDtPDc3lw0bNrBw4ULuueeew3xmkNS376kZaF8dIiIih3I44NSF8K+TwdcEX94Nw34I7shbfiQsemq64vF4SEtLO+w9DQ0N1NTUdDh6RfueGj1+EhGRUOM+EfLmmravCYqvj8i1a8Iy1GzatIkHHniA66677rD3LViwALfb3Xbk5OT0TkEH1VMjIiIhLm8e9Blh2lUfwNZn7K2nF9gaaubOnYvD4TjssX79+g6fU1paSkFBARdffDHXXnvtYb/+vHnz8Hg8bceOHTt65y/SoadGoUZEREJQTCKc2m4RvtV3QPNB++rpBQ7Lsm8Lz6qqKqqrqw97T25uLnFxcQCUlZUxY8YMTjvtNJ588kmczqPLZDU1NbjdbjweD6mpAVwgb8koqN1oFt272BO4rysiIhJo734byv5p2uPvgjHz7K2nB3r6/m3rQOH09HTS09N7dG9paSkzZ85k0qRJPPHEE0cdaHpVQ0swi+/Z30VERMQ2p/wRyv9lxtSsXWAW6EvIsLuqgAihZNC90tJSZsyYwZAhQ7jnnnuoqqqioqKCioqKI39yMHhb5vzHJNlbh4iIyJG482BEy/CN5lr44k576wmgsJjS/eabb7Jp0yY2bdrE4MGDO7xm49Oz1gLAW2/aLoUaEREJA+N+A1sXQXMdbHoYRv8MUkfbXdVxC4uemquuugrLsro8bNcaaABcifbVISIi0lOJmXDS7aZteaFkrr31BEhYhJqQ5msfahLsq0NERORonHQLJGaZ9s5XoPJ9W8sJBIWa4+Vr9redsfbVISIicjRikuHk3/rPV91qhlSEMYWa42W1CzUOl311iIiIHK3hV0LfcaZdXQzl4b2Lt0LN8bK8/rZCjYiIhBOnC8bO95+v+e+w7q1RqDle7ffOcOjbKSIiYSbnO+AeY9q7l8Oud+yt5zjoXfh4OeP9bW+DfXWIiIgcC4cTxvzSf77mv+2r5Tgp1ByvmGR/u7nOvjpERESO1ZBLIGWUaVe+C5Uf2lrOsVKoOV7t16ZRqBERkXDkdMGYO/znYdpbo1BzvJwu//o0rdsliIiIhJthl0HycNOueAN2F9tbzzFQqAmE1u0R1FMjIiLhyhnbccfuMOytUagJhNZxNQo1IiISzoZfCUk5pl32KtRssLeeo6RQEwitoUaPn0REJJy54mD0f/jPv3rIvlqOgUJNILR//BTGixaJiIiQe7V/EszXT0JTra3lHA2FmkCI62c+Wl5o8thbi4iIyGF4fRbLN1ezuKSU5Zur8foO+WU8Pg2G/dC0m2thy9PBL/IYxdhdQERIGuRvHyiFuL62lSIiItKdojXlFC5ZR7mnvu1aljuB+bPzKBib5b9x1I2w+f9M+6sH4YQ54HAEudqjp56aQEhsF2oOltpXh4iISDeK1pQzZ9HKDoEGoMJTz5xFKylaU+6/2G88pJ9h2jVfwq63g1jpsVOoCYSkwf72gZ321SEiItIFr8+icMk6uhr12XqtcMm6jo+iRt3ob3/1YG+WFzAKNYFw6OMnERGREFK8ZU+nHpr2LKDcU0/xlj3+iznfgcRs0y79B9Rt690iA0ChJhDa99To8ZOIiISYytruA0239zljYeT1pm35YNNjvVBZYCnUBEKiempERCR0ZaQkHNt9I39idvEG2PqMCTchTKEmEBIywNEykeygxtSIiEhoyR+eRpY7ge7mLzkws6Dyh6d1fCExCzLPNe26rVC1rBerPH4KNYHgcPrH1dRt0wJ8IiISUlxOB/Nn5wF0Cjat5/Nn5+FydhF7hl3ub29d1Cv1BYpCTaCknmg+Nu6F+gp7axERETlEwdgsFl4+kUx3x0dMme4EFl4+seM6Ne0Nvsi/cv7258Hb0LuFHgctvhco7jFQ/rppe9aaLjsREZEQUjA2i3PzMinesofK2noyUswjpy57aFrF9jEzobY+Y35xL/sX5FwUtJqPhnpqAsU91t/et8a+OkRERA7D5XQwdUR/LpwwiKkj+h8+0LQKk0dQCjWB0rddqPEo1IiISATJnAUJA027dAk07rO1nO4o1ARK6kn+9r619tUhIiISaM4YGHqpafsaYfuL9tbTDYWaQIntA8nDTduzRjOgREQksgy7zN/eudi+Og5DoSaQ3GPMx+b9cGC7vbWIiIgEUtpk/7YJu96C5gP21tMFhZpAaj+uZt8X9tUhIiISaA4HZH/LtL31sOsde+vpgkJNIPU7xd/e/bF9dYiIiPSGQd/2t0tfta+ObijUBFL66f521Yf21SEiItIbMs8BZ7xpl70acuNHFWoCKWmQf7Bw9QrwNtpbj4iISCDFJMPAs037wE7Y97m99RxCoSbQWntrvPWwd6W9tYiIiARaCD+CUqgJtPTp/rYeQYmISKQZ9C1/O9xDzZVXXsn777/fG7VEBoUaERGJZMlD/VsDVa8IqdWFjzrUeDweZs2axQknnMBdd91FaWlpb9QVvtwnQVw/0676MOQGUYmIiBy3zHNaGhbsXm5rKe0ddah55ZVXKC0tZc6cOTz33HMMGzaMb3zjG7z44os0NTX1Ro0AXHDBBQwZMoSEhASysrK44oorKCsr67U/75g5nDCgZVxNQzXUrLe3HhERkUAL0dm+xzSmJj09nVtuuYXVq1ezYsUKRo4cyRVXXEF2djY333wzGzduDHSdzJw5k+eff54NGzbw0ksvsXnzZr73ve8F/M8JiIx2j6AqltpXh4iISG8YEEGhplV5eTlvvvkmb775Ji6Xi29+85t88cUX5OXlcd999wWqRgBuvvlmTjvtNIYOHcq0adOYO3cuH3/88WF7hxoaGqipqelwBEXWN/ztstAaRCUiInLckrKhT65pVxeDt8HeelocdahpamripZde4tvf/jZDhw7lhRde4KabbqKsrIynnnqKt956i+eff54777yzN+oFYM+ePTzzzDNMmzaN2NjYbu9bsGABbre77cjJyem1mjroOw6SWv6sXe9A0/7g/LkiIiLB0joxxlsPe0JjCZOjDjVZWVlce+21DB06lOLiYj799FOuv/56UlNT2+6ZOXMmffv2DWSdANx+++0kJyfTv39/tm/fzuLFh98ldN68eXg8nrZjx44dAa+pSw6Hfx6/rxEq3grOnysiIhIsITjb96hDzX333UdZWRkPPfQQEyZM6PKevn37smXLliN+rblz5+JwOA57rF/vH2h76623smrVKt544w1cLhc/+tGPsA4zuyg+Pp7U1NQOR9Bkt1ucSI+gREQk0rQPNbs/sq+OdhzW4VJBL6uqqqK6uvqw9+Tm5hIXF9fp+s6dO8nJyWHZsmVMnTq1R39eTU0Nbrcbj8fT+wGn+SC81B+8ByEhE75TamZGiYiIRALLBy+lQ+MeiO8P/1ZlnlT0gp6+f8f0yp/eQ+np6aSnpx/T5/p8PsAMBg5JMYkw8BzTS1NfAXtXQdoku6sSEREJDIcT+udDeZFZwuRgmdkD0UZh0XWwYsUKHnzwQUpKSti2bRtvv/02l156KSNGjOhxL40tQnh/DBERkePWd6y/vW+NfXW0CItQk5SUxN///nfOOeccRo8ezTXXXMPJJ5/Me++9R3x8vN3lda/9/hg7XravDhERkd7gbhdqPGvtq6OFrY+femrcuHG8/fbbdpdx9JIGm6656mLYt9qk2PapVkREJJy5x/jbHvXURL5hl/vbW5+xrw4REZFAc58EtAwO1uOnKDD0++BwmfbWZ8xocRERkUgQk+xfWbhmne3vcQo1vS0hA7LON+0DO6DyA3vrERERCaTWR1DNdVC3zdZSFGqCYdgV/vbWRfbVISIiEmghNANKoSYYBl8AMX1Me/sLZp8MERGRMOT1WSzfXM3iklKWb67G2yfP/2LtBvsKI0xmP4W9mCTI+S5seQqaPGbNmiHfs7sqERGRo1K0ppzCJeso9/h/Oc/q04/5/adS4F4OB0ptrE49NcEzvN0sqC3/z746REREjkHRmnLmLFrZIdAAVOy3mLPtDoo8U+GgQk10yJgJiS3LR5f90/bBVCIiIj3l9VkULllHV5tFWi3/W1j2U7x1ZcEt7BAKNcHidMHI60zb8sHGh+2tR0REpIeKt+zp1EPTnoWT8qZ0iisSglhVZwo1wTTyWnDGmvbmxzRgWEREwkJlbc/eryrrfLauVaNQE0yJmTDkEtNuqIZtz9pbj4iISA9kpPSsByYjZjfUV/VyNd1TqAm2UTf62xseAKurJ5QiIiKhI394GlnuhNYNETpxYJEVW0V+8lpbBwsr1ARb/ymQNtm0966E3R/bW4+IiMgRuJwO5s8269EcGmxaz+dnP4rL4bN1WrdCTbA5HB17a7560L5aREREeqhgbBYLL59Iprvjo6hMdwILZ+4y69QANNXYUJ2hxffsMPT7sOo/oWE37HgBDv7JjLcREREJYQVjszg3L5PiLXuorK0nIyWB/OFpuLY8DtUtN3nrbKtPPTV2cCXAiGtN29cE6++ztx4REZEecjkdTB3RnwsnDGLqiP64nA5wJftvaFaoiT6jbgRnvGlvfMjMhhIREQlHMe1DzQHbylCosUtSNoy4xrSb62D9/baWIyIicsxikvxt9dREqbzb/YvxffVnaNxnazkiIiLHpH1PjVc9NdEpeQgMv9K0m2o0E0pERMJTjMbUCEDeXHC4THv9fdBUa289IiIiR0uhRgBIGQFDLzPtxj2wcaG99YiIiBwtZ5y/7Wu0rwzb/mTxG3MHbWsyrv+TrSPHRUREjprl9bcd9i2Bp1ATCtwn+je6rK+EDX+2tx4REZGj4Wv2t1uHVNhAoSZUjP01OFr+71i3wNZdTkVERI5Kh54ahRrpOwZyW9ataaqBNXfaW4+IiEhP+Zr8bacePwnAyYX+EeQbH4aar+ytR0REpCe8B/1tV6JtZSjUhJLELDjpVtO2mqFkrr31iIiI9IS33t9WqJE2J/2nCTcAO1+Gyg/trUdERORI1FMjXYpJhnHtxtOs+k+wLPvqEREROZL2C+613wcqyBRqQlHu1eAea9rVK2D7C/bWIyIicjj1u/zt+HTbylCoCUVOF5xyt/+85DZbl50WERE5rPahJjHTtjIUakJV1vmQeZ5p122DLzTFW0REQlR9hb+dMNC2MhRqQpXDAZMfBGe8OV9/L+z7wt6aREREutK+pyZBPTXSldQTWvaFwkzxLr4OLJ+9NYmIiBzqYPuemgzbylCoCXV5t0PKKNPevRw2/8XeekRERA7hPbiL5fvHsbj2Gyzfuh+vz55Zu2EXahoaGpgwYQIOh4OSkhK7y+l9rnjIf9h/XnK72fRSREQkBBStKWf6p7/m0q8X8B9bbuDSxz5m+h/epmhNedBrCbtQc9ttt5GdnW13GcE1cCYMu8K0G/fCyl/YW4+IiAgm0MxZtJLypv4drld46pmzaGXQg01YhZp//etfvPHGG9xzzz12lxJ8E++BuH6mvXURVCy1tx4REYlqXp9F4ZJ1mAdNjg6vtT58KlyyLqiPosIm1OzatYtrr72Wp59+mqSknq1W2NDQQE1NTYcjbCVkwIQ/+M+Lr4Om/fbVIyIiUa14yx7KPfXdvm4B5Z56irfsCVpNYRFqLMviqquu4vrrr2fy5Mk9/rwFCxbgdrvbjpycnF6sMghGXAPpp5v2/s2w6lZ76xERkahVWdt9oDmW+wLB1lAzd+5cHA7HYY/169fzwAMPUFtby7x5847q68+bNw+Px9N27Nixo5f+JkHicMJpT4Krpadq08NQ9i9bSxIRkeiUkZIQ0PsCwWFZ9u2WWFVVRXV19WHvyc3N5ZJLLmHJkiU4HP5ndl6vF5fLxQ9/+EOeeuqpHv15NTU1uN1uPB4Pqampx1W7rTY+DJ/MMe3ELPjmFxDf//CfIyIiEkBen8X0P7xNhecg1iFjasCMssl0J/Dh7WfjcnZ+/Wj09P3b1lDTU9u3b+8wHqasrIzzzz+fF198kSlTpjB48OAefZ2ICTWWBe9+E8qLzPmQ78P0Z+2tSUREok7R5zuY89fVgIXV7uFPa4RZePlECsZmHfef09P377AYUzNkyBDGjh3bdowaZRajGzFiRI8DTURxOGDKX/yzobY/B1v/Zm9NIiISdQqGeFg49C4yYzs+dcl0JwQs0ByNmKD+aRI4Sdlw6kL46Afm/JN/h4wzIWmQvXWJiEj08KylwL2cc1NXUDzwHioHXEJGSgL5w9OO+5HTsQjLUDNs2DDC4KlZ7xv6fdi5GLb9DZr2wcc/hplFpidHRESkt+1eDoDL4WPqSSfAIHt/sQ6Lx09yGJMfhMSWFZYr3oD199lbj4iIRI+qj/ztAdPsq6OFQk24i0+D057wn5fcDlXL7KtHRESiQ9N+2LvKtN1jzPuRzRRqIkHWeZDXsoaP1QwfXgL1VfbWJCIika26GCyvabcuDGszhZpIcfKdkHGWaR8sheVXgOWztyYREYlc7R89pU+3r452FGoihTMGTv8bJAw05+Wvw9q77K1JREQiV9WH/rZ6aiTgErNMsHG0/N/6xXyoeNvemkREJPL4vG0zn0jMguTh9tbTQqEm0gycCeMKTdvywbJL4UCZvTWJiEhk2bsKmmtNe8DpIbOUiEJNJBpzB2QVmHZ9pQk2vmZ7axIRkchRusTfzjzHvjoOoVATiRxOmPo0JLVsIVH5Pqy82d6aREQkcpT+w98edIF9dRxCoSZSJQyA058HZ6w5/+pB+Op/7a1JRETCX9122Fti2mmnmm17QoRCTSRLnwr5j/rPP/s5VLxlXz0iIhL+drbrpRkcOr00oFAT+XKvgpNuNW3LCx9cDDVf2VqSiIiEsRB99AQKNdFh/AIYNNu0m/bBe9+Gxr22liQiImGo0QOV75p28lDoO87Wcg6lUBMNnC6Y9oz/H1/tRtNj42uyty4REQkv5UX+945BF4bMVO5WCjXRIjYFzloC8enmfNdS+PTnYFn21iUiIuFj27P+doiNpwGFmuiSPBTOfBmcceZ808Ow4c/21iQiIuHh4C4ofdW0E7P8+w2GEIWaaJN+OuQ/5j9feTNsfbb7+0VERAC2LgKrZSHX4VeaPQdDjEJNNMr9EYz5VcuJBR//CMrfsLUkEREJYZYFXz/uP8+92r5aDkOhJlqdfCeMuNa0fU3wwb/B7mJ7axIRkdBUXQyedaadPh1SR9lbTzcUaqKVwwGn/i8M/o45b66D974JnvX21iUiIqGnQy/Nj+2r4wgUaqKZMwZO/6t/sFdDNbxzPhzYaW9dIiISOpoPwNa/mXZMMgy52N56DkOhJtq5EuDMxdBvgjk/sB3eKYCGPbaWJSIiIWLbc9Bca9pDvg+xfeyt5zAUagTi3DDjX9An15x71ppVh5sP2FuXiIjYy/LB+j/5z0f8xL5aekChRozETJj5BiQMNOe7l8P7F0LzQXvrEhER+5S9Zn7RBRgwzWyUHMIUasQvZQTMLILYVHNe8ZaCjYhINFv3R3877zYAvD6L5ZurWVxSyvLN1Xh9obMyfeitnCP26jcBZhSZAcPNtVDxJrx/EZz5CsQk2lyciIgETdVyqPrAtFNPhEGzKVpTTuGSdZR76ttuy3InMH92HgVjs2wq1E89NdJZ+lTTYxPTMhis4g344DvgrT/854mISOT4sl0vzUm3UrR2F3MWrewQaAAqPPXMWbSSojXlQS6wM4Ua6Vr6NJj5uj/YlL8O7yvYiIhEBc962LnYtBOz8Q65jMIl6+jqQVPrtcIl62x/FKVQI91Ln9axx6a8CN7/NwUbEZFI9+XdtMWV0TdRvL2uUw9NexZQ7qmneIu9y4Eo1MjhpZ9upnvHJJvz8n/BB98Fb0OH20J54JiIiBwFz5ew5SnTjk2FE66jsrZnv8z29L7eooHCcmQZ002wefcbZjuFstfMo6gzXoSYpJAfOCYiIkdh9TywvKZ90m0Qm0pGSlOPPjUjJaEXCzsy9dRIz2ScATNe69hj804BRSWbQ37gmIiI9FDlBx3G0nDizQDkD08jy52Ao5tPc2B+mc0fnhaUMrujUCM9l3Fmy6OoFAC8lR9R+PcVIT9wTEREesCyYNWt/vOT74SYJABcTgfzZ+cBdAo2refzZ+fhcnYXe4JDoUaOTsYZMOsdiO9Pcd0Yyhvd3d4aKgPHRESkB3a8CNUrTNs9BoZf1eHlgrFZLLx8Ipnujo+YMt0JLLx8YkgMN9CYGjl6aZNg1vtUPv/LHt1u98AxERE5Am8jlMzzn0/4IzhdnW4rGJvFuXmZFG/ZQ2VtPRkp5pGT3T00rRRq5Ni488iY+lvYvPWIt9o9cExERI5g08Owf7NpD5wJ2d/o9laX08HUEf2DVNjR0eMnOWb5eXlkpcbi6HJUTegMHBMRkcM4UAqrf+U/P+VucIRGz8vRCptQM2zYMBwOR4fj97//vd1lRTWX08H8C8YBjk7BpvU8FAaOiYjIYXx6o9nrDyD3x2aIQZgKm1ADcOedd1JeXt52/OxnP7O7pKjXNnAs9ZCBY7G7WXjGJgrGZNpUmYiIHNGOv8POV0w7IcP00oSxsBpTk5KSQmZmz98kGxoaaGjwr3xbU1PTG2VFvbaBY19XUbn6ETKqXyA/eS0ujw+WrYDTHgeXxtWIiISUxn2ml6bVpAcgPryHC4RVT83vf/97+vfvzymnnMLdd99Nc3PzYe9fsGABbre77cjJyQlSpdHH5XQwdWQGF/7br5h6xlW4HC2Po7b9DZaeDfWV9hYoIiIdldwOB1sWSM3+Ngy52N56AsBhWVZYrIx27733MnHiRNLS0li2bBnz5s3j6quv5t577+32c7rqqcnJycHj8ZCamhqMsqPXzn/AR5eC94A5Tx4GM/4J7jxbyxIREaDyfXjrLNOO6QPfWgfJofuLf01NDW63+4jv37aGmrlz5/KHP/zhsPd8+eWXnHjiiZ2uP/7441x33XXs37+f+Pj4Hv15Pf2mSIDsWQnvzYaDZeY8NhWmvwhZ59pbl4hINGs+AEUToWaDOZ/0Zxgd2mNUwyLUVFVVUV1dfdh7cnNziYuL63R97dq1jB07lvXr1zN69Oge/XkKNTY4UGqCzd5V5tzhgskPwQnX2VuXiEi0WvFT2PyYafefAud+1OVCe6Gkp+/ftg4UTk9PJz09/Zg+t6SkBKfTSUZGRoCrkoBKGgSz3odlP4TSf5idXz+5HjzrYOI94Iy1u0IRkeix7Tl/oHElwWlPhnygORphMftp+fLlrFixgpkzZ5KSksLy5cu5+eabufzyy+nXr5/d5cmRxPaBM/5uBqWt/5O59tWfTe/N9OchUdO+RUR63f4tUPxT//nkB8HdeXhHOAuL2U/x8fE8++yznHXWWYwZM4bf/e533HzzzTz66KN2lyY95XSZnpn8R/29M1UfmOe6VR/ZW5uISKTzNcFHP4CmlqVNhl4GuVfZWlJvCJvZT4GgMTUhYvcK+OC7cLDUnDtiYOK9eEfeQPHWvSG5SZqISLjy+iyK31xA5abXyYjdS/7Ag7i++ZmZvBEmwmJMjUSpAVPgGyvhw+9D5btgNVP09t8o/Fsm5fVJbbdluROYPzsvJLazFxEJR0Vryil8ZSXl+8cD4wHIqnQyf0QdBWPDJ9T0VFg8fpIIlJABZ78JJ91KkWcqc7bdQXl9YodbKjz1zFm0kqI15TYVKSISvorWlDNn0UrK93d8IFOx3xexP1sVasQ+zhi84/9A4e55LdtfdnzU1PqfYeGSdXh9UfOUVETkuHl9FoX/WIOFRTT9bFWoEVsVb9lDeZ2TQ/+ja2UB5Z56irfsCWpdIiLhrPjrKsprGom2n60KNWKrytr6gN4nIiJQ+cXjPbsvwn62KtSIrTJSerZ7d0btOxA9E/VERI7dpkfJqPprj27t6c/gcKFQI7bKH55Gljuhmw5ScOAjK7aK/O1Xmg0yG/cFszwRkfCy6x345Abyk9eSFVuFg65/GXRgZpjmD08Lbn29TKFGbOVyOpg/2+zcfWiwcbT87/zsR3E5fLD9OfjXBKhaFtwiRUTCgedLswaY1YzL4WP+5HLA0c3PVpg/Oy/i1gJTqBHbFYzNYuHlE8l0d+wGzXQnsPDySRR862aIdZuLddvgrTPhi/8Gn9eGakVEQlDtZnj7HGjca86zv0nB7NsO87N1YkSuAaYVhSVkeH0WxVv2dL2icN02sylm+y0V0qfDaU9Aykh7ChYRCQUHdsKb083PSYC0SXDO220rBh/2Z2uY6On7t0KNhA9fM6z5Laz9b7B85porCSb8HkbdAA51PIpIlDm4y/Re135lzt1jYNZ7EN/f3roCrKfv33oXkPDhjIGTfwPnvAvJw8w17wH47OewdKbpfm3H67NYvrmaxSWlLN9cHXGLTIlIlGvYA++c5w80fUaaldojLNAcDfXUSHhq2g8lt8PG//Vfa9drU7R2F4VL1lHu8a/BoL2kRCRiNNXC27OguticJ+XAuR9A8lB76+olevzUBYWaCFTxNqy4Buq2tl0q4hrmfP6dThMZW58gR+oAORGJEk218N63ofJ9c54wEGZ9AKkn2FtXL9LjJ4kOmWfDN7+AE/4dAK/lpPDL6S37nXQUyfudiEiUaKiGpef4A01cmnnkFMGB5mgo1Ej4i+0Dpz4EZy+l2HsO5U3pRNt+JyISBQ6UmkHBez4x53H9YObr0HecvXWFkBi7CxAJmMyzqRwzDNZ9ecRbI22/ExGJcLWb4O1z/Y/aEzLh7DcUaA6hUCMRJaOvu2f3Rdh+JyISwfZ+bmY51e8y58nD4Zy3oE+uvXWFIIUaiSite0lVeOq73PHEgY/MhAPkZ+wHonfao4iEpk4L5fVZj+v9b0PTPnODe6x55JSUbWudoUqhRiJK615ScxatxAEdgo0DH+Bg/sD/wfXaj2Hcb+DEm8AZa0utIiLtFa0p77wURexu5mefRIF7OfSfAjNeg/jI2oQykDRQWCJOt3tJJcPCkf9rfjh4D0DJbfCviVD5oU2ViogYRWvKmbNoZYdAA1DRlMacbXdQ5LgOzn5LgeYItE6NRKwu9ztp9sDqX7Us2tfun/7Qy2DCAkgeYlu9IhKdvD6L6X94u1OgaeXAItOdwIe3nxN2ezYFitapkajncjqYOqI/F04YxNQR/c0Pg7i+cOqDcH6x2fSt1ba/wqujYfWvzWrFIiJBUrxlT7eBBsDCQbmnQUtR9IBCjUSn/pPhvBUw+SH/Pineelj7W1hyAmx+HHxee2sUkajQ0yUmtBTFkSnUSPRyumDUv8PsjXDiLf4Bw/UVZuuF1yfDrnfsrVFEIl5GXc/G9WkpiiNTqBGJ6wcT/wTfWgeDv+O/vrcElp4N718ENRvtqk5EIlXzAVjxU/K3XkZWbFXLDM3OHJgNefOHa5DwkSjUiLRKGQln/h3OeRf6TfRf37kY/pkHn90E9ZU2FSciEWXfGnj9VNj8GC6Hj/nZjwKOThu8tJ7Pn50XtYOEj4ZCjcihBp4FBZ/AaU9CYssCV1YzbPgf+EculNwBjXttLVFEwpRlwcZHTKDxrDPXXEkUnPcTFv5wUuelKNwJLLx8IgVjs2woNvxoSrfI4TTXwbq74cs/gveg/3qsG078hVm8LzbFtvJEJIw07oUV18KOl/zX+p4Mpz8L7pOAbpaiUA9Nj9+/FWpEeuJgBay9CzY9Ar7GtsveuHSK+xdSmVpARt9U/QASka5VLYOPLoUD2/3XRt0Ip9wNLg0APhKFmi4o1Mhxq9sOa+6Er5+kaF8+hWU/pbwpve3lrNR45l8wRl3FImI0HzQ/M768G6yWZSLi0uC0x2HwhfbWFkYUarqgUCOBUlS8kjl/L2tZk9jfM+NoubLwsvEUnJxjS20iEiIqP4AVP4Har/zX0s+Aac9Asn4+HA2tKCzSS7w+i8Kle7FwwCFzFayWbTQLX3wX77r7zZgcEYkuTTXwyb/DW2f6A40zDk7+LZzztgJNL1KoETlKR17S3El5YxrFHz4Oi4fCF3dCg5Y3F4kKpa/BP8fAxoX+awOmwjdKYOwvwRljW2nRQKFG5Cj1eEnzpn7QUA1fzDfhZtWtcLC8l6sTkd7g9Vks31zN4pJSlm+uxus7ZORG/W5Ydjm89y04sNNci0mGSX+GWR+0zW6S3qXIKHKUerpUeUbOJKj7yAwObN4PX94DG/4MuVfBSbdByojeLVREAqJoTTmFS9Z16KHNcicwf3YeBXnpsOlR+PzX0NiuRzbzPMh/BPoMC37BUSysemr++c9/MmXKFBITE+nXrx8XXXSR3SVJFMofnkaWO6HTyp+t2pY0//afzL5SJ8wBZ7x50ddofgC+Ogo+/AFULTeLcYlISCpaU86cRSs7PXKu8NQzZ9FKihb9CD69wR9o4vqZhTtnFinQ2CBsQs1LL73EFVdcwdVXX83q1av56KOPuOyyy+wuS6KQy+lg/uw84NBhwl0sad5nOJz6v3DhVtM7E9OyUJ/lg+3PwZvT4PV82LIIvA3B+iuISA94fRaFS9bR1a8d5pqPwo3n4bVa3kqHXW72kMu9Ehxar8oOYTGlu7m5mWHDhlFYWMg111xzzF9HU7olkA7bJd3dOjWN++Crh8yWCw1VHV9LGGh6dUZeB4mZvVe4iPTI8s3VXPrYx0e872/jH2fqrF9A+rQgVBWdevr+HRZjalauXElpaSlOp5NTTjmFiooKJkyYwN13383YsWO7/byGhgYaGvy//dbU1ASjXIkSBWOzODcv8+iWNI/ra2ZAnPQL2PacCTd7V5nX6nfBF7+Btb+DIT+A0T+H/pOD8VcRkS70eFLAifdDuqZph4KwePz09ddfA/Cb3/yGX/3qV7z66qv069ePGTNmsGdP91NlFyxYgNvtbjtycvSPTgLL5XQwdUR/LpwwiKkj+vd8iwRXgumiLvjMzIwYcjE4XOY1XxNsfdpsePfG6bD1WT2aErFBRnLP3iIzUpN6uRLpKVtDzdy5c3E4HIc91q9fj8/nA+CXv/wl3/3ud5k0aRJPPPEEDoeDF154oduvP2/ePDweT9uxY8eOYP3VRHrG4YCM6TD9ebjga8iba5ZQb7V7GSy7FF4ZDJ/dAvvW2lerSLTwNsLGR8hfcxpZsVU48HV5W9ukgOFpXb4uwWfr46df/OIXXHXVVYe9Jzc3l/Jys7ZHXl5e2/X4+Hhyc3PZvn17d59KfHw88fHxAalVpNclD4EJC2Dsr2HrX82jKc8a81rDbthwnzkGTIURP4Ehl0BsH3trFokkvib4+knzCLhuGy5gfvajzNl2Bw6slhXDjU6TAiQk2Bpq0tPTSU9PP+J9kyZNIj4+ng0bNjB9+nQAmpqa2Lp1K0OHDu3tMkWCKyYJRv4ERlwDle/Bpsdgx0vga3kEtXu5OT77Dxh6qQk4/U/VbAuRY+VtNI981/wW6rZ2eKngpP4sPKUfhe/Ud5gUkHmkSQFii7CY/QRw00038eKLL/L4448zdOhQ7r77bpYsWcL69evp169fj76GZj9J2GrYA1ufgc2Pwb4vOr/uHmvCzbDLIOHIvyiICGYV4E2PwFcPQn1Fx9eyvgHjfgMD8gEzvfuoJgVIQEXcLt1NTU3MmzePp59+moMHDzJlyhTuv/9+xowZ0+OvoVAjYc+yYM+nsPn/zCOq5v0dX3e4zEqmwy6DwRdCbIo9dYqEMs+XsOF+2PL/wHvIDKes81vCzGl2VCbdiLhQEwgKNRJRmvbD9hdMwNm9rPPrrkQYNNsEnKwCcGl8mUQxy4KKt2D9vVBe1PE1hxMGfwdO/AWkT7WnPjkshZouKNRIxPKsg6+fgm3PwgEzeN5rOSmuG0NlUz8yEpvIz8vDNewyyDgLnC6bCxYJkqZa89/Fhj/7B963ikkxj21H/8ys/i0hS6GmCwo1EvEsH1Qto2jZ2xR+NpzyJv9U06zYKuZnP0pB5lYY8n3I+TcYME0BRyKPZUHVR/D1X2Db8+A90PH15GFmccsR10Cs3gvCgUJNFxRqJBq0bsB36H/YZq0NBwuH3kWBe7m5GD/APKIafBFkngsxiUGuViSADpabcTKbH4farzq/PmAanHiLGW/mDIsF9aWFQk0XFGok0nl9FtP/8HanHYVbObDIjN3Nhydeg8txyIJiriQzSHLwRTDoWxDfv/cLFjkGHWYi9YkhP+5jXFv+AmWvgeXteHOs24wrG3ENpE2yp2A5bhG195OI9Ezxlj3dBhoACwflTekUD3uGqb4XzYDJ5jrzovcA7HzZHA4XpJ9hAs7gCzTeQEJG0ZpyCv+xlvIa/9YhWbF7mZ+9mwJ3u0AzcCbkXgM53zFrP0lUUKgRiSA93oAv+QyY8AMznbViKex8BXYu9u8cbnmh8l1zrLwJUk4wU8WzzoOBMzQOQYLP10TRh28x5zVvy6NV/xoxFU39mbPtDhaOeoSCSRNgxNXQJ9emQsVOCjUiESQjJeHo7nMlmEdNg74Fpz4M1R+bgLPjZdi/2f8JtRvNsfEhcMSYrRqyzjNBJ22SBhtL7/A2wq6lsP1FvDsWU/j5PVgMoH2gAbBw4gAKd93EuePO1qJ4UUyhRiSC5A9PI8udQIWnvtNAYTBvBZndbcDndEH66eaY8EczTbx0MZS/DlXLwGo291nNUPWBOT7/NcT1g8xZJuBknmNmlmjLBjlWB0rNv7nyIih/E5r2AVC8fxzlTd2vlm0B5Z56irfsYeoIjQeLVgo1IhHE5XQwf3YecxatxAEdgs1RbcDncEDfMeYYc4dZ62PXu1DxBpS/0XFmSeNeswjg9hfMeeIgSJ8OGWeYj+6xeHFqiXnpmrfBTL8uLzJHV9uAAJW+nu2x1NNHsBKZFGpEIkzB2CwWXj6RwiXrArcBX2wKDJ5tDoD9W6HiTRNwKt5q+20agIOlsP05cwBF+8+hsPQnlDf4t2zICpHNAENhP59QqCGoLMuE4oqlJsTsets/WP1Qcf3MHkxDLibjYD48vuqIX76nj2AlMmlKt0iECtqbpc9r9qMqf8M8ktq9rO1NqsgzlTnb7ug0sNPRcmXhefUUTB4PKSPMUvVBVLSmvFPwC3bYCoUaep23EfaugqoP/UfD7m5udpgd57MKzNH/1Lb1ZFqXKzjSo9UPb9eYmkikdWq6oFAjEgS+Zti3Gm/FB0x/Ppvy+iQOHdgJZjHAzNhqs2ZOXAr0O8UMOm49Ukb2WtDpfoFCY+HlE3s9VIRCDb2i0QO7P/YHmOoV4D3Y/f0JA836SFkFZgHIhAHd3tr6PYOuH62G7fdMjkjr1IiIPZwxkDaJ4r3DKK//uNvbLJxmzZy6MUx1fOGfQt4qNtUfdPqdAu4xkHrica967PVZFC5Z1+Vv+xbmDbJwyTrOzcvstd/4Q6GGgKivhL0lLccq87FmA3T5N2sR27dlQPp0E2b6je9xeO2VR6sSURRqRKRX9HjNnJRZkLDHjMVpr6kGKt8zRyuHE5JzzQBm9xhIPckEndTRZtxPDxx5gcLen0UTCjUcFZ8X6rZ0DC97S+Bg2ZE/N3mYCTDp002YcecdVw9cwdgszs3LjK5xSNJjCjUi0it6vGZO/i9hxL1wcBfs+cwce1s+HtjZ8WbLB/s3mWPn4o6vJQ5qCTijIHk49Blm3lCTh5k9rlqmmfc4bPXiLJpQqKETyzIhpXajGcjbujZRzVdmzSJf45G/hjMO3GMhfZo/xCQNDnipLqcjNMKehByFGhHpFUe9Zk7iQBj0TXO0qq+EPSth3+fgWdtyfNl512UwPT0HS81ibYdyJUHyUEgeRsbBicDUI9bfm7NojnqRxEDwNkJ9uVkH5mCp/+P+rf4A09X3tTtx/aDfBOg7wXxMO8WESmds4GoWOUoKNSLSKwKyZk5CBmQXmKOV5YO6rbBvLdRugJr1LceG7mfVeA9AzZdQ8yX51utkxf6Fiqb+WHR+DOLAIjOhjvzdd0DdQIjPMHUkZJjNEWP6mEddsSkmLB3DQoPHtUgimF6V5jozlb5xLzQe+nFv5wDTugXG0XLGm0HbKaPM+Jd+E8yRNESLLErI0ewnEelVQZ22XL/b9DjUbYMD20wvRF27w2tqaJ1qDlaHYOPABzhYOPQuCtzLj/znOZwm5MS0hJy2j33MgGmHy2wr4XD5j5brRWXZzFk+ATAbjfpraJnuPuF1CjLWm8XpfA3m8U9znT+8tK7wHAiOGLNpacoos89XygnmMV7KCZCUE/Tp9iKH0pTuLijUiNgjJBaYsyzzOKsl4BSt20NhcTrlB/2zqbJiq5mf/XDPAk0AFHmmUlj20w7L/2fFVjE/+9HA1eCIgcQsM+YoKbvl4yD/x6Qc82hOj40khCnUdEGhRkTa6zJsWU3mUU39LhOCWj82VJkZWU210FwLzfv97aaWo3k/h53O3FUNlpPiujFUNvUjI3Yv+clrcTl8nW90xoMr0Yxlievb8WNsF+eJA01wSchQT4uEPa1TIyJyBF3Poolr6cEYdPRf0PJB8wHzmMhqBstrPvq8LW2v/1pL22X5mOqKNzOHnC0f25+74lseYWn8isiRKNSIiASKwwmxfcwhIkGnPkkRERGJCAo1IiIiEhEUakRERCQiKNSIiIhIRFCoERERkYigUCMiIiIRQaFGREREIoJCjYiIiEQEhRoRERGJCFG1onDrNlc1NTU2VyIiIiI91fq+faTtKqMq1NTW1gKQk5NjcyUiIiJytGpra3G73d2+HlW7dPt8PsrKykhJScER5pvD1dTUkJOTw44dO6J6x3F9H/Q9aKXvg74HrfR9iLzvgWVZ1NbWkp2djdPZ/ciZqOqpcTqdDB482O4yAio1NTUi/sEeL30f9D1ope+Dvget9H2IrO/B4XpoWmmgsIiIiEQEhRoRERGJCAo1YSo+Pp758+cTHx9vdym20vdB34NW+j7oe9BK34fo/R5E1UBhERERiVzqqREREZGIoFAjIiIiEUGhRkRERCKCQo2IiIhEBIWaCNPQ0MCECRNwOByUlJTYXU5QXXDBBQwZMoSEhASysrK44oorKCsrs7usoNq6dSvXXHMNw4cPJzExkREjRjB//nwaGxvtLi2ofve73zFt2jSSkpLo27ev3eUEzUMPPcSwYcNISEhgypQpFBcX211SUL3//vvMnj2b7OxsHA4Hr7zyit0lBd2CBQs49dRTSUlJISMjg4suuogNGzbYXVbQKNREmNtuu43s7Gy7y7DFzJkzef7559mwYQMvvfQSmzdv5nvf+57dZQXV+vXr8fl8PPLII6xdu5b77ruPhx9+mDvuuMPu0oKqsbGRiy++mDlz5thdStA899xz3HLLLcyfP5+VK1cyfvx4zj//fCorK+0uLWjq6uoYP348Dz30kN2l2Oa9997jhhtu4OOPP+bNN9+kqamJ8847j7q6OrtLCw5LIsZrr71mnXjiidbatWstwFq1apXdJdlq8eLFlsPhsBobG+0uxVZ//OMfreHDh9tdhi2eeOIJy+12211GUOTn51s33HBD27nX67Wys7OtBQsW2FiVfQDr5ZdftrsM21VWVlqA9d5779ldSlCopyZC7Nq1i2uvvZann36apKQku8ux3Z49e3jmmWeYNm0asbGxdpdjK4/HQ1pamt1lSC9qbGzks88+Y9asWW3XnE4ns2bNYvny5TZWJnbzeDwAUfMzQKEmAliWxVVXXcX111/P5MmT7S7HVrfffjvJycn079+f7du3s3jxYrtLstWmTZt44IEHuO666+wuRXrR7t278Xq9DBw4sMP1gQMHUlFRYVNVYjefz8dNN93E6aefztixY+0uJygUakLY3LlzcTgchz3Wr1/PAw88QG1tLfPmzbO75IDr6feg1a233sqqVat44403cLlc/OhHP8KKgEWzj/b7AFBaWkpBQQEXX3wx1157rU2VB86xfA9EotkNN9zAmjVrePbZZ+0uJWi0TUIIq6qqorq6+rD35Obmcskll7BkyRIcDkfbda/Xi8vl4oc//CFPPfVUb5faa3r6PYiLi+t0fefOneTk5LBs2TKmTp3aWyUGxdF+H8rKypgxYwannXYaTz75JE5n+P/+ciz/Fp588kluuukm9u3b18vV2auxsZGkpCRefPFFLrroorbrV155Jfv27YvKHkuHw8HLL7/c4fsRTW688UYWL17M+++/z/Dhw+0uJ2hi7C5Aupeenk56evoR7/vzn//Mb3/727bzsrIyzj//fJ577jmmTJnSmyX2up5+D7ri8/kAM8093B3N96G0tJSZM2cyadIknnjiiYgINHB8/xYiXVxcHJMmTWLp0qVtb+I+n4+lS5dy44032lucBJVlWfzsZz/j5Zdf5t13342qQAMKNRFhyJAhHc779OkDwIgRIxg8eLAdJQXdihUr+OSTT5g+fTr9+vVj8+bN/PrXv2bEiBFh30tzNEpLS5kxYwZDhw7lnnvuoaqqqu21zMxMGysLru3bt7Nnzx62b9+O1+ttW7Np5MiRbf99RJpbbrmFK6+8ksmTJ5Ofn8/9999PXV0dV199td2lBc3+/fvZtGlT2/mWLVsoKSkhLS2t08/JSHXDDTfw17/+lcWLF5OSktI2psrtdpOYmGhzdUFg69wr6RVbtmyJuindn3/+uTVz5kwrLS3Nio+Pt4YNG2Zdf/311s6dO+0uLaieeOIJC+jyiCZXXnlll9+Dd955x+7SetUDDzxgDRkyxIqLi7Py8/Otjz/+2O6Sguqdd97p8v/3K6+80u7Sgqa7//6feOIJu0sLCo2pERERkYgQGQ/bRUREJOop1IiIiEhEUKgRERGRiKBQIyIiIhFBoUZEREQigkKNiIiIRASFGhEREYkICjUiIiISERRqREREJCIo1IiIiEhEUKgRERGRiKBQIyJhq6qqiszMTO666662a8uWLSMuLo6lS5faWJmI2EEbWopIWHvttde46KKLWLZsGaNHj2bChAlceOGF3HvvvXaXJiJBplAjImHvhhtu4K233mLy5Ml88cUXfPLJJ8THx9tdlogEmUKNiIS9gwcPMnbsWHbs2MFnn33GuHHj7C5JRGygMTUiEvY2b95MWVkZPp+PrVu32l2OiNhEPTUiEtYaGxvJz89nwoQJjB49mvvvv58vvviCjIwMu0sTkSBTqBGRsHbrrbfy4osvsnr1avr06cNZZ52F2+3m1Vdftbs0EQkyPX4SkbD17rvvcv/99/P000+TmpqK0+nk6aef5oMPPmDhwoV2lyciQaaeGhEREYkI6qkRERGRiKBQIyIiIhFBoUZEREQigkKNiIiIRASFGhEREYkICjUiIiISERRqREREJCIo1IiIiEhEUKgRERGRiKBQIyIiIhFBoUZEREQiwv8HS6Rmdh5DvFcAAAAASUVORK5CYII=", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGwCAYAAABRgJRuAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVm1JREFUeJzt3Xl8VPW9//HXzGQPyUAgIQmEJSBIQEHAIIgKipouqL2ttm5Va61ytb0uVcG2l6a3lbba6q16cWldfmLrWqVYjQvugkSFqIAgIGsWEgJMQiDbzPn98U0yCUkgwGTOLO/n43Ga7zlzEj6kmHnne76Lw7IsCxEREZEw57S7ABEREZFAUKgRERGRiKBQIyIiIhFBoUZEREQigkKNiIiIRASFGhEREYkICjUiIiISEWLsLiCYfD4fZWVlpKSk4HA47C5HREREesCyLGpra8nOzsbp7L4/JqpCTVlZGTk5OXaXISIiIkdh+/btDB48uNvXoyrUpKSkAOabkpqaanM1IiIi0hM1NTXk5OS0vY93J6pCTesjp9TUVIUaERGRMHO4oSMaKCwiIiIRQaFGREREIoJCjYiIiEQEhRoRERGJCAo1IiIiEhEUakRERCQiKNSIiIhIRFCoERERkYgQVYvviUQly4LmWmj0gK8BvA3gazRtX2PLeft2I1jN4IgBh8t/ONuft2u74iEmBWJT/B9d8Xb/rUUkCinUiIQTywf1VXCgDA6Umo8N1dC4p+XYaz427fWfN+01nxdMzliI6dM57MT1g4QMiM+AxIHmY0IGJAw0H2P6gDabFZGjpFAjEkoaqqHmK9j3tQkt+0v94WV/KdSXg6/J7ioPz9fkD1pHwpXYLvRkQfIw6DPMfEweaj7G9VPwEZEuKdSIBFtTDdRugJoNUPuVabd+PNIQ0B2HC+L6Qmw/8zGuH8SmmtDgjDOPh5xx4Ixvd96u7XCB5TWHr9nftrzm0VT7c289NNWaR1ztP7a2m/dBc13P6vYegLqt5uhOTEq7oNNy9BkOqaOhzwg9+hKJYgo1Ir3F1wQ162BPScuxCjxroX7n0X29+HRIzIakQZA4qOVjtunZiOsHsX3Nx7i+ofcYx+cFbx007DZ///pKaKj0t9s+trQbdgFW11+ruRb2fmGOgzmckJwLqceD+/iWj2PBnWdCnYhENIUakUBo2gd7P+sYYPauNgNweyppCKQc13KMhKQcf4BJzArvHginC5ypJlj0GXb4+33NJtzUbYF9W8zH1mPfFti/zQxoPpjlg30bzVH2csfXknJMwOk7DvpOgLRJkDrKBCERiQgKNSJHyrKgbjNUfgBVLUfNOrrtWWgvYaB5TNIWXkaZj31GQExir5ceNpwxJtAlDYL0Uzu/bvngQIU/6NRuhJr15v+HmnXg3d/5c/ZvN0d5kf9aTB/od5IJOK1HyigTwkQk7CjUiByOr9n0wrQPMfUVh/kkhwkv/SaYN81+E6DveDPjR46dwwlJ2eZIn9bxNctnBlXXrDOP+zxr/EeTp+O9zfug6n1ztIpJ9gedfhNhwFTTcxZKj/NEpEsOy7J68OtlZKipqcHtduPxeEhN1fN16YblM4+Qyotg59uwa/mhB7o6YlreBCe2hJgJ0PcE8+YoocOyzCyyvV/A7k/9x/5th//chIGQPh3ST4OM6SagOvU7oUiw9PT9W6FGBKB+F1S8DmVFUPGaGbDandhUGDDNPBZJnw798yEmKXi1SmDVV8HulbCnXdA51OwrMI+tBkw1//9nnAb9p+jfgEgvUqjpgkKNtPE1Q3Wx6Y0pK4Ldn9DtmJjEQeaNK326OdzjNOYi0jVUm6BTXQxVH8KuD81U/O44YszjqoEzIescE3rDeWC3SIhRqOmCQk2U8zZCxZuw/TnYsbj7NWFikmHgWZBdAJnnQJ9cjaeIdj4veFZD5fst46reN4+yuhOTDBkzTMDJPMeMr9K/IZGjplDTBYWaKORtgIo3YFtLkDl4oGirvidCVoEJMgNOBVdccOuU8GJZZtZV1Qf+oFPzZff3J+X4A07mWRDfP2ilikQChZouKNRECW89lLcEmdJ/dR1kYvpA9jcg6xuQda6ZRSNyLOorTU9g+etmfNaB8m5udEDaZBj0LRh8gQnU6sUROSSFmi4o1EQwyzKzlDb9DbY/3/X4h5gUGHweDLnQBBlXQvDrlOhgWWYKeWvAqXzXhO2uJA8z4Wbw+WbMlmZViXSiUNMFhZoIdKACNj8JXz/asgDeQWJTYdD5LUHmbAUZsYe33jyiKn/dHHs/6/q+uDQYNNuEnKxzNKNKpIVCTRcUaiKErxnKXjFBpvRls6liezEpkPMdE2Qyz9YsFAk9+0vNo9HtL8HOt8wmoQdzJZgxOIMvgJwLzL5eIlFKoaYLCjVhrnYTbHwYNv+/rlf0zTgDcn8EQ76rhe8kfDTuhbJXYcdLJqw37+t8jzPWjP8adonpyVEPjkQZhZouKNSEIcsy02fX3WNmLx28lkxiNuReCblXmaXsRcKZt8H03OxYbI6uwntMsum9GXqJeaTqjA16mSLBplDTBYWaMOJthG3PmjCzZ2XH15yxMOg80yuTdY4GVkpksnywa4UZ+L71H13PporvDzkXmh6c9FO147hELIWaLijUhIGG3bDxIfjq/s6LmyVmw6gbYMSPISHdnvpE7ODzQtV7sOXvsO15aNrb+Z6kHBNucn8EqaOCXqJIb1Ko6YJCTQjbtwW+/CN8/Th4D3R8LW0SjL7JDPzVongS7bwNUP6aCTil/+r83wtAxukm/Od8V+NvJCIo1HRBoSYE1W2F1b+Drx87aAaIw6zbcfzNZu0OLU4m0lnTPjP2ZuvfTdA5eCZgrBuGXWoCTtpJ9tQoEgAKNV1QqAkhddtgzZ1mWravyX89po/pPh/9M0gZYV99IuGmvsqs2bTpr11v2dBvIoz8MQy9GOL6Br08kWOhUNMFhZoQsH8HrFkAmx7pGGZiU2H0jXD8jVqPQ+RYtK2u/VfY+gx493d83ZVoHuWOvA4GnKJeUAkLCjVdUKix0YFyE2Y2PgS+Rv/1mBQY/V9w/E0Qn2ZffSKRqKkGtj4NG/8Kuz/u/HrayaZXdMhFGq8mIU2hpgsKNTZoroMv74a1f+z4G2NMMoz6GYy5RTsWiwTDns/M3mibn+w8eyohE467DkZeC4mZtpQncigKNV1QqAkiy2dW/v3sFx2nZruSYPRP4fifQ8IA++oTiVbNB2DbM7D+f2FPScfXnLEw5Aem96b/ZFvKE+mKQk0XFGqCpOpD+ORnHRfNc8TAcXNg3C8hIcO+2kTEsCyzyeb6v8COf5pfRNobMM2Em5z/0KrFYruevn9rKVYJnP2lUHI7bHmq4/XB58OEP2pBMJFQ4nBAxmnmqNsGG/4PNj4CjbvN67uWmSMpB8b83EwL15o3EuLUUyPHztdktjNY/RszhqZV3/Ew6R4YONO+2kSk55r3m19K1v8FPKs7vhY/wMxQHHW9poRL0PX0/TvsNgp54IEHGDZsGAkJCUyZMoXi4mK7S4puu1fBa/mmh6Y10MSlwckLoeBTBRqRcBKTBCOvgW9+Dme9Bdnf9r/WsAs+/yW8NARK5sKBLjbbFLFZWIWaZ555hptvvpn58+ezcuVKxo8fz7nnnktlZaXdpUUfbz2UzIPXTvYPNnQ44bjrYfYGM5PC6bK1RBE5Sg6H+YVkxhL4xmdmwb7WzTKba2HtH2DxMPj4P2HfZltLFWkvrB4/TZkyhZNPPpn7778fAJ/PR05ODj/96U+ZO3dup/sbGhpoaGhoO6+pqSEnJ0ePn45V5Qew4mqo/cp/re8JMOVv0P9k++oSkd5Tuwm+vMtsadJ+rSmHy4SevLnQd6x99UlEi7jHT42NjXz66afMmjWr7ZrT6WTWrFksX768y89ZsGABbre77cjJyQlWuZGpqRY+vgHePM0faJyxcEIhnPuJAo1IJEsZAfkPwnmbzcDhmGRz3fLClkXwygnw4aVQs8HeOiWqhU2o2bVrF16vl4EDB3a4PnDgQCoqun62O2/ePDweT9uxffv2YJQamSqWwr/HwoYH/Nf6T4GCVXDCf2s1UpFokZQNJ90F528zv9DEta4EbuHd8jTLn7mQxc//muVr1uL1hc2DAIkQET2lOz4+nvj4eLvLCG++Jvj8v80zdFp+QLmSYPzvYNRPNW5GJFrFp5lfaMbcAhsWUvTBGxRuu5jypnTz+iebyUr8kvmz8yiYOMbeWiVqhE2oGTBgAC6Xi507d3a4vnPnTjIztax3r9i3BT68GKo/8l8beCZMeQT65NpWloiEkJhkiryXMmfTGCw69sxUHIhjzrObWFi6hIJzrtaWKNLrwubxU1xcHJMmTWLp0qVt13w+H0uXLmXq1Kk2Vhahtj0Hr07wBxpHjFlA78w3FGhEpI3XZ1G4ZG1LnOm447eFE7AoLM7A+1IufP5rs8mmSC8Jm1ADcPPNN/PII4/wxBNP8OWXXzJnzhzq6uq46qqr7C4tcjTvhxU/gQ8ugiaPudYnF87+EPJu9U/rFBEBijfvptxT3+3rFk7Km9Ip9gyF1YXwrxGwYSH4moNYpUSLsHn8BPD973+fqqoq/vu//5uKigomTJhAUVFRp8HDcpT2roYPvw+etf5rQ38AJz8IcW776hKRkFVZ232g6XBfc8sGtg27zPo2X90PJ/0Jsgt6sTqJNmG1Ts2x0jYJh7D9RVh+uX9VYFcSTL4Pcq8yC3GJiHRh+aZqLn7ko8Pe94/LcpjqWQBbn+74QlYBnHS31riRQ4q4dWqkl1gWrP4tvP8f/kDT90Qo+ARG/EiBRkQOKX94GlnuBLr7SeEAstwJ5OedAKf+A85eZpaDaFVeBK+ON7039VXBKFkimEJNNGveD8sugc9/5b829GI45yNwawqmiByey+lg/uw84OBhwv7z+bPzcDlbztKnwjnLYdrfIWmIuWZ5zTibJSNh7V3gbUDkaCjURKv9pfDm6e26gh0w/k6Y9hTEJNpamoiEl4JxWSy8bCKZ7oQO1zPdCSy8bCIF47I6foLDAcMuhm+vM2texfQx15tqoOQ2eHkM7PhXkKqXSKIxNdFo1wp47wKob1mJOaaPCTODz7O1LBEJb16fRfHm3VTW1pORkkD+8DR/D82hHKgwPcab/gbt17oZdB5M/gskD+21miU89PT9W6Em2mx9Fpb/EHwt3bvJw+CMf5kNKUVE7LTnM1h5E+x823/NlWRWLj7+ZrPXnEQlDRSWzjY8BB/+wB9oMs6Acz9WoBGR0NBvPJy5FKb9AxJaVor37oeSuWYx0Mr3bC1PQp9CTbRY83v4+DraunZzfwQzX4eEAbaWJSLSgcMBw35gxtuMusG/4KdnLbx5Bnx0lWZJSbcUaiKdZcGq2+Gzef5rY26DKX/VztoiErri3GatrHOLIW2y//rXj8PLo2HjI2D5bCtPQpNCTSTzeaH4J/DlH/3XJvweTvqD1p8RkfCQNsksMzH5AYhtWdm8cY/52fbm6VDzlb31SUhRqIlU3gYzfmbTX1suOCD/Ici73dayRESOmNMFo/7TPJIadqn/etWHZuG+dfeYX+Ik6inURCJvA7z3Hdj+vDl3xpqVPEf+xN66RESORWImTFtkBhP3GWGueeth5c2w9Az12ohCTcTxNZkemvJXzbkrEU7/Fwz9vr11iYgESuaZ8M3PYNTP/NfUayMo1EQWy2dmBux4yZy7kmDma9oFV0QiT0wyTP5fmPVuN702G+ytT2yhUBMpLAs+ngNbnjLnzjg4YzFknGZvXSIivSnj9G56bU5Ur00UUqiJBJYFq34OGx825w4XTH8OMmfZW5eISDAcstdmJtRttbc+CRqFmkjwRSGs+3PLiQOmLtI+TiISfbrstXkfXhkP256zry4JGoWacLfuHlhd6D+f8ohZjVNEJBq177VJHmauNXngg4tgxY+hua7Tp3h9Fss3VbO4pJTlm6rx+qJmS8SIE2N3AXIMdvwLVt7iP594L4y42rZyRERCRsbp8I0Ssz3M1qfNtU1/g6oPzN5SaScBULS6nMIlayn31Ld9apY7gfmz8ygYl2VD4XIs1FMTrvZ8BssuoW0vp3Hz4fj/srUkEZGQEueGaX+HUx43PTgANevh9VNg3T0UfVHGnEUrOwQagApPPXMWraRodXnwa5ZjolATjg5UwLuz/d2oQ38AJ8y3tyYRkVDkcEDuFVCw0my5AOBrxPvpzyl8/j26etDUeq1wyVo9igozCjXhxltvVgvev92c98+HKY9qLycRkUNJHQVnL4MxPweguG4s5Q0p3d5uAeWeeoo37w5SgRIICjXhxLLgo6uh+iNznjQYTn8JYhJtLUtEJCy44uCku2Dma1Q6hvXoUypr6w9/k4QMDRQOJ2vuhK1/N21XEpyxBBI1kE1E5IhknUPG9FzYtP6wt2akJAShIAkU9dSEi7LX4PNf+s+nPQX9JthWjohIOMsfPYIsdwKOLkfVgAMzCyp/eFpwC5NjolATDg6Uw/LL/efj74ScC2wrR0Qk3LmcDubPzgMcHDwi0YEPsJg/Ow+XU+MVw4lCTaizfLDscmioMufZ34K8ufbWJCISAQrGZbHwsolkujs+YsqMrWbh0DspqP89+Jpsqk6OhsbUhLq1f4CdS007MRtOeUwznUREAqRgXBZn52VSvHk3lXtryChbSH7Nn3A5fPDVctizEqY/q/GLYcJhWVbUTMKvqanB7Xbj8XhITU21u5zDq1oGb54OlhdwwFlvwcAZdlclIhK5LAs2PQKf/BR8jeZaQqaZaTpgiq2lRbOevn/r8VOoatwDH17cEmiAcb9SoBER6W0OB4z8Ccx63yybAVBfAUtnaFPMMKBQE4osC1ZcA/u3mfP000yoERGR4BiQDwWfmj2kwCx8+sFFsGaB+RktIUmhJhRtew62v2Dacf3M9G2nhj+JiARVQgbMfAOGX+G/9tkd8NFV4G20ry7plkJNqGnYDZ/+1H+e/zAk59hXj4hINHPFmQka4+/0X9v8BLx9DjRU21eXdEmhJtSsuhXqK0178AWQ811byxERiXoOB4ydZ2ZBuVqmf1e+C69PhZoN9tYmHSjUhJKdb8PXj5p2bCpMvl/Tt0VEQsWQC+GsdyBhoDmv3QCvnwI737WzKmlHoSZUNB+AFT/xn0/4PSQNsq8eERHpbMAUOHcFuMeZ88bd8PbZsOXv9tYlgEJN6Fj9P7Bvo2mnnwojr7W3HhER6VryUDjnQ8gqMOe+Jlh2Kaz730N+mtdnsXxTNYtLSlm+qRqvT7OoAi1sptT87ne/49///jclJSXExcWxd+9eu0sKnL1fwJd3mbYzzgwOdihvioiErNhUOGMJfHIDbHzIXFt5I9TvhPG/6zR0oGh1OYVL1lLuqW+7luVOYP7sPArGabXiQAmbd87GxkYuvPBC5syZY3cpgbfyFrCaTXvsHeDOs7ceERE5PGcMnLwQxv23/9raBbDix+BrbrtUtLqcOYtWdgg0ABWeeuYsWknR6vJgVRzxwibUFBYWctNNN3HCCSfYXUpglb8OFW+YdvJwbVYpIhJOHA44sdBM7Gjd7/vrR81Cfd4GvD6LwiVr6epBU+u1wiVr9SgqQMIm1ByNhoYGampqOhwhxeeFVbf5z8ffCa54++oREZGjM+p6OPVpcMaa8x0vwruzKd5Y2qmHpj0LKPfUU7x5d3DqjHARHWoWLFiA2+1uO3JyQmwRuy2LYO9npp02GYZeZG89IiJy9IZeBGe8DK4kc17xBpUf/bpHn1pZ233wkZ6zNdTMnTsXh8NxyGPdunVH/fXnzZuHx+NpO7Zv3x7A6o9R8wH4/Jf+85P+qMHBIiLhLuscOPMNiHUDkFFf3KNPy0hJ6M2qooats59uueUWrrzyykPek5ube9RfPz4+nvj4EH2c89VfYP8O087+FgycaW89IiISGOnTYNY78Pa55FtryIqtoqJpABadF1N1AJnuBPKHpwW9zEhka6hJT08nPT3dzhLs0VBtdnoF0zsz4ff21iMiIoHVbwLMeh/X0jOZn/0wc7begQOrQ7Bpbc2fnYfLqdXjAyFsnnds27aNkpIStm3bhtfrpaSkhJKSEvbt22d3aUdu3b3Q5DHt3Kug7zhbyxERkV6QOgpmvUNB5jYWDr2TzNhdHV7OdCew8LKJWqcmgByWZYXFPLIrr7ySJ554otP1t99+mxkzZvToa9TU1OB2u/F4PKSmpga4wh5q2geLh0DjHnDEwHlfaxduEZFIVrMBls7Au7+C4rqxVMZNIOOUX5M/arh6aHqop+/fYRNqAiEkQs26e2HlTaY9/AqY+rg9dYiISPC0BBsOlJnzvuPhrKUQ39/WssJFT9+/w+bxU0TwNcG6P/vPx9xqXy0iIhI8qcfBWW9DYsujpr2fwVuzzBhLCRiFmmDa+jTsb5lWnv1t6DvW3npERCR4UkfBWe/4g82eEnjrbGjQwnuBolATLJYFa//oP8+7rft7RUQkMqWOMj02CZnmfM8qeOeb0Fxnb10RQqEmWMpeAc9q0x4wFdKn21uPiIjYI3V0x2BTvQLe/54ZoiDHRKEmWDqMpbmt07b0IiISRdzHw8zX2lYeprwIProKLJ+9dYU5hZpg2Pc17HzLtPuMhMHn2VuPiIjYr9+JcMYScLVskbDlKVh5ixmuIEdFoSYYvm63vs6IH2mPJxERMTJOg1OfAYfLnK+/F9b+wdaSwpneXXubzwtfP2baDicM/6G99YiISGgZfB7kP+I//2webPpbl7d6fRbLN1WzuKSU5Zuq8frUq9OerXs/RYWdb/mncWcVQNIge+sREZHQM+IqaKiCktvNefFPIK4/5FzQdkvR6nIKl6yl3FPfdi3LncD82XnaaqGFemp629eP+tu5P7KvDhERCW1jboXjbzZtywcf/gAqPwBMoJmzaGWHQANQ4alnzqKVFK0uD3a1IUmhpjc17IbtL5p2/AAYNNveekREJHQ5HHDSXTDscnPua4D3v4O35msKl6ylqwdNrdcKl6zVoygUanrX1n+Yf5QAwy4DV5y99YiISGhzOOGUv0HmOea8YRfFS27s1EPTngWUe+op3qyViRVqetPWZ/zt3Kvsq0NERMKHMxamP2MW6QMq93p69GmVtd0Hn2ihUNNb6nfBrg9NO3W0WY9ARESkJ+L6whkvQ1w/MmL39OhTMlISeremMKBQ01vK/u1fGXLQ+fbWIiIi4SdlJJz2Avl91pMVW4WDrlcbdmBmQeUPTwtufSFIoaa3lP7L39YKwiIicjQGzsSVfx/zsx8GHJ2CTeuGO/Nn5+FyavsdhZre4K2H8tdMOz4d+p9ibz0iIhK+Rv6Egvx8Fg69k8zY6g4vZboTWHjZRK1T00KL7/WGirf828gP+jY4XfbWIyIi4e2kuymomc3ZqVdTXDeWytgTyDj1D+SPHKQemnYUanpD+0dPg/ToSUREjpEzBk59GtfrU5nq+AL4Ana54LgnDvup0USPnwLNsqB0iWm7EiDrbHvrERGRyBDnhtNfhJg+5nzLk7Dpr/bWFGIUagKt9is4UGbaGTMhJtneekREJHKkjoYp7YLMJz+F3avsqyfEKNQEWtWH/nbG6fbVISIikWno9+G4603b1wAfXAiNPVugL9Ip1ARa1Qf+dvp0++oQEZHINfFPkHayae/bBCt+ZIY/RDmFmkBr7alxxkH/yfbWIiIikckVD9Ofhdi+5nz7P2H9/9paUihQqAmk+kozpgYgbZIZKCwiItIb+gyDqf/Pf77qVqhabls5oUChJpCqlvnbevQkIiK9bfBsGHObaVvN8OEPonp8jUJNIHUYT3OqfXWIiEj0GP9b/y/S+7fByhttLcdOCjWBtKtdT82AafbVISIi0cMZC1OfhJgUc/7147D9JTsrso1CTaBYPtj7uWn3yYWEdHvrERGR6NFnGExqN1C4+CdmnGeUUagJlP3b/fs9ucfaW4uIiESf3Cth8Pmm3VBlgk2UTfNWqAkUz5f+duoY++oQEZHo5HBA/sMQ3/KkYMdi2Bxde0Mp1ARKTbtQ41aoERERGyRkwJRH/Oef/AzqttpXT5Ap1ARKh56aPPvqEBGR6Db4fPMoCqC5FpZfacZ9RgGFmkCpWetvu4+3rw4REZGJ90LSENOufAc2LLSzmqBRqAkEy/L31CQOgthUe+sREZHoFueGqe3G03x2B+wvs6+eIFGoCYSGKmjcbdoaTyMiIqFg4AwYcbVpN9XAyptsLScYFGoCof0grJTj7KtDRESkvQl/gPgBpr3tWSgrsreeXhYWoWbLli1cffXVDB8+nMTEREaMGMH8+fNpbGy0uzTjQIW/nZBpXx0iIiLtxfeHk+72n3/8n9C83756ellYhJp169bh8/l46KGHWLNmDffccw8PPvggd9xxh92lGfU7/e1EhRoREQkhw38IGTNMu24zrP6treX0phi7C+iJgoICCgoK2s5zc3NZv349Cxcu5O677z7EZwZJffuemoH21SEiInIwhwNOXgivngi+JvjyLhh2KfSNvNXvw6Knpisej4e0tLRD3tPQ0EBNTU2Ho1e076nR4ycREQk17uMhb65pW83w8XURuXZNWIaajRs3ct9993Httdce8r4FCxbgdrvbjpycnN4p6IB6akREJMTlzYM+I0y76gPYvMjeenqBraFm7ty5OByOQx7r1q3r8DmlpaUUFBRw4YUXcs011xzy68+bNw+Px9N2bN++vXf+Ih16ahRqREQkBMUkwsn/5z///BcRN2jYYVn2beFZVVVFdXX1Ie/Jzc0lLi4OgLKyMmbMmMEpp5zC448/jtN5ZJmspqYGt9uNx+MhNTWAC+QtGQW1G8yiexd6Avd1RUREAu2db0PZv017/O9gbIhMujmEnr5/2zpQOD09nfT09B7dW1paysyZM5k0aRKPPfbYEQeaXtXQEszie/Z3ERERsc1Jf4TyV82YmjW/hxE/NhthRoAQSgbdKy0tZcaMGQwZMoS7776bqqoqKioqqKioOPwnB4O3pfsuJsneOkRERA7HnQcjWoZvNNfCF4X21hNAYTGl+4033mDjxo1s3LiRwYMHd3jNxqdnrQWAt960XQo1IiISBk74NWxZBM11sPEhGP0zSB1td1XHLCx6aq688kosy+rysF1roAFwJdpXh4iISE8lZsKY203b8kLJ7fbWEyBhEWpCmq99qEmwrw4REZEjMeZmSMwy7R2LofI9e+sJAIWaY+Vr9redsfbVISIiciRikuHEdlsmrLrVDKkIYwo1x8pqF2ocLvvqEBEROVLDr4C+J5h2dTGUh/cu3go1x8ry+tsKNSIiEk6cLhg333+++n/CurdGoeZYtd87w6Fvp4iIhJmc74C7ZXPLXcth59v21nMM9C58rJzx/ra3wb46REREjobDCWN/4T9f/T/21XKMFGqOVUyyv91cZ18dIiIiR2vIRZAyyrQr34HKD2wt52gp1Byr9mvTKNSIiEg4cro67gEVpr01CjXHyunyr0/jjazdTkVEJIoMuwSSh5t2xeuwq9jeeo6CQk0gtG6PoJ4aEREJV85YGDvXfx6GvTUKNYHQOq5GoUZERMLZ8CsgqWWPxbKXoWa9vfUcIYWaQGgNNXr8JCIi4cwVD6P/y3/+1QP21XIUFGoCof3jpzBetEhERITcH/knwXz9ODTV2lrOkVCoCYS4fuaj5YUmj721iIiIHILXZ7F8UzWLS0pZvqkar++gX8bj02DYpabdXAubnwx+kUcpxu4CIkLSIH97fynE9bWtFBERke4UrS6ncMlayj31bdey3AnMn51Hwbgs/42jboBNfzXtr+6H4+aAwxHkao+cemoCIbFdqDlQal8dIiIi3ShaXc6cRSs7BBqACk89cxatpGh1uf9iv/GQfppp13wJO98KYqVHT6EmEFpHigPs32FfHSIiIl3w+iwKl6ylq1GfrdcKl6zt+Chq1A3+9lf392Z5AaNQEwgHP34SEREJIcWbd3fqoWnPAso99RRv3u2/mPMdSMw27dJ/Qd3W3i0yABRqAqF9T40eP4mISIiprO0+0HR7nzMWRl5n2pYPNj7SC5UFlkJNICSqp0ZEREJXRkrC0d038sdmF2+ALU+ZcBPCFGoCISEDHC0TyQ5oTI2IiISW/OFpZLkT6G7+kgMzCyp/eFrHFxKzIPNs067bAlXLerHKY6dQEwgOp39cTd1WLcAnIiIhxeV0MH92HkCnYNN6Pn92Hi5nF7Fn2GX+9pZFvVJfoCjUBErq8eZj4x6or7C3FhERkYMUjMti4WUTyXR3fMSU6U5g4WUTO65T097gC/wr5297FrwNvVvoMdDie4HiHgvlr5m2Z43pshMREQkhBeOyODsvk+LNu6msrScjxTxy6rKHplVsHzMTastT5hf3slch54Kg1Xwk1FMTKO5x/vbe1fbVISIicggup4OpI/pz/oRBTB3R/9CBplWYPIJSqAmUvu1CjUehRkREIkjmLEgYaNqlS6Bxr63ldEehJlBSx/jbe9fYV4eIiEigOWNg6MWm7WuEbc/bW083FGoCJbYPJA83bc9qzYASEZHIMuwSf3vHYvvqOASFmkByjzUfm/fB/m321iIiIhJIaZP92ybsfBOa99tbTxcUagKp/biavV/YV4eIiEigORyQ/S3T9tbDzrftracLCjWB1O8kf3vXR/bVISIi0hsGfdvfLn3Zvjq6oVATSOmn+ttVH9hXh4iISG/IPAuc8aZd9nLIjR9VqAmkpEH+wcLVK8DbaG89IiIigRSTDAPPNO39O2Dv5/bWcxCFmkBr7a3x1sOelfbWIiIiEmgh/AhKoSbQ0qf723oEJSIikWbQt/ztcA81V1xxBe+9915v1BIZFGpERCSSJQ/1bw1UvSKkVhc+4lDj8XiYNWsWxx13HHfeeSelpaW9UVf4co+BuH6mXfVByA2iEhEROWaZZ7U0LNi13NZS2jviUPPSSy9RWlrKnDlzeOaZZxg2bBjf+MY3eP7552lqauqNGgE477zzGDJkCAkJCWRlZXH55ZdTVlbWa3/eUXM4YUDLuJqGaqhZZ289IiIigRais32PakxNeno6N998M5999hkrVqxg5MiRXH755WRnZ3PTTTexYcOGQNfJzJkzefbZZ1m/fj0vvPACmzZt4nvf+17A/5yAyGj3CKpiqX11iIiI9IYBERRqWpWXl/PGG2/wxhtv4HK5+OY3v8kXX3xBXl4e99xzT6BqBOCmm27ilFNOYejQoUybNo25c+fy0UcfHbJ3qKGhgZqamg5HUGR9w98uC61BVCIiIscsKRv65Jp2dTF4G+ytp8URh5qmpiZeeOEFvv3tbzN06FCee+45brzxRsrKynjiiSd48803efbZZ/nNb37TG/UCsHv3bp566immTZtGbGxst/ctWLAAt9vdduTk5PRaTR30PQGSWv6snW9D077g/LkiIiLB0joxxlsPu0NjCZMjDjVZWVlcc801DB06lOLiYj755BOuu+46UlNT2+6ZOXMmffv2DWSdANx+++0kJyfTv39/tm3bxuLFh94ldN68eXg8nrZj+/btAa+pSw6Hfx6/rxEq3gzOnysiIhIsITjb94hDzT333ENZWRkPPPAAEyZM6PKevn37snnz5sN+rblz5+JwOA55rFvnH2h76623smrVKl5//XVcLhc//OEPsQ4xuyg+Pp7U1NQOR9Bkt1ucSI+gREQk0rQPNbs+tK+OdhzWoVJBL6uqqqK6uvqQ9+Tm5hIXF9fp+o4dO8jJyWHZsmVMnTq1R39eTU0Nbrcbj8fT+wGn+QC80B+8ByAhE75TamZGiYiIRALLBy+kQ+NuiO8P/1FlnlT0gp6+f8f0yp/eQ+np6aSnpx/V5/p8PsAMBg5JMYmQOQtKl0B9BexZBWmT7K5KREQkMBxO6J8P5UVmCZMDZWYPRBuFRdfBihUruP/++ykpKWHr1q289dZbXHzxxYwYMaLHvTS2COH9MURERI5Z33H+9t7V9tXRIixCTVJSEv/85z8566yzGD16NFdffTUnnngi7777LvHx8XaX173sdvtjbH/RvjpERER6g7tdqPGssa+OFrY+fuqpE044gbfeesvuMo5c0iDTNVddDHs/Mym2faoVEREJZ+6x/rZHPTWRb9hl/vaWp+yrQ0REJNDcY4CWwcF6/BQFhn4fHC7T3vKUGS0uIiISCWKS/SsL16y1/T1Ooaa3JWRAVoFp798Ole/bW4+IiEggtT6Caq6Duq22lqJQEwwdHkEtsq8OERGRQAuhGVAKNcEw+DyI6WPa254z+2SIiIiEIa/PYvmmahaXlLJ8UzXePnn+F2vX21cYYTL7KezFJEHOd2HzE9DkMWvWDPme3VWJiIgckaLV5RQuWUu5x//LeVaffszvP5UC93LYX2pjdeqpCZ7h7R5Bbf5/9tUhIiJyFIpWlzNn0coOgQagYp/FnK13UOSZCgcUaqJDxkxIbFk+uuzftg+mEhER6Smvz6JwyVq62izSavnfwrKf4K0rC25hB1GoCRanC0Zea9qWDzY8aG89IiIiPVS8eXenHpr2LJyUN6VTXJEQxKo6U6gJppHXgDPWtDc9ogHDIiISFipre/Z+VVnns3WtGoWaYErMhCEXmXZDNWx92t56REREeiAjpWc9MBkxu6C+qper6Z5CTbCNusHfXn8fWF09oRQREQkd+cPTyHIntG6I0IkDi6zYKvKT19g6WFihJtj6T4G0yaa9ZyXs+sjeekRERA7D5XQwf7ZZj+bgYNN6Pj/7YVwOn63TuhVqgs3h6Nhb89X99tUiIiLSQwXjslh42UQy3R0fRWW6E1g4c6dZpwagqcaG6gwtvmeHod+HVT+Hhl2w/Tk48Ccz3kZERCSEFYzL4uy8TIo376aytp6MlATyh6fh2vwoVLfc5K2zrT711NjBlQAjrjFtXxOsu8feekRERHrI5XQwdUR/zp8wiKkj+uNyOsCV7L+hWaEm+oy6AZzxpr3hATMbSkREJBzFtA81+20rQ6HGLknZMOJq026ug3X32lqOiIjIUYtJ8rfVUxOl8m73L8b31V+gca+t5YiIiByV9j01XvXURKfkITD8CtNuqtFMKBERCU8xGlMjAHlzweEy7XX3QFOtvfWIiIgcKYUaASBlBAy9xLQbd8OGhfbWIyIicqSccf62r9G+Mmz7k8Vv7B20rcm47k+2jhwXERE5YpbX33bYtwSeQk0ocB/v3+iyvhLW/8XeekRERI6Er9nfbh1SYQOFmlAx7lfgaPm/Y+0CW3c5FREROSIdemoUaqTvWMhtWbemqQZW/8beekRERHrK1+RvO/X4SQBOLPSPIN/wINR8ZW89IiIiPeE94G+7Em0rQ6EmlCRmwZhbTdtqhpK59tYjIiLSE956f1uhRtqM+bkJNwA7XoTKD+ytR0RE5HDUUyNdikmGE9qNp1n1c7As++oRERE5nPYL7rXfByrIFGpCUe5V4B5n2tUrYNtz9tYjIiJyKPU7/e34dNvKUKgJRU4XnHSX/7zkNluXnRYRETmk9qEmMdO2MhRqQlXWuZB5jmnXbYUvNMVbRERCVH2Fv50w0LYyFGpClcMBk+8HZ7w5X/dn2PuFvTWJiIh0pX1PTYJ6aqQrqce17AuFmeJdfC1YPntrEhEROdiB9j01GbaVoVAT6vJuh5RRpr1rOWz6m731iIiIHMR7YCfL953A4tpvsHzLPrw+e2bthl2oaWhoYMKECTgcDkpKSuwup/e54iH/Qf95ye1m00sREZEQULS6nOmf/IqLv17Af22+nosf+Yjpf3iLotXlQa8l7ELNbbfdRnZ2tt1lBNfAmTDsctNu3AMrb7G3HhEREUygmbNoJeVN/Ttcr/DUM2fRyqAHm7AKNa+++iqvv/46d999t92lBN/EuyGun2lvWQQVS+2tR0REoprXZ1G4ZC3mQZOjw2utD58Kl6wN6qOosAk1O3fu5JprruHJJ58kKalnqxU2NDRQU1PT4QhbCRkw4Q/+8+JroWmfffWIiEhUK968m3JPfbevW0C5p57izbuDVlNYhBrLsrjyyiu57rrrmDx5co8/b8GCBbjd7rYjJyenF6sMghFXQ/qppr1vE6y61d56REQkalXWdh9ojua+QLA11MydOxeHw3HIY926ddx3333U1tYyb968I/r68+bNw+PxtB3bt2/vpb9JkDiccMrj4Grpqdr4IJS9amtJIiISnTJSEgJ6XyA4LMu+3RKrqqqorq4+5D25ublcdNFFLFmyBIfD/8zO6/Xicrm49NJLeeKJJ3r059XU1OB2u/F4PKSmph5T7bba8CB8PMe0E7Pgm19AfP9Df46IiEgAeX0W0//wFhWeA1gHjakBM8om053AB7eficvZ+fUj0dP3b1tDTU9t27atw3iYsrIyzj33XJ5//nmmTJnC4MGDe/R1IibUWBa8800oLzLnQ74P05+2tyYREYk6RZ9vZ87fPwMsrHYPf1ojzMLLJlIwLuuY/5yevn+HxZiaIUOGMG7cuLZj1CizGN2IESN6HGgiisMBU/7mnw217RnY8g97axIRkahTMMTDwqF3khnb8alLpjshYIHmSMQE9U+TwEnKhpMXwoc/MOcf/ydknA5Jg+ytS0REoodnDQXu5ZyduoLigXdTOeAiMlISyB+edsyPnI5GWIaaYcOGEQZPzXrf0O/DjsWw9R/QtBc++hHMLDI9OSIiIr1t13IAXA4fU8ccB4Ps/cU6LB4/ySFMvh8SW1ZYrngd1t1jbz0iIhI9qj70twdMs6+OFgo14S4+DU55zH9ecjtULbOvHhERiQ5N+2DPKtN2jzXvRzZTqIkEWedAXssaPlYzfHAR1FfZW5OIiES26mKwvKbdujCszRRqIsWJv4GMM0z7QCksvxwsn701iYhI5Gr/6Cl9un11tKNQEymcMXDqPyBhoDkvfw3W3GlvTSIiErmqPvC31VMjAZeYZYKNo+X/1i/mQ8Vb9tYkIiKRx+dtm/lEYhYkD7e3nhYKNZFm4Ew4odC0LR8suxj2l9lbk4iIRJY9q6C51rQHnBoyS4ko1ESisXdAVoFp11eaYONrtrcmERGJHKVL/O3Ms+yr4yAKNZHI4YSpT0JSyxYSle/BypvsrUlERCJH6b/87UHn2VfHQRRqIlXCADj1WXDGmvOv7oev/s/emkREJPzVbYM9JaaddrLZtidEKNREsvSpkP+w//zTn0HFm/bVIyIi4W9Hu16awaHTSwMKNZEv90oYc6tpW154/0Ko+crWkkREJIyF6KMnUKiJDuMXwKDZpt20F979NjTusbUkEREJQ40eqHzHtJOHQt8TbC3nYAo10cDpgmlP+f/x1W4wPTa+JnvrEhGR8FJe5H/vGHR+yEzlbqVQEy1iU+CMJRCfbs53LoVPfgaWZW9dIiISPrY+7W+H2HgaUKiJLslD4fQXwRlnzjc+COv/Ym9NIiISHg7shNKXTTsxy7/fYAhRqIk26adC/iP+85U3wZanu79fREQEYMsisFoWch1+hdlzMMQo1ESj3B/C2F+2nFjw0Q+h/HVbSxIRkRBmWfD1o/7z3Kvsq+UQFGqi1Ym/gRHXmLavCd7/D9hVbG9NIiISmqqLwbPWtNOnQ+ooe+vphkJNtHI44OT/g8HfMefNdfDuN8Gzzt66REQk9HTopfmRfXUchkJNNHPGwKl/9w/2aqiGt8+F/TvsrUtEREJH837Y8g/TjkmGIRfaW88hKNREO1cCnL4Y+k0w5/u3wdsF0LDb1rJERCREbH0GmmtNe8j3IbaPvfUcgkKNQJwbZrwKfXLNuWeNWXW4eb+9dYmIiL0sH6z7k/98xI/tq6UHFGrESMyEma9DwkBzvms5vHc+NB+wty4REbFP2SvmF12AAdPMRskhTKFG/FJGwMwiiE015xVvKtiIiESztX/0t/NuA8Drs1i+qZrFJaUs31SN1xc6K9OH3so5Yq9+E2BGkRkw3FwLFW/AexfA6S9BTKLNxYmISNBULYeq90079XgYNJui1eUULllLuae+7bYsdwLzZ+dRMC7LpkL91FMjnaVPNT02MS2DwSpeh/e/A976Q3+eiIhEji/b9dKMuZWiNTuZs2hlh0ADUOGpZ86ilRStLg9ygZ0p1EjX0qfBzNf8wab8NXhPwUZEJCp41sGOxaadmI13yCUULllLVw+aWq8VLllr+6MohRrpXvq0jj025UXw3n8o2IiIRLov76Itroy+keJtdZ16aNqzgHJPPcWb7V0ORKFGDi39VDPdOybZnJe/Cu9/F7wNHW4L5YFjIiJyBDxfwuYnTDs2FY67lsranv0y29P7eosGCsvhZUw3weadb5jtFMpeMY+iTnseYpJCfuCYiIgcgc/mgeU17TG3QWwqGSlNPfrUjJSEXizs8NRTIz2TcRrMeKVjj83bBRSVbAr5gWMiItJDle93GEvD8TcBkD88jSx3Ao5uPs2B+WU2f3haUMrsjkKN9FzG6S2PolIA8FZ+SOE/V4T8wDEREekBy4JVt/rPT/wNxCQB4HI6mD87D6BTsGk9nz87D5ezu9gTHAo1cmQyToNZb0N8f4rrxlLe6O721lAZOCYiIj2w/XmoXmHa7rEw/MoOLxeMy2LhZRPJdHd8xJTpTmDhZRNDYriBxtTIkUubBLPeo/LZX/TodrsHjomIyGF4G6Fknv98wh/B6ep0W8G4LM7Oy6R4824qa+vJSDGPnOzuoWmlUCNHx51HxtTfwqYth73V7oFjIiJyGBsfhH2bTHvgTMj+Rre3upwOpo7oH6TCjoweP8lRy8/LIys1FkeXo2pCZ+CYiIgcwv5S+OyX/vOT7gJHaPS8HKmwCTXDhg3D4XB0OH7/+9/bXVZUczkdzD/vBMDRKdi0nofCwDERETmET24we/0B5P7IDDEIU2ETagB+85vfUF5e3nb89Kc/tbukqNc2cCz1oIFjsbtYeNpGCsZm2lSZiIgc1vZ/wo6XTDshw/TShLGwGlOTkpJCZmbP3yQbGhpoaPCvfFtTU9MbZUW9toFjX1dR+fnDZOx6lvzkNbg8Pli2Ak55FFwaVyMiElIa95pemlaT7oP48B4uEFY9Nb///e/p378/J510EnfddRfNzc2HvH/BggW43e62IycnJ0iVRh+X08HUkRmc/51fMPW0K3E5Wh5Hbf0HLD0T6ivtLVBERDoquR0OtCyQmv1tGHKhvfUEgMOyrLBYGe3Pf/4zEydOJC0tjWXLljFv3jyuuuoq/vznP3f7OV311OTk5ODxeEhNTQ1G2dFrx7/gw4vBu9+cJw+DGf8Gd56tZYmICFD5Hrx5hmnH9IFvrYXk0P3Fv6amBrfbfdj3b1tDzdy5c/nDH/5wyHu+/PJLjj/++E7XH330Ua699lr27dtHfHx8j/68nn5TJEB2r4R3Z8OBMnMemwrTn4ess+2tS0QkmjXvh6KJULPenE/6C4wO7TGqYRFqqqqqqK6uPuQ9ubm5xMXFdbq+Zs0axo0bx7p16xg9enSP/jyFGhvsLzXBZs8qc+5wweQH4Lhr7a1LRCRarfgJbHrEtPtPgbM/7HKhvVDS0/dvWwcKp6enk56eflSfW1JSgtPpJCMjI8BVSUAlDYJZ78GyS6H0X2bn14+vA89amHg3OGPtrlBEJHpsfcYfaFxJcMrjIR9ojkRYzH5avnw5K1asYObMmaSkpLB8+XJuuukmLrvsMvr162d3eXI4sX3gtH+aQWnr/mSuffUX03sz/VlI1LRvEZFet28zFP/Efz75fnB3Ht4RzsJi9lN8fDxPP/00Z5xxBmPHjuV3v/sdN910Ew8//LDdpUlPOV2mZyb/YX/vTNX75rlu1Yf21iYiEul8TfDhD6CpZWmToZdA7pW2ltQbwmb2UyBoTE2I2LUC3v8uHCg1544YmPhnvCOvp3jLnpDcJE1EJFx5fRbFbyygcuNrZMTuIX/gAVzf/NRM3ggTYTGmRqLUgCnwjZXwwfeh8h2wmil66x8U/iOT8vqkttuy3AnMn50XEtvZi4iEo6LV5RS+tJLyfeOB8QBkVTqZP6KOgnHhE2p6KiweP0kESsiAM9+AMbdS5JnKnK13UF6f2OGWCk89cxatpGh1uU1FioiEr6LV5cxZtJLyfR0fyFTs80Xsz1aFGrGPMwbv+D9QuGtey/aXHR81tf5nWLhkLV5f1DwlFRE5Zl6fReG/VmNhEU0/WxVqxFbFm3dTXufk4P/oWllAuaee4s27g1qXiEg4K/66ivKaRqLtZ6tCjdiqsrY+oPeJiAhUfvFoz+6LsJ+tCjViq4yUnu3enVH7NkTPRD0RkaO38WEyqv7eo1t7+jM4XCjUiK3yh6eR5U7opoMUHPjIiq0if9sVZoPMxr3BLE9EJLzsfBs+vp785DVkxVbhoOtfBh2YGab5w9OCW18vU6gRW7mcDubPNjt3HxxsHC3/Oz/7YVwOH2x7Bl6dAFXLglukiEg48Hxp1gCzmnE5fMyfXA44uvnZCvNn50XcWmAKNWK7gnFZLLxsIpnujt2gme4EFl42iYJv3QSxbnOxbiu8eTp88T/g89pQrYhICKrdBG+dBY17zHn2NymYfdshfrZOjMg1wLSisIQMr8+iePPurlcUrttqNsVsv6VC+nQ45TFIGWlPwSIioWD/Dnhjuvk5CdBvIsx6u23F4EP+bA0TPX3/VqiR8OFrhtW/hTX/A5bPXHMlwYTfw6jrwaGORxGJMvWVpve6Zr05d+fBWe9CwgB76wqwnr5/611AwoczBk78NZz1DiQPM9e8++HTn8HSmab7tR2vz2L5pmoWl5SyfFN1xC0yJSJRrmE3vHW2P9D0GQFnvhlxgeZIqKdGwlPTPii5HTb8n/9au16bojU7KVyylnKPfw0G7SUlIhGjqRbemgXVxeY8KQfOfh+Sh9pbVy/R46cuKNREoIq3YMXVULel7VIRVzPn8+90msjY+gQ5UgfIiUiUaNoH734bKt815wkDYdZ7kDrK3rp6kR4/SXTIPBO++QUc958AeC0nhV9Ob9nvpKNI3u9ERKJEQ7WZ5dQaaOLSzObAERxojoRCjYS/2D5w8gNw5lKKvWdR3pROtO13IiJRYH8ZvHmG/5FTXD+Y+Rr0PcHeukJIjN0FiARM5plUjh0Ga7887K2Rtt+JiES42k1mDE3ro/aETDjzdQWagyjUSETJ6Ovu2X0Rtt+JiESwPZ/D2+dCfYU5Tx4OZ70JfXLtrSsEKdRIRGndS6rCU9/ljicOfGQm7Cc/Yx/QP9jliYgcUqeF8vqsw/Xet6Fpr7nBPc48ckrKtrXOUKVQIxGldS+pOYtW4oAOwcaBD3Awf+D/4nrlR3DCr+H4G8EZa0utIiLtFa0u77wURewu5mePocC9HPpPgRmvQHxkbUIZSBooLBGn272kkmHhyP8zPxy8+6HkNnh1IlR+YFOlIiJG0epy5ixa2SHQAFQ0pTFn6x0UOa41C+sp0ByS1qmRiNXlfifNHvjsly2L9rX7pz/0EpiwAJKH2FaviEQnr89i+h/e6hRoWjmwyHQn8MHtZ4Xdnk2BonVqJOq5nA6mjujP+RMGMXVEf/PDIK4vnHw/nFsMaZP8N2/9O7w8Gj77lVnYSkQkSIo37+420ABYOCj3NGgpih5QqJHo1H8ynLMCJj8A8S0Dhr31sOa3sOQ42PQo+Lz21igiUaGnS0xoKYrDU6iR6OV0waj/hNkb4Pib/QOG6yvM1guvTYadb9tbo4hEvIy6no3r01IUh6dQIxLXDyb+Cb61FgZ/x399TwksPRPeuwBqNthVnYhEqub9sOIn5G+5hKzYqpYZmp05MBvy5g/XIOHDUagRaZUyEk7/J5z1DvSb6L++YzH8Ow8+vRHqK20qTkQiyt7V8NrJsOkRXA4f87MfBhydNnhpPZ8/Oy9qBwkfCYUakYMNPAMKPoZTHofElgWurGZY/7/wr1wouQMa99haooiEKcuCDQ+ZQONZa665kig458csvHRS56Uo3AksvGwiBeOybCg2/GhKt8ihNNfB2rvgyz+C94D/eqwbjr/FLN4Xm2JbeSISRhr3wIprYPsL/mt9T4RTnwb3GKCbpSjUQ9Pj92+FGpGeOFABa+6EjQ+Br7HtsjcuneL+hVSmFpDRN1U/gESka1XL4MOLYf82/7VRN8BJd4FLA4APR6GmCwo1cszqtsHq/4GvH6Nobz6FZT+hvCm97eWs1HjmnzdWXcUiYjQfgNW/gS/vAqtlmYi4NDjlURh8vr21hRGFmi4o1EigFBWvZM4/y1rWJPb3zDhariy8ZDwFJ+bYUpuIhIjK92HFj6H2K/+19NNg2lOQrJ8PR0IrCov0Eq/PonDpHiwccNBcBatlG83C59/Bu/ZeMyZHRKJLUw18/J/w5un+QOOMgxN/C2e9pUDTixRqRI7Q4Zc0d1LemEbxB4/C4qHwxW+gQcubi0SF0lfg32Nhw0L/tQFT4RslMO4X4IyxrbRooFAjcoR6vKR5Uz9oqIYv5ptws+pWOFDey9WJSG/w+iyWb6pmcUkpyzdV4/UdNHKjfhcsuwze/Rbs32GuxSTDpL/ArPfbZjdJ71JkFDlCPV2qPCNnMtR9aAYHNu+DL++G9X+B3CthzG2QMqJ3CxWRgChaXU7hkrUdemiz3AnMn51HQV46bHwYPv8VNLbrkc08B/Ifgj7Dgl9wFAurnpp///vfTJkyhcTERPr168cFF1xgd0kShfKHp5HlTui08mertiXNv3232VfquDngjDcv+hrND8CXR8EHP4Cq5WYxLhEJSUWry5mzaGWnR84VnnrmLFpJ0aIfwifX+wNNXD+zcOfMIgUaG4RNqHnhhRe4/PLLueqqq/jss8/48MMPueSSS+wuS6KQy+lg/uw84OBhwl0sad5nOJz8f3D+Fsi7HWJaFuqzfLDtGXhjGryWD5sXgbchWH8FEekBr8+icMlauvq1w1zzUbjhHLxWy1vpsMvMHnK5V4BD61XZISymdDc3NzNs2DAKCwu5+uqrj/rraEq3BNIhu6S7W6emcS989YDZcqGhquNrCQNNr87IayExs/cKF5EeWb6pmosf+eiw9/1j/KNMnXULpE8LQlXRqafv32ExpmblypWUlpbidDo56aSTqKioYMKECdx1112MGzeu289raGigocH/229NTU0wypUoUTAui7PzMo9sSfO4vmYGxJhbYOvTJtzsKTGv1e+EL34Na34HQ34Ao38G/ScH4W8iIl3p8aSA4++FdE3TDgVh8fjp66+/BuDXv/41v/zlL3n55Zfp168fM2bMYPfu7qfKLliwALfb3Xbk5OgfnQSWy+lg6oj+nD9hEFNH9O/5FgmuBDNguGClmRmR8z1wuMxrvibY8qTZ8O71U2HL03o0JWKDjOSevUVmpCb1ciXSU7aGmrlz5+JwOA55rFu3Dp/PB8AvfvELvvvd7zJp0iQee+wxHA4Hzz33XLdff968eXg8nrZj+/btwfqrifSMwwEZ0+G05+C8ryFvrllCvdWuZbDsYnhpMHx6M+xdY1+tItHC2wgbHiJ/9SlkxVbhwNflbW2TAoandfm6BJ+tj59uueUWrrzyykPek5ubS3m5WdsjLy+v7Xp8fDy5ubls27atu08lPj6e+Pj4gNQq0uuSh8CEBTDuV7Dl7+bRlGe1ea1hF6y/xxwDpsKIH8OQiyC2j701i0QSXxN8/bh5BFy3FRcwP/th5my9AwdWy4rhRqdJARISbA016enppKenH/a+SZMmER8fz/r165k+fToATU1NbNmyhaFDh/Z2mSLBFZMEI38MI66Gyndh4yOw/QXwtTyC2rXcHJ/+Fwy92ASc/idrtoXI0fI2mke+q38LdVs6vFQwpj8LT+pH4dv1HSYFZB5uUoDYIixmPwHceOONPP/88zz66KMMHTqUu+66iyVLlrBu3Tr69evXo6+h2U8Sthp2w5anYNMjsPeLzq+7x5lwM+wSSDj8LwoiglkFeOND8NX9UF/R8bWsb8AJv4YB+YCZ3n1EkwIkoCJul+6mpibmzZvHk08+yYEDB5gyZQr33nsvY8eO7fHXUKiRsGdZsPsT2PRX84iqeV/H1x0us5LpsEtg8PkQm2JPnSKhzPMlrL8XNv8/8B40wynr3JYwc4odlUk3Ii7UBIJCjUSUpn2w7TkTcHYt6/y6KxEGzTYBJ6sAXBpfJlHMsqDiTVj3Zygv6viawwmDvwPH3wLpU+2pTw5JoaYLCjUSsTxr4esnzNo3+83gea/lpLhuLJVN/chIbCI/Lw/XsEsg4wxwumwuWCRImmpb1oT6i3/gfauYFDN2bfTPzOrfErIUarqgUCMRz/JB1TKKlr1F4afDKW/yTzXNiq1ifvbDFGRugSHfh5z/gAHTFHAk8lgWVH0IX/8Ntj4L3v0dX08eCqP/ywSaWL0XhAOFmi4o1Eg0aN2A7+D/sM1aGw4WDr2TAvdyczF+gHlENfgCyDwbYhKDXK1IAB0oN+NkNj0KtV91fn3ANDj+JvPv3RkWC+pLC4WaLijUSKTz+iym/+GtTjsKt3JgkRm7iw+OvxqX46AFxVxJZpDk4Atg0Lcgvn/vFyxyFDrMROoTQ37cR7g2/w3KXgHL2/Hm2FQYeonpldG2I2ErovZ+EpGeKd68u9tAA2DhoLwpneJhTzHV97wZMNlcZ1707ocdL5rD4YL000zAGXyexhtIyChaXU7hv9ZQXuPfOiQrdg/zs3dR4G4XaAbOhNwfmcesMdrGIFoo1IhEkB5vwJd8Gkz4gZnOWrEUdrwEOxb7dw63vFD5jjlW3ggpx5mp4lnnwMAZGocgwedrouiDN5nzirfl0ap/jZiKpv7M2XoHC0c9RMGkCTDiKuiTa1OhYieFGpEIkpGScGT3uRLMo6ZB34KTH4Tqj0zA2f4i7Nvk/4TaDebY8AA4YsxWDVnnmKCTNkmDjaV3eBth51LY9jze7Ysp/PxuLAbQPtAAWDhxAIU7b+TsE87UonhRTKFGJILkD08jy51Ahae+00BhMG8Fmd1twOd0Qfqp5pjwRzNNvHQxlL8GVcvAajb3Wc1Q9b45Pv8VxPWDzFkm4GSeBcnDtGWDHL39pebfXHkRlL8BTXsBKN53AuVN3a+WbQHlnnqKN+9m6giNB4tWCjUiEcTldDB/dh5zFq3EAR2CzRFtwOdwQN+x5hh7h1nrY+c7UPE6lL/ecWZJ4x6zCOC258x54iBInw4Zp5mP7nF4cWqJeemat8FMvy4vMkdX24AAlb6e7bHU00ewEpkUakQiTMG4LBZeNpHCJWsDtwFfbAoMnm0OgH1boOINE3Aq3mz7bRqAA6Ww7RlzAEX7zqKw9MeUN/i3bMgKkc0AQ2E/n1CoIagsy4TiiqUmxOx8yz9Y/WBx/cweTEMuJONAPjy66rBfvqePYCUyaUq3SIQK2pulz2v2oyp/3TyS2rWs7U2qyDOVOVvv6DSw09FyZeE59RRMHg8pI8xS9UFUtLq8U/ALdtgKhRp6nbcR9qyCqg/8R8Oubm52mB3nswrM0f/ktvVkWpcrONyj1Q9u15iaSKR1arqgUCMSBL5m2PsZ3or3mf5sNuX1SRw8sBPMYoCZsdVmzZy4FOh3khl03HqkjOy1oNP9AoXGwssm9nqoCIUaekWjB3Z95A8w1SvAe6D7+xMGmvWRsgrMApAJA7q9tfV7Bl0/Wg3b75kcltapERF7OGMgbRLFe4ZRXv9Rt7dZOM2aOXVjmer4wj+FvFVsqj/o9DsJ3GMh9fhjXvXY67MoXLK2y9/2LcwbZOGStZydl9lrv/GHQg0BUV8Je0pajlXmY8166PJv1iK2b8uA9OkmzPQb3+Pw2iuPViWiKNSISK/o8Zo5KbMgYbcZi9NeUw1UvmuOVg4nJOeaAczusZA6xgSd1NFm3E8PHH6Bwt6fRRMKNRwRnxfqNncML3tK4EDZ4T83eZgJMOnTTZhx5x1TD1zBuCzOzsuMrnFI0mMKNSLSK3q8Zk7+L2DEn+HATtj9qTn2tHzcv6PjzZYP9m00x47FHV9LHNQScEZB8nDoM8y8oSYPM3tctUwz73HY6sVZNKFQQyeWZUJK7QYzkLd1baKar8yaRb7Gw38NZxy4x0H6NH+ISRoc8FJdTkdohD0JOQo1ItIrjnjNnMSBMOib5mhVXwm7V8Lez8GzpuX4svOuy2B6eg6UmsXaDuZKMjszJw8j48BEYOph6+/NWTRHvEhiIHgbob7crANzoNT/cd8Wf4Dp6vvandi+0G+CeTTYb4I53GPAGRu4mkWOkEKNiPSKgKyZk5AB2QXmaGX5oG4L7F0DteuhZp3/aKju+ut490PNl1DzJfnWa2TF/o2Kpv5YdH4M4sAiM6GO/F13QN1AiM8wdSRkQKwbYvqYR12xKSYsHcVCg8e0SCKYXpXmOjOVvnEPNB78cU/nANO6BcaRcsab2Wkpo6DveBNe0k6CpCFaZFFCjmY/iUivCuq05fpdpsehbivs32p6IeraHV5TQ+tUc7A6BBsHPsDBwqF3UuBefvg/z+E0ISemJeS0fexjBkw7XGZbCYfLf7RcLyrLZs7yCYDZaNRfQ8t09wmvUZCxzixO52swj3+a6/zhpXWF50BwxJhNS1OOM+El5ThzpI6CxMHaBkNspyndXVCoEbFHSCwwZ1nmcVZLwClau5vC4nTKD/hnU2XFVjM/+8GeBZoAKPJMpbDsJx2W/8+KrWJ+9sOBq8ERA4lZZsxRUnbLx0HtPg4244/02EhCmEJNFxRqRKS9LsOW1WQe1dTvNCGo9WNDlZmR1VQLzbXQvM/fbmo5mvdxyOnMXdVgOSmuG0tlUz8yYveQn7wGl8PX+UZnPLgSzSq7cX07fozt4jxxoAkuCRlBX9hQJNC0To2IyGF0PYsmzvRgJA068i9o+aB5v3lMZDWD5TUffd6Wttd/raXtsnxMdcWbmUPOlo/tz13xLY+wNH5F5HAUakREAsXhhNg+5hCRoFOfpIiIiEQEhRoRERGJCAo1IiIiEhEUakRERCQiKNSIiIhIRFCoERERkYigUCMiIiIRQaFGREREIoJCjYiIiESEqFpRuHWbq5qaGpsrERERkZ5qfd8+3HaVURVqamtrAcjJybG5EhERETlStbW1uN3ubl+Pql26fT4fZWVlpKSk4AjzzeFqamrIyclh+/btUb3juL4P+h600vdB34NW+j5E3vfAsixqa2vJzs7G6ex+5ExU9dQ4nU4GDx5sdxkBlZqaGhH/YI+Vvg/6HrTS90Hfg1b6PkTW9+BQPTStNFBYREREIoJCjYiIiEQEhZowFR8fz/z584mPj7e7FFvp+6DvQSt9H/Q9aKXvQ/R+D6JqoLCIiIhELvXUiIiISERQqBEREZGIoFAjIiIiEUGhRkRERCKCQk2EaWhoYMKECTgcDkpKSuwuJ6jOO+88hgwZQkJCAllZWVx++eWUlZXZXVZQbdmyhauvvprhw4eTmJjIiBEjmD9/Po2NjXaXFlS/+93vmDZtGklJSfTt29fucoLmgQceYNiwYSQkJDBlyhSKi4vtLimo3nvvPWbPnk12djYOh4OXXnrJ7pKCbsGCBZx88smkpKSQkZHBBRdcwPr16+0uK2gUaiLMbbfdRnZ2tt1l2GLmzJk8++yzrF+/nhdeeIFNmzbxve99z+6ygmrdunX4fD4eeugh1qxZwz333MODDz7IHXfcYXdpQdXY2MiFF17InDlz7C4laJ555hluvvlm5s+fz8qVKxk/fjznnnsulZWVdpcWNHV1dYwfP54HHnjA7lJs8+6773L99dfz0Ucf8cYbb9DU1MQ555xDXV2d3aUFhyUR45VXXrGOP/54a82aNRZgrVq1yu6SbLV48WLL4XBYjY2Ndpdiqz/+8Y/W8OHD7S7DFo899pjldrvtLiMo8vPzreuvv77t3Ov1WtnZ2daCBQtsrMo+gPXiiy/aXYbtKisrLcB699137S4lKNRTEyF27tzJNddcw5NPPklSUpLd5dhu9+7dPPXUU0ybNo3Y2Fi7y7GVx+MhLS3N7jKkFzU2NvLpp58ya9astmtOp5NZs2axfPlyGysTu3k8HoCo+RmgUBMBLMviyiuv5LrrrmPy5Ml2l2Or22+/neTkZPr378+2bdtYvHix3SXZauPGjdx3331ce+21dpcivWjXrl14vV4GDhzY4frAgQOpqKiwqSqxm8/n48Ybb+TUU09l3LhxdpcTFAo1IWzu3Lk4HI5DHuvWreO+++6jtraWefPm2V1ywPX0e9Dq1ltvZdWqVbz++uu4XC5++MMfYkXAotlH+n0AKC0tpaCggAsvvJBrrrnGpsoD52i+ByLR7Prrr2f16tU8/fTTdpcSNNomIYRVVVVRXV19yHtyc3O56KKLWLJkCQ6Ho+261+vF5XJx6aWX8sQTT/R2qb2mp9+DuLi4Ttd37NhBTk4Oy5YtY+rUqb1VYlAc6fehrKyMGTNmcMopp/D444/jdIb/7y9H82/h8ccf58Ybb2Tv3r29XJ29GhsbSUpK4vnnn+eCCy5ou37FFVewd+/eqOyxdDgcvPjiix2+H9HkhhtuYPHixbz33nsMHz7c7nKCJsbuAqR76enppKenH/a+v/zlL/z2t79tOy8rK+Pcc8/lmWeeYcqUKb1ZYq/r6fegKz6fDzDT3MPdkXwfSktLmTlzJpMmTeKxxx6LiEADx/ZvIdLFxcUxadIkli5d2vYm7vP5WLp0KTfccIO9xUlQWZbFT3/6U1588UXeeeedqAo0oFATEYYMGdLhvE+fPgCMGDGCwYMH21FS0K1YsYKPP/6Y6dOn069fPzZt2sSvfvUrRowYEfa9NEeitLSUGTNmMHToUO6++26qqqraXsvMzLSxsuDatm0bu3fvZtu2bXi93rY1m0aOHNn230ekufnmm7niiiuYPHky+fn53HvvvdTV1XHVVVfZXVrQ7Nu3j40bN7adb968mZKSEtLS0jr9nIxU119/PX//+99ZvHgxKSkpbWOq3G43iYmJNlcXBLbOvZJesXnz5qib0v35559bM2fOtNLS0qz4+Hhr2LBh1nXXXWft2LHD7tKC6rHHHrOALo9ocsUVV3T5PXj77bftLq1X3XfffdaQIUOsuLg4Kz8/3/roo4/sLimo3n777S7/f7/iiivsLi1ouvvv/7HHHrO7tKDQmBoRERGJCJHxsF1ERESinkKNiIiIRASFGhEREYkICjUiIiISERRqREREJCIo1IiIiEhEUKgRERGRiKBQIyIiIhFBoUZEREQigkKNiIiIRASFGhEREYkICjUiEraqqqrIzMzkzjvvbLu2bNky4uLiWLp0qY2ViYgdtKGliIS1V155hQsuuIBly5YxevRoJkyYwPnnn8+f//xnu0sTkSBTqBGRsHf99dfz5ptvMnnyZL744gs+/vhj4uPj7S5LRIJMoUZEwt6BAwcYN24c27dv59NPP+WEE06wuyQRsYHG1IhI2Nu0aRNlZWX4fD62bNlidzkiYhP11IhIWGtsbCQ/P58JEyYwevRo7r33Xr744gsyMjLsLk1EgkyhRkTC2q233srzzz/PZ599Rp8+fTjjjDNwu928/PLLdpcmIkGmx08iErbeeecd7r33Xp588klSU1NxOp08+eSTvP/++yxcuNDu8kQkyNRTIyIiIhFBPTUiIiISERRqREREJCIo1IiIiEhEUKgRERGRiKBQIyIiIhFBoUZEREQigkKNiIiIRASFGhEREYkICjUiIiISERRqREREJCIo1IiIiEhE+P/p7WR3SjN7HgAAAABJRU5ErkJggg==", "text/plain": [ "
" ] @@ -231,7 +234,7 @@ "fig, ax = plt.subplots()\n", "\n", "# Plot observed data\n", - "ax.plot(*X, 'o')\n", + "ax.plot(*Xdata, 'o')\n", "\n", "# Plot fitted ellipse\n", "x0, y0, a, b, theta = sol.beta\n", diff --git a/docs/reference/index.md b/docs/reference/index.md index e8c4ffd..6abd8c7 100644 --- a/docs/reference/index.md +++ b/docs/reference/index.md @@ -1,6 +1,6 @@ :::odrpack options: members: - - odr + - odr_fit - OdrResult - OdrStop \ No newline at end of file diff --git a/meson.build b/meson.build index 4cac98e..ffecbca 100644 --- a/meson.build +++ b/meson.build @@ -1,7 +1,7 @@ project( 'odrpack', ['cpp'], - version : '0.2.0', + version : '0.3.0', meson_version: '>=1.5.0', default_options: [ 'buildtype=release', @@ -36,11 +36,12 @@ nanobind_ext = python.extension_module( python.install_sources( [ 'src/odrpack/__init__.py', - 'src/odrpack/driver.py', + 'src/odrpack/__odrpack.pyi', + 'src/odrpack/odr_fortran.py', + 'src/odrpack/odr_scipy.py', 'src/odrpack/exceptions.py', 'src/odrpack/result.py', 'src/odrpack/py.typed', - 'src/odrpack/__odrpack.pyi', ], subdir: 'odrpack' ) diff --git a/src/odrpack/__init__.py b/src/odrpack/__init__.py index 9351592..ca698c2 100644 --- a/src/odrpack/__init__.py +++ b/src/odrpack/__init__.py @@ -9,8 +9,9 @@ to the observations of the dependent variable. """ -__version__ = '0.2.0' +__version__ = '0.3.0' -from odrpack.driver import * from odrpack.exceptions import * +from odrpack.odr_fortran import * +from odrpack.odr_scipy import * from odrpack.result import * diff --git a/src/odrpack/__odrpack.py.cpp b/src/odrpack/__odrpack.py.cpp index 674b2ef..2981063 100644 --- a/src/odrpack/__odrpack.py.cpp +++ b/src/odrpack/__odrpack.py.cpp @@ -29,20 +29,18 @@ class SelfCleaningPyObject { }; /* -Wrapper for the C-interface of the Fortran ODR routine. This wrapper is -intentionally very thin, with all argument checks and array dimension -calculations delegated to the companion Python caller, which serves as the entry -point for all function calls. - -Some arguments have a default value of `nullptr` — this is by design, as the -Fortran code automatically interprets `nullptr` as an absent optional argument. -This approach avoids the redundant definition of default values in multiple -places. +Wrapper for the C-interface of the Fortran ODR routine. This wrapper is intentionally very +thin, with all argument checks and array dimension calculations delegated to the companion +Python caller, which serves as the entry point for all function calls. + +Some arguments have a default value of `nullptr` — this is by design, as the Fortran code +automatically interprets `nullptr` as an absent optional argument. This approach avoids the +redundant definition of default values in multiple places. */ int odr_wrapper(int n, int m, + int q, int npar, - int nq, int ldwe, int ld2we, int ldwd, @@ -67,7 +65,7 @@ int odr_wrapper(int n, std::optional> scld, std::optional> lower, std::optional> upper, - std::optional> work, + std::optional> rwork, std::optional> iwork, std::optional job, std::optional ndigit, @@ -101,7 +99,7 @@ int odr_wrapper(int n, auto lower_ptr = lower ? lower.value().data() : nullptr; auto upper_ptr = upper ? upper.value().data() : nullptr; - auto work_ptr = work ? work.value().data() : nullptr; + auto rwork_ptr = rwork ? rwork.value().data() : nullptr; auto iwork_ptr = iwork ? iwork.value().data() : nullptr; auto job_ptr = job ? &job.value() : nullptr; @@ -112,9 +110,9 @@ int odr_wrapper(int n, auto maxit_ptr = maxit ? &maxit.value() : nullptr; auto iprint_ptr = iprint ? &iprint.value() : nullptr; - int lwork = 1; + int lrwork = 1; int liwork = 1; - if (work) lwork = work.value().size(); + if (rwork) lrwork = rwork.value().size(); if (iwork) liwork = iwork.value().size(); // Build static pointers to the Python functions @@ -134,8 +132,7 @@ int odr_wrapper(int n, // Define the overall user-supplied model function 'fcn' odrpack_fcn_t fcn = nullptr; - fcn = [](const int *n, const int *m, const int *npar, const int *nq, - const int *ldn, const int *ldm, const int *ldnp, const double beta[], + fcn = [](const int *n, const int *m, const int *q, const int *npar, const double beta[], const double xplusd[], const int ifixb[], const int ifixx[], const int *ldifx, const int *ideval, double f[], double fjacb[], double fjacd[], int *istop) { @@ -153,7 +150,7 @@ int odr_wrapper(int n, auto f_object = fcn_f_holder(beta_ndarray, xplusd_ndarray); auto f_ndarray = nb::cast>(f_object); auto f_ndarray_ptr = f_ndarray.data(); - for (auto i = 0; i < (*nq) * (*ldn); i++) { + for (auto i = 0; i < (*q) * (*n); i++) { f[i] = f_ndarray_ptr[i]; } } @@ -163,7 +160,7 @@ int odr_wrapper(int n, auto fjacb_object = fcn_fjacb_holder(beta_ndarray, xplusd_ndarray); auto fjacb_ndarray = nb::cast>(fjacb_object); auto fjacb_ndarray_ptr = fjacb_ndarray.data(); - for (auto i = 0; i < (*nq) * (*ldnp) * (*ldn); i++) { + for (auto i = 0; i < (*q) * (*npar) * (*n); i++) { fjacb[i] = fjacb_ndarray_ptr[i]; } } @@ -173,7 +170,7 @@ int odr_wrapper(int n, auto fjacd_object = fcn_fjacd_holder(beta_ndarray, xplusd_ndarray); auto fjacd_ndarray = nb::cast>(fjacd_object); auto fjacd_ndarray_ptr = fjacd_ndarray.data(); - for (auto i = 0; i < (*nq) * (*ldnp) * (*ldn); i++) { + for (auto i = 0; i < (*q) * (*npar) * (*n); i++) { fjacd[i] = fjacd_ndarray_ptr[i]; } } @@ -213,10 +210,10 @@ int odr_wrapper(int n, // Call the C function int info = -1; - odr_long_c(fcn, &n, &m, &npar, &nq, &ldwe, &ld2we, &ldwd, &ld2wd, &ldifx, - &ldstpd, &ldscld, &lwork, &liwork, beta_ptr, y_ptr, x_ptr, we_ptr, + odr_long_c(fcn, &n, &m, &q, &npar, &ldwe, &ld2we, &ldwd, &ld2wd, &ldifx, + &ldstpd, &ldscld, &lrwork, &liwork, beta_ptr, y_ptr, x_ptr, we_ptr, wd_ptr, ifixb_ptr, ifixx_ptr, stpb_ptr, stpd_ptr, sclb_ptr, - scld_ptr, delta_ptr, lower_ptr, upper_ptr, work_ptr, iwork_ptr, + scld_ptr, delta_ptr, lower_ptr, upper_ptr, rwork_ptr, iwork_ptr, job_ptr, ndigit_ptr, taufac_ptr, sstol_ptr, partol_ptr, maxit_ptr, iprint_ptr, &lunerr, &lunrpt, &info); @@ -245,14 +242,14 @@ n : int Number of observations. m : int Number of columns in the independent variable data. +q : int + Number of responses per observation. npar : int Number of function parameters. -nq : int - Number of responses per observation. ldwe : int Leading dimension of the `we` array, must be in `{1, n}`. ld2we : int - Second dimension of the `we` array, must be in `{1, nq}`. + Second dimension of the `we` array, must be in `{1, q}`. ldwd : int Leading dimension of the `wd` array, must be in `{1, n}`. ld2wd : int @@ -274,13 +271,13 @@ fjacd : Callable beta : np.ndarray[float64] Array of function parameters with shape `(npar)`. y : np.ndarray[float64] - Dependent variables with shape `(nq, n)`. Ignored for implicit models. + Dependent variables with shape `(q, n)`. Ignored for implicit models. x : np.ndarray[float64] Explanatory variables with shape `(m, n)`. delta : np.ndarray[float64] Initial errors in `x` data with shape `(m, n)`. we : np.ndarray[float64], optional - Weights for `epsilon` with shape `(nq, ld2we, ldwe)`. Default is None. + Weights for `epsilon` with shape `(q, ld2we, ldwe)`. Default is None. wd : np.ndarray[float64], optional Weights for `delta` with shape `(m, ld2wd, ldwd)`. Default is None. ifixb : np.ndarray[int32], optional @@ -299,7 +296,7 @@ lower : np.ndarray[float64], optional Lower bounds for `beta`. Default is None. upper : np.ndarray[float64], optional Upper bounds for `beta`. Default is None. -work : np.ndarray[float64], optional +rwork : np.ndarray[float64], optional Real work space. Default is None. iwork : np.ndarray[int32], optional Integer work space. Default is None. @@ -332,7 +329,7 @@ Notes - Ensure all array dimensions and functions are consistent with the provided arguments. - Input arrays will automatically be made contiguous and cast to the correct type if necessary. )doc", - nb::arg("n"), nb::arg("m"), nb::arg("npar"), nb::arg("nq"), + nb::arg("n"), nb::arg("m"), nb::arg("q"), nb::arg("npar"), nb::arg("ldwe"), nb::arg("ld2we"), nb::arg("ldwd"), nb::arg("ld2wd"), nb::arg("ldifx"), nb::arg("ldstpd"), nb::arg("ldscld"), nb::arg("f"), @@ -352,7 +349,7 @@ Notes nb::arg("scld").none() = nullptr, nb::arg("lower").none() = nullptr, nb::arg("upper").none() = nullptr, - nb::arg("work").none() = nullptr, + nb::arg("rwork").none() = nullptr, nb::arg("iwork").none() = nullptr, nb::arg("job").none() = nullptr, nb::arg("ndigit").none() = nullptr, @@ -367,11 +364,11 @@ Notes // Calculate the dimensions of the workspace arrays m.def( "workspace_dimensions", - [](int n, int m, int npar, int nq, bool isodr) { - int lwork = 0; + [](int n, int m, int q, int npar, bool isodr) { + int lrwork = 0; int liwork = 0; - workspace_dimensions_c(&n, &m, &npar, &nq, &isodr, &lwork, &liwork); - return nb::make_tuple(lwork, liwork); + workspace_dimensions_c(&n, &m, &q, &npar, &isodr, &lrwork, &liwork); + return nb::make_tuple(lrwork, liwork); }, R"doc( Calculate the dimensions of the workspace arrays. @@ -382,51 +379,51 @@ n : int Number of observations. m : int Number of columns of data in the explanatory variable. +q : int + Number of responses per observation. npar : int Number of function parameters. -nq : int - Number of responses per observation. isodr : bool Variable designating whether the solution is by ODR (`True`) or by OLS (`False`). Returns ------- tuple[int, int] - A tuple containing the lengths of the work arrays (`lwork`, `liwork`). + A tuple containing the lengths of the work arrays (`lrwork`, `liwork`). )doc", - nb::arg("n"), nb::arg("m"), nb::arg("npar"), nb::arg("nq"), + nb::arg("n"), nb::arg("m"), nb::arg("q"), nb::arg("npar"), nb::arg("isodr")); // Get storage locations within the integer work space m.def( - "diwinf", - [](int m, int npar, int nq) { - iworkidx_t iworkidx = {}; - diwinf_c(&m, &npar, &nq, &iworkidx); + "loc_iwork", + [](int m, int q, int npar) { + iworkidx_t iwi = {}; + loc_iwork_c(&m, &q, &npar, &iwi); std::map result; - result["msgb"] = iworkidx.msgb; - result["msgd"] = iworkidx.msgd; - result["ifix2"] = iworkidx.ifix2; - result["istop"] = iworkidx.istop; - result["nnzw"] = iworkidx.nnzw; - result["npp"] = iworkidx.npp; - result["idf"] = iworkidx.idf; - result["job"] = iworkidx.job; - result["iprin"] = iworkidx.iprin; - result["luner"] = iworkidx.luner; - result["lunrp"] = iworkidx.lunrp; - result["nrow"] = iworkidx.nrow; - result["ntol"] = iworkidx.ntol; - result["neta"] = iworkidx.neta; - result["maxit"] = iworkidx.maxit; - result["niter"] = iworkidx.niter; - result["nfev"] = iworkidx.nfev; - result["njev"] = iworkidx.njev; - result["int2"] = iworkidx.int2; - result["irank"] = iworkidx.irank; - result["ldtt"] = iworkidx.ldtt; - result["bound"] = iworkidx.bound; - result["liwkmn"] = iworkidx.liwkmn; + result["msgb"] = iwi.msgb; + result["msgd"] = iwi.msgd; + result["ifix2"] = iwi.ifix2; + result["istop"] = iwi.istop; + result["nnzw"] = iwi.nnzw; + result["npp"] = iwi.npp; + result["idf"] = iwi.idf; + result["job"] = iwi.job; + result["iprin"] = iwi.iprin; + result["luner"] = iwi.luner; + result["lunrp"] = iwi.lunrp; + result["nrow"] = iwi.nrow; + result["ntol"] = iwi.ntol; + result["neta"] = iwi.neta; + result["maxit"] = iwi.maxit; + result["niter"] = iwi.niter; + result["nfev"] = iwi.nfev; + result["njev"] = iwi.njev; + result["int2"] = iwi.int2; + result["irank"] = iwi.irank; + result["ldtt"] = iwi.ldtt; + result["bound"] = iwi.bound; + result["liwkmn"] = iwi.liwkmn; return result; }, R"doc( @@ -436,77 +433,77 @@ Parameters ---------- m : int Number of columns of data in the explanatory variable. +q : int + Number of responses per observation. npar : int Number of function parameters. -nq : int - Number of responses per observation. Returns ------- dict[str, int] A dictionary containing the 0-based indexes of the integer work array. )doc", - nb::arg("m"), nb::arg("npar"), nb::arg("nq")); + nb::arg("m"), nb::arg("q"), nb::arg("npar")); // Get storage locations within the real work space m.def( - "dwinf", - [](int n, int m, int npar, int nq, int ldwe, int ld2we, bool isodr) { - workidx_t workidx = {}; - dwinf_c(&n, &m, &npar, &nq, &ldwe, &ld2we, &isodr, &workidx); + "loc_rwork", + [](int n, int m, int q, int npar, int ldwe, int ld2we, bool isodr) { + rworkidx_t rwi = {}; + loc_rwork_c(&n, &m, &q, &npar, &ldwe, &ld2we, &isodr, &rwi); std::map result; - result["delta"] = workidx.delta; - result["eps"] = workidx.eps; - result["xplus"] = workidx.xplus; - result["fn"] = workidx.fn; - result["sd"] = workidx.sd; - result["vcv"] = workidx.vcv; - result["rvar"] = workidx.rvar; - result["wss"] = workidx.wss; - result["wssde"] = workidx.wssde; - result["wssep"] = workidx.wssep; - result["rcond"] = workidx.rcond; - result["eta"] = workidx.eta; - result["olmav"] = workidx.olmav; - result["tau"] = workidx.tau; - result["alpha"] = workidx.alpha; - result["actrs"] = workidx.actrs; - result["pnorm"] = workidx.pnorm; - result["rnors"] = workidx.rnors; - result["prers"] = workidx.prers; - result["partl"] = workidx.partl; - result["sstol"] = workidx.sstol; - result["taufc"] = workidx.taufc; - result["epsma"] = workidx.epsma; - result["beta0"] = workidx.beta0; - result["betac"] = workidx.betac; - result["betas"] = workidx.betas; - result["betan"] = workidx.betan; - result["s"] = workidx.s; - result["ss"] = workidx.ss; - result["ssf"] = workidx.ssf; - result["qraux"] = workidx.qraux; - result["u"] = workidx.u; - result["fs"] = workidx.fs; - result["fjacb"] = workidx.fjacb; - result["we1"] = workidx.we1; - result["diff"] = workidx.diff; - result["delts"] = workidx.delts; - result["deltn"] = workidx.deltn; - result["t"] = workidx.t; - result["tt"] = workidx.tt; - result["omega"] = workidx.omega; - result["fjacd"] = workidx.fjacd; - result["wrk1"] = workidx.wrk1; - result["wrk2"] = workidx.wrk2; - result["wrk3"] = workidx.wrk3; - result["wrk4"] = workidx.wrk4; - result["wrk5"] = workidx.wrk5; - result["wrk6"] = workidx.wrk6; - result["wrk7"] = workidx.wrk7; - result["lower"] = workidx.lower; - result["upper"] = workidx.upper; - result["lwkmn"] = workidx.lwkmn; + result["delta"] = rwi.delta; + result["eps"] = rwi.eps; + result["xplus"] = rwi.xplus; + result["fn"] = rwi.fn; + result["sd"] = rwi.sd; + result["vcv"] = rwi.vcv; + result["rvar"] = rwi.rvar; + result["wss"] = rwi.wss; + result["wssde"] = rwi.wssde; + result["wssep"] = rwi.wssep; + result["rcond"] = rwi.rcond; + result["eta"] = rwi.eta; + result["olmav"] = rwi.olmav; + result["tau"] = rwi.tau; + result["alpha"] = rwi.alpha; + result["actrs"] = rwi.actrs; + result["pnorm"] = rwi.pnorm; + result["rnors"] = rwi.rnors; + result["prers"] = rwi.prers; + result["partl"] = rwi.partl; + result["sstol"] = rwi.sstol; + result["taufc"] = rwi.taufc; + result["epsma"] = rwi.epsma; + result["beta0"] = rwi.beta0; + result["betac"] = rwi.betac; + result["betas"] = rwi.betas; + result["betan"] = rwi.betan; + result["s"] = rwi.s; + result["ss"] = rwi.ss; + result["ssf"] = rwi.ssf; + result["qraux"] = rwi.qraux; + result["u"] = rwi.u; + result["fs"] = rwi.fs; + result["fjacb"] = rwi.fjacb; + result["we1"] = rwi.we1; + result["diff"] = rwi.diff; + result["delts"] = rwi.delts; + result["deltn"] = rwi.deltn; + result["t"] = rwi.t; + result["tt"] = rwi.tt; + result["omega"] = rwi.omega; + result["fjacd"] = rwi.fjacd; + result["wrk1"] = rwi.wrk1; + result["wrk2"] = rwi.wrk2; + result["wrk3"] = rwi.wrk3; + result["wrk4"] = rwi.wrk4; + result["wrk5"] = rwi.wrk5; + result["wrk6"] = rwi.wrk6; + result["wrk7"] = rwi.wrk7; + result["lower"] = rwi.lower; + result["upper"] = rwi.upper; + result["lrwkmn"] = rwi.lrwkmn; return result; }, R"doc( @@ -518,10 +515,10 @@ n : int Number of observations. m : int Number of columns of data in the explanatory variable. +q : int + Number of responses per observation. npar : int Number of function parameters. -nq : int - Number of responses per observation. ldwe : int Leading dimension of the `we` array. ld2we : int @@ -534,6 +531,6 @@ Returns dict[str, int] A dictionary containing the 0-based indexes of the real work array. )doc", - nb::arg("n"), nb::arg("m"), nb::arg("npar"), nb::arg("nq"), + nb::arg("n"), nb::arg("m"), nb::arg("q"), nb::arg("npar"), nb::arg("ldwe"), nb::arg("ld2we"), nb::arg("isodr")); } \ No newline at end of file diff --git a/src/odrpack/__odrpack.pyi b/src/odrpack/__odrpack.pyi index 8936ed9..3548f61 100644 --- a/src/odrpack/__odrpack.pyi +++ b/src/odrpack/__odrpack.pyi @@ -4,7 +4,7 @@ from typing import Annotated from numpy.typing import ArrayLike -def diwinf(m: int, npar: int, nq: int) -> dict[str, int]: +def loc_iwork(m: int, q: int, npar: int) -> dict[str, int]: """ Get storage locations within the integer work space. @@ -12,10 +12,10 @@ def diwinf(m: int, npar: int, nq: int) -> dict[str, int]: ---------- m : int Number of columns of data in the explanatory variable. + q : int + Number of responses per observation. npar : int Number of function parameters. - nq : int - Number of responses per observation. Returns ------- @@ -23,7 +23,7 @@ def diwinf(m: int, npar: int, nq: int) -> dict[str, int]: A dictionary containing the 0-based indexes of the integer work array. """ -def dwinf(n: int, m: int, npar: int, nq: int, ldwe: int, ld2we: int, isodr: bool) -> dict[str, int]: +def loc_rwork(n: int, m: int, q: int, npar: int, ldwe: int, ld2we: int, isodr: bool) -> dict[str, int]: """ Get storage locations within the real work space. @@ -33,10 +33,10 @@ def dwinf(n: int, m: int, npar: int, nq: int, ldwe: int, ld2we: int, isodr: bool Number of observations. m : int Number of columns of data in the explanatory variable. + q : int + Number of responses per observation. npar : int Number of function parameters. - nq : int - Number of responses per observation. ldwe : int Leading dimension of the `we` array. ld2we : int @@ -50,7 +50,7 @@ def dwinf(n: int, m: int, npar: int, nq: int, ldwe: int, ld2we: int, isodr: bool A dictionary containing the 0-based indexes of the real work array. """ -def odr(n: int, m: int, npar: int, nq: int, ldwe: int, ld2we: int, ldwd: int, ld2wd: int, ldifx: int, ldstpd: int, ldscld: int, f: Callable, fjacb: Callable, fjacd: Callable, beta: Annotated[ArrayLike, dict(dtype='float64', order='C')], y: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)], x: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)], delta: Annotated[ArrayLike, dict(dtype='float64', order='C')], we: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, wd: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, ifixb: Annotated[ArrayLike, dict(dtype='int32', order='C', writable=False)] | None = None, ifixx: Annotated[ArrayLike, dict(dtype='int32', order='C', writable=False)] | None = None, stpb: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, stpd: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, sclb: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, scld: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, lower: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, upper: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, work: Annotated[ArrayLike, dict(dtype='float64', order='C')] | None = None, iwork: Annotated[ArrayLike, dict(dtype='int32', order='C')] | None = None, job: int | None = None, ndigit: int | None = None, taufac: float | None = None, sstol: float | None = None, partol: float | None = None, maxit: int | None = None, iprint: int | None = None, errfile: str | None = None, rptfile: str | None = None) -> int: +def odr(n: int, m: int, q: int, npar: int, ldwe: int, ld2we: int, ldwd: int, ld2wd: int, ldifx: int, ldstpd: int, ldscld: int, f: Callable, fjacb: Callable, fjacd: Callable, beta: Annotated[ArrayLike, dict(dtype='float64', order='C')], y: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)], x: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)], delta: Annotated[ArrayLike, dict(dtype='float64', order='C')], we: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, wd: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, ifixb: Annotated[ArrayLike, dict(dtype='int32', order='C', writable=False)] | None = None, ifixx: Annotated[ArrayLike, dict(dtype='int32', order='C', writable=False)] | None = None, stpb: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, stpd: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, sclb: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, scld: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, lower: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, upper: Annotated[ArrayLike, dict(dtype='float64', order='C', writable=False)] | None = None, rwork: Annotated[ArrayLike, dict(dtype='float64', order='C')] | None = None, iwork: Annotated[ArrayLike, dict(dtype='int32', order='C')] | None = None, job: int | None = None, ndigit: int | None = None, taufac: float | None = None, sstol: float | None = None, partol: float | None = None, maxit: int | None = None, iprint: int | None = None, errfile: str | None = None, rptfile: str | None = None) -> int: """ C++ wrapper for the Orthogonal Distance Regression (ODR) routine. @@ -60,14 +60,14 @@ def odr(n: int, m: int, npar: int, nq: int, ldwe: int, ld2we: int, ldwd: int, ld Number of observations. m : int Number of columns in the independent variable data. + q : int + Number of responses per observation. npar : int Number of function parameters. - nq : int - Number of responses per observation. ldwe : int Leading dimension of the `we` array, must be in `{1, n}`. ld2we : int - Second dimension of the `we` array, must be in `{1, nq}`. + Second dimension of the `we` array, must be in `{1, q}`. ldwd : int Leading dimension of the `wd` array, must be in `{1, n}`. ld2wd : int @@ -89,13 +89,13 @@ def odr(n: int, m: int, npar: int, nq: int, ldwe: int, ld2we: int, ldwd: int, ld beta : np.ndarray[float64] Array of function parameters with shape `(npar)`. y : np.ndarray[float64] - Dependent variables with shape `(nq, n)`. Ignored for implicit models. + Dependent variables with shape `(q, n)`. Ignored for implicit models. x : np.ndarray[float64] Explanatory variables with shape `(m, n)`. delta : np.ndarray[float64] Initial errors in `x` data with shape `(m, n)`. we : np.ndarray[float64], optional - Weights for `epsilon` with shape `(nq, ld2we, ldwe)`. Default is None. + Weights for `epsilon` with shape `(q, ld2we, ldwe)`. Default is None. wd : np.ndarray[float64], optional Weights for `delta` with shape `(m, ld2wd, ldwd)`. Default is None. ifixb : np.ndarray[int32], optional @@ -114,7 +114,7 @@ def odr(n: int, m: int, npar: int, nq: int, ldwe: int, ld2we: int, ldwd: int, ld Lower bounds for `beta`. Default is None. upper : np.ndarray[float64], optional Upper bounds for `beta`. Default is None. - work : np.ndarray[float64], optional + rwork : np.ndarray[float64], optional Real work space. Default is None. iwork : np.ndarray[int32], optional Integer work space. Default is None. @@ -148,7 +148,7 @@ def odr(n: int, m: int, npar: int, nq: int, ldwe: int, ld2we: int, ldwd: int, ld - Input arrays will automatically be made contiguous and cast to the correct type if necessary. """ -def workspace_dimensions(n: int, m: int, npar: int, nq: int, isodr: bool) -> tuple: +def workspace_dimensions(n: int, m: int, q: int, npar: int, isodr: bool) -> tuple: """ Calculate the dimensions of the workspace arrays. @@ -158,15 +158,15 @@ def workspace_dimensions(n: int, m: int, npar: int, nq: int, isodr: bool) -> tup Number of observations. m : int Number of columns of data in the explanatory variable. + q : int + Number of responses per observation. npar : int Number of function parameters. - nq : int - Number of responses per observation. isodr : bool Variable designating whether the solution is by ODR (`True`) or by OLS (`False`). Returns ------- tuple[int, int] - A tuple containing the lengths of the work arrays (`lwork`, `liwork`). + A tuple containing the lengths of the work arrays (`lrwork`, `liwork`). """ diff --git a/src/odrpack/exceptions.py b/src/odrpack/exceptions.py index 5c4790b..3891635 100644 --- a/src/odrpack/exceptions.py +++ b/src/odrpack/exceptions.py @@ -3,9 +3,9 @@ class OdrStop(Exception): """ - Exception to stop the fitting process. + Exception to stop the regression. This exception can be raised in the model function or its Jacobians to - instruct the `odr` routine to stop fitting. + stop the regression process. """ pass diff --git a/src/odrpack/driver.py b/src/odrpack/odr_fortran.py similarity index 85% rename from src/odrpack/driver.py rename to src/odrpack/odr_fortran.py index ae9af93..4f53218 100644 --- a/src/odrpack/driver.py +++ b/src/odrpack/odr_fortran.py @@ -1,9 +1,10 @@ +import warnings from typing import Callable import numpy as np from numpy.typing import NDArray -from odrpack.__odrpack import diwinf, dwinf +from odrpack.__odrpack import loc_iwork, loc_rwork from odrpack.__odrpack import odr as _odr from odrpack.__odrpack import workspace_dimensions from odrpack.result import OdrResult @@ -40,7 +41,7 @@ def odr(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float stpd: NDArray[np.float64] | None = None, sclb: NDArray[np.float64] | None = None, scld: NDArray[np.float64] | None = None, - work: NDArray[np.float64] | None = None, + rwork: NDArray[np.float64] | None = None, iwork: NDArray[np.int32] | None = None, ) -> OdrResult: r"""Solve a weighted orthogonal distance regression (ODR) problem, also @@ -56,7 +57,7 @@ def odr(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float within the bounds specified by arguments `lower` and `upper` (if they are specified). y : NDArray[np.float64] - Array of shape `(n,)` or `(nq, n)` containing the values of the response + Array of shape `(n,)` or `(q, n)` containing the values of the response variable(s). When the model is explicit, the user must specify a value for each element of `y`. If some responses of some observations are actually missing, then the user can set the corresponding weight in @@ -69,13 +70,13 @@ def odr(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float we : float | NDArray[np.float64] | None Scalar or array specifying how the errors on `y` are to be weighted. If `we` is a scalar, then it is used for all data points. If `we` is - an array of shape `(n,)` and `nq==1`, then `we[i]` represents the weight - for `y[i]`. If `we` is an array of shape `(nq)`, then it represents the + an array of shape `(n,)` and `q==1`, then `we[i]` represents the weight + for `y[i]`. If `we` is an array of shape `(q)`, then it represents the diagonal of the covariant weighting matrix for all data points. If `we` - is an array of shape `(nq, nq)`, then it represents the full covariant + is an array of shape `(q, q)`, then it represents the full covariant weighting matrix for all data points. If `we` is an array of shape - `(nq, n)`, then `we[:, i]` represents the diagonal of the covariant - weighting matrix for `y[:, i]`. If `we` is an array of shape `(nq, nq, n)`, + `(q, n)`, then `we[:, i]` represents the diagonal of the covariant + weighting matrix for `y[:, i]`. If `we` is an array of shape `(q, q, n)`, then `we[:, :, i]` represents the full covariant weighting matrix for `y[:, i]`. For a comprehensive description of the options, refer to page 25 of the ODRPACK95 guide. By default, `we` is set to one for all `y` @@ -97,13 +98,13 @@ def odr(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float fjacb : Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float64]] | None Jacobian of the function to be fitted with respect to `beta`, with the signature `fjacb(beta, x)`. It must return an array with shape - `(n, npar, nq)` or compatible. To activate this option, `job` must be + `(n, npar, q)` or compatible. To activate this option, `job` must be set accordingly. By default, the Jacobian is evaluated numerically according to the finite difference scheme defined in `job`. fjacd : Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float64]] | None Jacobian of the function to be fitted with respect to `delta`, with the signature `fjacd(beta, x)`. It must return an array with shape - `(n, m, nq)` or compatible. To activate this option, `job` must be + `(n, m, q)` or compatible. To activate this option, `job` must be set accordingly. By default, the Jacobian is evaluated numerically according to the finite difference scheme defined in `job`. ifixb : NDArray[np.int32] | None @@ -139,7 +140,8 @@ def odr(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float forward finite difference, and covariance matrix computed using Jacobian matrices recomputed at the final solution. Another common option is 20, corresponding to an explicit orthogonal distance regression with - user-supplied jacobians `fjacb` and `fjacd`. To initialize `delta0` with + user-supplied jacobians `fjacb` and `fjacd`. For an implicit orthogonal + distance regression, `job` must be set to 1. To initialize `delta0` with the user supplied values, the 4th digit of `job` must be set to 1, e.g. 1000. To restart a previous run, the 5th digit of `job` must be set to 1, e.g. 10000. For a comprehensive description of the options, refer to @@ -211,7 +213,7 @@ def odr(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float should not be confused with the weighting matrices `we` and `wd`. By default, `scld` is set internally based on the relative magnitudes of `x`. For further details, refer to pages 32 and 85 of the ODRPACK95 guide. - work : NDArray[np.float64] | None + rwork : NDArray[np.float64] | None Array containing the real-valued internal state of the odrpack solver. It is only required for a restart (see `job`), in which case it must be set to the state of the previous run. @@ -254,6 +256,13 @@ def odr(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float array([1.63337057, 0.9 ]) """ + # Future deprecation warning + warnings.warn( + "This function is deprecated and will be removed in a future version. " + "Please use `odr_fit` instead", + FutureWarning + ) + # Interpret job is_odr = _get_digit(job, 1) < 2 has_jac = _get_digit(job, 2) > 1 @@ -270,12 +279,12 @@ def odr(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float f"`x` must be a rank-1 array of shape `(n,)` or a rank-2 array of shape `(m, n)`, but has shape {x.shape}.") if y.ndim == 1: - nq = 1 + q = 1 elif y.ndim == 2: - nq = y.shape[0] + q = y.shape[0] else: raise ValueError( - f"`y` must be a rank-1 array of shape `(n,)` or a rank-2 array of shape `(nq, n)`, but has shape {y.shape}.") + f"`y` must be a rank-1 array of shape `(n,)` or a rank-2 array of shape `(q, n)`, but has shape {y.shape}.") if x.shape[-1] == y.shape[-1]: n = x.shape[-1] @@ -318,7 +327,7 @@ def odr(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float raise ValueError("`delta0` must have the same shape as `x`.") delta = delta0.copy() elif not has_delta0 and delta0 is None: - delta = np.zeros_like(x) + delta = np.zeros(x.shape, dtype=np.float64) else: raise ValueError("Inconsistent arguments for `job` and `delta0`.") @@ -372,23 +381,23 @@ def odr(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float if isinstance(we, (float, int)): ldwe = 1 ld2we = 1 - we = np.full((nq,), we, dtype=np.float64) + we = np.full((q,), we, dtype=np.float64) elif isinstance(we, np.ndarray): - if we.shape == (nq,): + if we.shape == (q,): ldwe = 1 ld2we = 1 - elif we.shape == (nq, nq): + elif we.shape == (q, q): ldwe = 1 - ld2we = nq - elif we.shape == (nq, n) or (we.shape == (n,) and nq == 1): + ld2we = q + elif we.shape == (q, n) or (we.shape == (n,) and q == 1): ldwe = n ld2we = 1 - elif we.shape in ((nq, 1, 1), (nq, 1, n), (nq, nq, 1), (nq, nq, n)): + elif we.shape in ((q, 1, 1), (q, 1, n), (q, q, 1), (q, q, n)): ldwe = we.shape[2] ld2we = we.shape[1] else: raise ValueError( - r"`we` must be a array of shape `(nq,)`, `(n,)`, `(nq, nq)`, `(nq, n)`, `(nq, 1, 1)`, `(nq, 1, n)`, `(nq, nq, 1)`, or `(nq, nq, n)`. See page 25 of the ODRPACK95 User Guide.") + r"`we` must be a array of shape `(q,)`, `(n,)`, `(q, q)`, `(q, n)`, `(q, 1, 1)`, `(q, 1, n)`, `(q, q, 1)`, or `(q, q, n)`. See page 25 of the ODRPACK95 User Guide.") else: raise TypeError("`we` must be a float or an array.") else: @@ -433,9 +442,9 @@ def fdummy(beta, x): return np.array([np.nan]) # will never be called if has_jac and fjacb is not None: fjacb0 = fjacb(beta0, x) - if fjacb0.shape[-1] != n or fjacb0.size != n*npar*nq: + if fjacb0.shape[-1] != n or fjacb0.size != n*npar*q: raise ValueError( - "Function `fjacb` must return an array with shape `(n, npar, nq)` or compatible.") + "Function `fjacb` must return an array with shape `(n, npar, q)` or compatible.") elif not has_jac and fjacb is None: fjacb = fdummy else: @@ -443,21 +452,21 @@ def fdummy(beta, x): return np.array([np.nan]) # will never be called if has_jac and fjacd is not None: fjacd0 = fjacd(beta0, x) - if fjacd0.shape[-1] != n or fjacd0.size != n*m*nq: + if fjacd0.shape[-1] != n or fjacd0.size != n*m*q: raise ValueError( - "Function `fjacd` must return an array with shape `(n, m, nq)` or compatible.") + "Function `fjacd` must return an array with shape `(n, m, q)` or compatible.") elif not has_jac and fjacd is None: fjacd = fdummy else: raise ValueError("Inconsistent arguments for `job` and `fjacd`.") # Check/allocate work arrays - lwork, liwork = workspace_dimensions(n, m, npar, nq, is_odr) - if (not is_restart) and (work is None) and (iwork is None): - work = np.zeros(lwork, dtype=np.float64) + lwork, liwork = workspace_dimensions(n, m, q, npar, is_odr) + if (not is_restart) and (rwork is None) and (iwork is None): + rwork = np.zeros(lwork, dtype=np.float64) iwork = np.zeros(liwork, dtype=np.int32) - elif is_restart and (work is not None) and (iwork is not None): - if work.size != lwork: + elif is_restart and (rwork is not None) and (iwork is not None): + if rwork.size != lwork: raise ValueError( "Work array `work` does not have the correct length.") if iwork.size != liwork: @@ -469,7 +478,7 @@ def fdummy(beta, x): return np.array([np.nan]) # will never be called # Call the ODRPACK95 routine # Note: beta, delta, work, and iwork are modified in place - info = _odr(n=n, m=m, npar=npar, nq=nq, + info = _odr(n=n, m=m, q=q, npar=npar, ldwe=ldwe, ld2we=ld2we, ldwd=ldwd, ld2wd=ld2wd, ldifx=ldifx, @@ -479,27 +488,27 @@ def fdummy(beta, x): return np.array([np.nan]) # will never be called delta=delta, we=we, wd=wd, ifixb=ifixb, ifixx=ifixx, lower=lower, upper=upper, - work=work, iwork=iwork, + rwork=rwork, iwork=iwork, job=job, ndigit=ndigit, taufac=taufac, sstol=sstol, partol=partol, maxit=maxit, iprint=iprint, errfile=errfile, rptfile=rptfile ) # Indexes of integer and real work arrays - iwork_idx: dict[str, int] = diwinf(m, npar, nq) - work_idx: dict[str, int] = dwinf(n, m, npar, nq, ldwe, ld2we, is_odr) + iwork_idx: dict[str, int] = loc_iwork(m, q, npar) + rwork_idx: dict[str, int] = loc_rwork(n, m, q, npar, ldwe, ld2we, is_odr) # Return the result # Extract results without messing up the original work arrays - i0_eps = work_idx['eps'] - eps = work[i0_eps:i0_eps+y.size].copy() + i0_eps = rwork_idx['eps'] + eps = rwork[i0_eps:i0_eps+y.size].copy() eps = np.reshape(eps, y.shape) - i0_sd = work_idx['sd'] - sd_beta = work[i0_sd:i0_sd+beta.size].copy() + i0_sd = rwork_idx['sd'] + sd_beta = rwork[i0_sd:i0_sd+beta.size].copy() - i0_vcv = work_idx['vcv'] - cov_beta = work[i0_vcv:i0_vcv+beta.size**2].copy() + i0_vcv = rwork_idx['vcv'] + cov_beta = rwork[i0_vcv:i0_vcv+beta.size**2].copy() cov_beta = np.reshape(cov_beta, (beta.size, beta.size)) result = OdrResult( @@ -510,41 +519,24 @@ def fdummy(beta, x): return np.array([np.nan]) # will never be called yest=y+eps, sd_beta=sd_beta, cov_beta=cov_beta, - res_var=work[work_idx['rvar']], + res_var=rwork[rwork_idx['rvar']], info=info, - stopreason=_interpret_info(info), success=info < 4, nfev=iwork[iwork_idx['nfev']], njev=iwork[iwork_idx['njev']], niter=iwork[iwork_idx['niter']], irank=iwork[iwork_idx['irank']], - inv_condnum=work[work_idx['rcond']], - sum_square=work[work_idx['wss']], - sum_square_delta=work[work_idx['wssde']], - sum_square_eps=work[work_idx['wssep']], + inv_condnum=rwork[rwork_idx['rcond']], + sum_square=rwork[rwork_idx['wss']], + sum_square_delta=rwork[rwork_idx['wssde']], + sum_square_eps=rwork[rwork_idx['wssep']], iwork=iwork, - work=work, + rwork=rwork, ) return result def _get_digit(number: int, ndigit: int) -> int: - """Return the `ndigit`-th digit from the right of `number`.""" + """Return the n-th digit from the right of `number`.""" return (number // 10**(ndigit-1)) % 10 - - -def _interpret_info(info: int) -> str: - """Return a message corresponding to the value of `info`.""" - message = "" - if info == 1: - message = "Sum of squares convergence." - elif info == 2: - message = "Parameter convergence." - elif info == 3: - message = "Sum of squares and parameter convergence." - elif info == 4: - message = "Iteration limit reached." - elif info >= 5: - message = "Questionable results or fatal errors detected. See report and error message." - return message diff --git a/src/odrpack/odr_scipy.py b/src/odrpack/odr_scipy.py new file mode 100644 index 0000000..68dd9b7 --- /dev/null +++ b/src/odrpack/odr_scipy.py @@ -0,0 +1,547 @@ +from typing import Callable, Literal + +import numpy as np +from numpy.typing import NDArray + +from odrpack.__odrpack import loc_iwork, loc_rwork +from odrpack.__odrpack import odr as _odr +from odrpack.__odrpack import workspace_dimensions +from odrpack.result import OdrResult + +__all__ = ['odr_fit'] + + +def odr_fit(f: Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float64]], + xdata: NDArray[np.float64], + ydata: NDArray[np.float64], + beta0: NDArray[np.float64], + *, + weight_x: float | NDArray[np.float64] | None = None, + weight_y: float | NDArray[np.float64] | None = None, + bounds: tuple[NDArray[np.float64] | None, + NDArray[np.float64] | None] | None = None, + task: Literal['explicit-ODR', 'implicit-ODR', 'OLS'] = 'explicit-ODR', + fix_beta: NDArray[np.bool] | None = None, + fix_x: NDArray[np.bool] | None = None, + jac_beta: Callable[[NDArray[np.float64], NDArray[np.float64]], + NDArray[np.float64]] | None = None, + jac_x: Callable[[NDArray[np.float64], NDArray[np.float64]], + NDArray[np.float64]] | None = None, + delta0: NDArray[np.float64] | None = None, + diff_scheme: Literal['forward', 'central'] = 'forward', + report: Literal['none', 'short', 'long', 'iteration'] = 'none', + maxit: int = 50, + ndigit: int | None = None, + taufac: float | None = None, + sstol: float | None = None, + partol: float | None = None, + step_beta: NDArray[np.float64] | None = None, + step_delta: NDArray[np.float64] | None = None, + scale_beta: NDArray[np.float64] | None = None, + scale_delta: NDArray[np.float64] | None = None, + rptfile: str | None = None, + errfile: str | None = None, + ) -> OdrResult: + r"""Solve a weighted orthogonal distance regression (ODR) problem, also + known as errors-in-variables regression. + + Parameters + ---------- + f : Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float64]] + Function to be fitted, with the signature `f(beta, x)`. It must return + an array with the same shape as `y`. + xdata : NDArray[np.float64] + Array of shape `(n,)` or `(m, n)` containing the observed values of the + explanatory variable(s). + ydata : NDArray[np.float64] + Array of shape `(n,)` or `(q, n)` containing the observed values of the + response variable(s). When the model is explicit, `ydata` must contain + a value for each observation. To ignore specific values (e.g., missing + data), set the corresponding entry in `weight_y` to zero. When the model + is implicit, `ydata` is not used (but must be defined). + beta0 : NDArray[np.float64] + Array of shape `(npar,)` with the initial guesses of the model parameters, + within the optional bounds specified by `bounds`. + weight_x : float | NDArray[np.float64] | None + Scalar or array specifying how the errors on `xdata` are to be weighted. + If `weight_x` is a scalar, then it is used for all data points. If + `weight_x` is an array of shape `(n,)` and `m==1`, then `weight_x[i]` + represents the weight for `xdata[i]`. If `weight_x` is an array of shape + `(m)`, then it represents the diagonal of the covariant weighting matrix + for all data points. If `weight_x` is an array of shape `(m, m)`, then + it represents the full covariant weighting matrix for all data points. + If `weight_x` is an array of shape `(m, n)`, then `weight_x[:, i]` + represents the diagonal of the covariant weighting matrix for `xdata[:, i]`. + If `weight_x` is an array of shape `(m, m, n)`, then `weight_x[:, :, i]` + represents the full covariant weighting matrix for `xdata[:, i]`. For a + comprehensive description of the options, refer to page 26 of the + ODRPACK95 guide. By default, `weight_x` is set to one for all `xdata` + points. + weight_y : float | NDArray[np.float64] | None + Scalar or array specifying how the errors on `ydata` are to be weighted. + If `weight_y` is a scalar, then it is used for all data points. If + `weight_y` is an array of shape `(n,)` and `q==1`, then `weight_y[i]` + represents the weight for `ydata[i]`. If `weight_y` is an array of shape + `(q)`, then it represents the diagonal of the covariant weighting matrix + for all data points. If `weight_y` is an array of shape `(q, q)`, then + it represents the full covariant weighting matrix for all data points. + If `weight_y` is an array of shape `(q, n)`, then `weight_y[:, i]` + represents the diagonal of the covariant weighting matrix for `ydata[:, i]`. + If `weight_y` is an array of shape `(q, q, n)`, then `weight_y[:, :, i]` + represents the full covariant weighting matrix for `ydata[:, i]`. For a + comprehensive description of the options, refer to page 25 of the + ODRPACK95 guide. By default, `weight_y` is set to one for all `ydata` + points. + bounds : tuple[NDArray[np.float64] | None, NDArray[np.float64] | None] | None + Tuple of arrays with the same shape as `beta0`, specifying the lower and + upper bounds of the model parameters. The first array contains the lower + bounds, and the second contains the upper bounds. By default, the bounds + are set to negative and positive infinity, respectively, for all elements + of `beta`. + task : Literal['explicit-ODR', 'implicit-ODR', 'OLS'] + Specifies the regression task to be performed. `'explicit-ODR'` solves + an orthogonal distance regression problem with an explicit model. + `'implicit-ODR'` handles models defined implicitly. `'OLS'` performs + ordinary least squares fitting. + fix_beta : NDArray[np.bool] | None + Array with the same shape as `beta0`, specifying which elements of `beta` + are to be held fixed. `True` means the parameter is fixed; `False` means + it is adjustable. By default, all elements of `beta` are set to `False`. + fix_x : NDArray[np.bool] | None + Array with the same shape as `xdata`, specifying which elements of `xdata` + are to be held fixed. Alternatively, it can be a rank-1 array of shape + `(m,)` or `(n,)`, in which case it will be broadcast along the other + axis. `True` means the element is fixed; `False` means it is adjustable. + By default, in orthogonal distance regression mode, all elements of + `fix_x` are set to `False`. In ordinary least squares mode (`task='OLS'`), + all `xdata` values are automatically treated as fixed. + jac_beta : Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float64]] | None + Jacobian of the function to be fitted with respect to `beta`, with the + signature `jac_beta(beta, x)`. It must return an array with shape + `(n, npar, q)` or a compatible shape. By default, the Jacobian is + approximated numerically using the finite difference scheme specified + by `diff_scheme`. + jac_x : Callable[[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float64]] | None + Jacobian of the function to be fitted with respect to `x`, with the + signature `jac_x(beta, x)`. It must return an array with shape + `(n, m, q)` or a compatible shape. By default, the Jacobian is approximated + numerically using the finite difference scheme specified by `diff_scheme`. + delta0 : NDArray[np.float64] | None + Array with the same shape as `xdata`, containing the initial guesses of + the errors in the explanatory variable. By default, `delta0` is set to + zero for all elements of `xdata`. + diff_scheme : Literal['forward', 'central'] + Finite difference scheme used to approximate the Jacobian matrices when + the user does not provide `jac_beta` and `jac_x`. The default method is + forward differences. Central differences are generally more accurate but + require one additional `f` evaluation per partial derivative. + report : Literal['none', 'short', 'long', 'iteration'] + Specifies the level of computation reporting. `'none'` disables all output. + `'short'` prints a brief initial and final summary. `'long'` provides a + detailed initial and final summary. `'iteration'` outputs information at + each iteration step in addition to the detailed summaries. This is + useful for debugging or monitoring progress. + maxit : int | None + Maximum number of allowed iterations. + ndigit : int | None + Number of reliable decimal digits in the values computed by the model + function `f` and its Jacobians `jac_beta`, and `jac_x`. By default, + the value is numerically determined by evaluating `f`. + taufac : float | None + Factor ranging from 0 to 1 to initialize the trust region radius. The + default value is 1. Reducing `taufac` may be appropriate if, at the + first iteration, the computed results for the full Gauss-Newton step + cause an overflow, or cause `beta` and/or `delta` to leave the region + of interest. + sstol : float | None + Factor ranging from 0 to 1 specifying the stopping tolerance for the + sum of the squares convergence. The default value is `eps**(1/2)`, + where `eps` is the machine precision in `float64`. + partol : float | None + Factor ranging from 0 to 1 specifying the stopping tolerance for + parameter convergence (i.e., `beta` and `delta`). When the model is + explicit, the default value is `eps**(2/3)`, and when the model is + implicit, the default value is `eps**(1/3)`, where `eps` is the machine + precision in `float64`. + step_beta : NDArray[np.float64] | None + Array with the same shape as `beta0` containing the _relative_ step + sizes used to compute the finite difference derivatives with respect + to the model parameters. By default, the step size is set internally + based on the value of `ndigit` and the type of finite differences + specified by `diff_scheme`. For additional details, refer to pages 31 + and 78 of the ODRPACK95 guide. + step_delta : NDArray[np.float64] | None + Array with the same shape as `xdata`, containing the _relative_ step + sizes used to compute the finite difference derivatives with respect to + the errors in the explanatory variable. Alternatively, it can be a rank-1 + array of shape `(m,)` or `(n,)`, in which case it will be broadcast along + the other axis. By default, step size is set internally based on the value + of `ndigit` and the type of finite differences specified by `diff_scheme`. + For additional details, refer to pages 31 and 78 of the ODRPACK95 guide. + scale_beta : NDArray[np.float64] | None + Array with the same shape as `beta0` containing the scale values of the + model parameters. Scaling is used to improve the numerical stability + of the regression, but does not affect the problem specification. Scaling + should not be confused with the weighting matrices `weight_x` and + `weight_y`. By default, the scale is set internally based on the relative + magnitudes of `beta`. For further details, refer to pages 32 and 84 of + the ODRPACK95 guide. + scale_delta : NDArray[np.float64] | None + Array with the same shape as `xdata`, containing the scale values of the + errors in the explanatory variable. Alternatively, it can be a rank-1 + array of shape `(m,)` or `(n,)`, in which case it will be broadcast along + the other axis. Scaling is used to improve the numerical stability of + the regression, but does not affect the problem specification. Scaling + should not be confused with the weighting matrices `weight_x` and + `weight_y`. By default, the scale is set internally based on the relative + magnitudes of `xdata`. For further details, refer to pages 32 and 85 of + the ODRPACK95 guide. + rptfile : str | None + File name for storing the computation reports, as defined by `report`. + By default, the reports are sent to standard output. + errfile : str | None + File name for storing the error reports, as defined by `report`. By + default, the reports are sent to standard output. + + Returns + ------- + OdrResult + An object containing the results of the regression. + + + References + ---------- + + [1] Jason W. Zwolak, Paul T. Boggs, and Layne T. Watson. + Algorithm 869: ODRPACK95: A weighted orthogonal distance regression code + with bound constraints. ACM Trans. Math. Softw. 33, 4 (August 2007), 27-es. + https://doi.org/10.1145/1268776.1268782 + + [2] Jason W. Zwolak, Paul T. Boggs, and Layne T. Watson. User's Reference + Guide for ODRPACK95, 2005. + https://github.com/HugoMVale/odrpack95/blob/main/original/Doc/guide.pdf + + Examples + -------- + >>> import numpy as np + >>> from odrpack import odr_fit + >>> xdata = np.array([0.982, 1.998, 4.978, 6.01]) + >>> ydata = np.array([2.7, 7.4, 148.0, 403.0]) + >>> beta0 = np.array([2., 0.5]) + >>> bounds = (np.array([0., 0.]), np.array([10., 0.9])) + >>> def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + ... return beta[0] * np.exp(beta[1]*x) + >>> sol = odr_fit(f, xdata, ydata, beta0, bounds=bounds) + >>> sol.beta + array([1.63336897, 0.9 ]) + """ + + # Check xdata and ydata + if xdata.ndim == 1: + m = 1 + elif xdata.ndim == 2: + m = xdata.shape[0] + else: + raise ValueError( + f"`xdata` must be a rank-1 array of shape `(n,)` or a rank-2 array of shape `(m, n)`, but has shape {xdata.shape}.") + + if ydata.ndim == 1: + q = 1 + elif ydata.ndim == 2: + q = ydata.shape[0] + else: + raise ValueError( + f"`ydata` must be a rank-1 array of shape `(n,)` or a rank-2 array of shape `(q, n)`, but has shape {ydata.shape}.") + + if xdata.shape[-1] == ydata.shape[-1]: + n = xdata.shape[-1] + else: + raise ValueError( + f"The last dimension of `xdata` and `ydata` must be identical, but x.shape={xdata.shape} and y.shape={ydata.shape}.") + + # Check beta0 + if beta0.ndim == 1: + npar = beta0.size + beta = beta0.copy() + else: + raise ValueError( + f"`beta0` must be a rank-1 array of shape `(npar,)`, but has shape {beta0.shape}.") + + # Check p bounds + if bounds is not None: + lower, upper = bounds + if lower is not None: + if lower.shape != beta0.shape: + raise ValueError( + "The lower bound `bounds[0]` must have the same shape as `beta0`.") + if np.any(lower >= beta0): + raise ValueError( + "The lower bound `bounds[0]` must be less than `beta0`.") + if upper is not None: + if upper.shape != beta0.shape: + raise ValueError( + "The upper bound `bounds[1]` must have the same shape as `beta0`.") + if np.any(upper <= beta0): + raise ValueError( + "The upper bound `bounds[1]` must be greater than `beta0`.") + else: + lower, upper = None, None + + # Check other beta related arguments + if fix_beta is not None and fix_beta.shape != beta0.shape: + raise ValueError("`fix_beta` must have the same shape as `beta0`.") + + if step_beta is not None and step_beta.shape != beta0.shape: + raise ValueError("`step_beta` must have the same shape as `beta0`.") + + if scale_beta is not None and scale_beta.shape != beta0.shape: + raise ValueError("`scale_beta` must have the same shape as `beta0`.") + + # Check delta0 + if delta0 is not None: + if delta0.shape != xdata.shape: + raise ValueError("`delta0` must have the same shape as `xdata`.") + delta = delta0.copy() + has_delta0 = True + else: + delta = np.zeros(xdata.shape, dtype=np.float64) + has_delta0 = False + + # Check fix_x + if fix_x is not None: + if fix_x.shape == xdata.shape: + ldifx = n + elif fix_x.shape == (m,) and m > 1 and n != m: + ldifx = 1 + elif fix_x.shape == (n,) and m > 1 and n != m: + ldifx = n + fix_x = np.tile(fix_x, (m, 1)) + else: + raise ValueError( + "`fix_x` must either have the same shape as `xdata` or be a rank-1 array of shape `(m,)` or `(n,)`. See page 26 of the ODRPACK95 User Guide.") + else: + ldifx = 1 + + # Check step_delta + if step_delta is not None: + if step_delta.shape == xdata.shape: + ldstpd = n + elif step_delta.shape == (m,) and m > 1 and n != m: + ldstpd = 1 + elif step_delta.shape == (n,) and m > 1 and n != m: + ldstpd = n + step_delta = np.tile(step_delta, (m, 1)) + else: + raise ValueError( + "`step_delta` must either have the same shape as `xdata` or be a rank-1 array of shape `(m,)` or `(n,)`. See page 31 of the ODRPACK95 User Guide.") + else: + ldstpd = 1 + + # Check scale_delta + if scale_delta is not None: + if scale_delta.shape == xdata.shape: + ldscld = n + elif scale_delta.shape == (m,) and m > 1 and n != m: + ldscld = 1 + elif scale_delta.shape == (n,) and m > 1 and n != m: + ldscld = n + scale_delta = np.tile(scale_delta, (m, 1)) + else: + raise ValueError( + "`scale_delta` must either have the same shape as `xdata` or be a rank-1 array of shape `(m,)` or `(n,)`. See page 32 of the ODRPACK95 User Guide.") + else: + ldscld = 1 + + # Check weight_x + if weight_x is not None: + if isinstance(weight_x, (float, int)): + ldwd = 1 + ld2wd = 1 + weight_x = np.full((m,), weight_x, dtype=np.float64) + elif isinstance(weight_x, np.ndarray): + if weight_x.shape == (m,): + ldwd = 1 + ld2wd = 1 + elif weight_x.shape == (m, m): + ldwd = 1 + ld2wd = m + elif weight_x.shape == (m, n) or (weight_x.shape == (n,) and m == 1): + ldwd = n + ld2wd = 1 + elif weight_x.shape in ((m, 1, 1), (m, 1, n), (m, m, 1), (m, m, n)): + ldwd = weight_x.shape[2] + ld2wd = weight_x.shape[1] + else: + raise ValueError( + r"`weight_x` must be a array of shape `(m,)`, `(n,)`, `(m, m)`, `(m, n)`, `(m, 1, 1)`, `(m, 1, n)`, `(m, m, 1)`, or `(m, m, n)`. See page 26 of the ODRPACK95 User Guide.") + else: + raise TypeError("`weight_x` must be a float or an array.") + else: + ldwd = 1 + ld2wd = 1 + + # Check weight_y + if weight_y is not None: + if isinstance(weight_y, (float, int)): + ldwe = 1 + ld2we = 1 + weight_y = np.full((q,), weight_y, dtype=np.float64) + elif isinstance(weight_y, np.ndarray): + if weight_y.shape == (q,): + ldwe = 1 + ld2we = 1 + elif weight_y.shape == (q, q): + ldwe = 1 + ld2we = q + elif weight_y.shape == (q, n) or (weight_y.shape == (n,) and q == 1): + ldwe = n + ld2we = 1 + elif weight_y.shape in ((q, 1, 1), (q, 1, n), (q, q, 1), (q, q, n)): + ldwe = weight_y.shape[2] + ld2we = weight_y.shape[1] + else: + raise ValueError( + r"`weight_y` must be a array of shape `(q,)`, `(n,)`, `(q, q)`, `(q, n)`, `(q, 1, 1)`, `(q, 1, n)`, `(q, q, 1)`, or `(q, q, n)`. See page 25 of the ODRPACK95 User Guide.") + else: + raise TypeError("`weight_y` must be a float or an array.") + else: + ldwe = 1 + ld2we = 1 + + # Check model function + f0 = f(beta0, xdata) + if f0.shape != ydata.shape: + raise ValueError( + "Function `f` must return an array with the same shape as `ydata`.") + + # Check model jacobians + def fdummy(beta, x): return np.array([np.nan]) # will never be called + + if jac_beta is None and jac_x is None: + has_jac = False + jac_beta = fdummy + jac_x = fdummy + elif jac_beta is not None and jac_x is not None: + has_jac = True + jac0_beta = jac_beta(beta0, xdata) + if jac0_beta.shape[-1] != n or jac0_beta.size != n*npar*q: + raise ValueError( + "Function `jac_beta` must return an array with shape `(n, npar, q)` or compatible.") + jac0_x = jac_x(beta0, xdata) + if jac0_x.shape[-1] != n or jac0_x.size != n*m*q: + raise ValueError( + "Function `jac_x` must return an array with shape `(n, m, q)` or compatible.") + else: + raise ValueError( + "Inconsistent arguments for `jac_beta` and `jac_x`. Either both must be provided or none.") + + # Set iprint flag + iprint_mapping = { + 'none': 0, + 'short': 1001, + 'long': 2002, + 'iteration': 2212 + } + iprint = iprint_mapping[report] + + # Set job flag + jobl = [0, 0, 0, 0, 0] + + if task == "explicit-ODR": + jobl[-1] = 0 + is_odr = True + elif task == "implicit-ODR": + jobl[-1] = 1 + is_odr = True + elif task == "OLS": + jobl[-1] = 2 + is_odr = False + else: + raise ValueError( + f"Invalid value for `task`: {task}.") + + if has_jac: + jobl[-2] = 2 + else: + if diff_scheme == "forward": + jobl[-2] = 0 + elif diff_scheme == "central": + jobl[-2] = 1 + else: + raise ValueError( + f"Invalid value for `diff_scheme`: {diff_scheme}.") + + if has_delta0: + jobl[-4] = 1 + + job = sum(jobl[i] * 10 ** (len(jobl) - 1 - i) for i in range(len(jobl))) + + # Allocate work arrays (drop restart possibility) + lrwork, liwork = workspace_dimensions(n, m, q, npar, is_odr) + rwork = np.zeros(lrwork, dtype=np.float64) + iwork = np.zeros(liwork, dtype=np.int32) + + # Convert fix to ifix + ifixb = (~fix_beta).astype(np.int32) if fix_beta is not None else None + ifixx = (~fix_x).astype(np.int32) if fix_x is not None else None + + # Call the ODRPACK95 routine + # Note: beta, delta, work, and iwork are modified in place + info = _odr( + n=n, m=m, q=q, npar=npar, + ldwe=ldwe, ld2we=ld2we, + ldwd=ldwd, ld2wd=ld2wd, + ldifx=ldifx, + ldstpd=ldstpd, ldscld=ldscld, + f=f, fjacb=jac_beta, fjacd=jac_x, + beta=beta, y=ydata, x=xdata, + delta=delta, + we=weight_y, wd=weight_x, ifixb=ifixb, ifixx=ifixx, + lower=lower, upper=upper, + rwork=rwork, iwork=iwork, + job=job, + ndigit=ndigit, taufac=taufac, sstol=sstol, partol=partol, maxit=maxit, + iprint=iprint, errfile=errfile, rptfile=rptfile + ) + + # Indexes of integer and real work arrays + iwork_idx: dict[str, int] = loc_iwork(m, q, npar) + rwork_idx: dict[str, int] = loc_rwork(n, m, q, npar, ldwe, ld2we, is_odr) + + # Return the result + # Extract results without messing up the original work arrays + i0_eps = rwork_idx['eps'] + eps = rwork[i0_eps:i0_eps+ydata.size].copy() + eps = np.reshape(eps, ydata.shape) + + i0_sd = rwork_idx['sd'] + sd_beta = rwork[i0_sd:i0_sd+beta.size].copy() + + i0_vcv = rwork_idx['vcv'] + cov_beta = rwork[i0_vcv:i0_vcv+beta.size**2].copy() + cov_beta = np.reshape(cov_beta, (beta.size, beta.size)) + + result = OdrResult( + beta=beta, + delta=delta, + eps=eps, + xplus=xdata+delta, + yest=ydata+eps, + sd_beta=sd_beta, + cov_beta=cov_beta, + res_var=rwork[rwork_idx['rvar']], + info=info, + success=info < 4, + nfev=iwork[iwork_idx['nfev']], + njev=iwork[iwork_idx['njev']], + niter=iwork[iwork_idx['niter']], + irank=iwork[iwork_idx['irank']], + inv_condnum=rwork[rwork_idx['rcond']], + sum_square=rwork[rwork_idx['wss']], + sum_square_delta=rwork[rwork_idx['wssde']], + sum_square_eps=rwork[rwork_idx['wssep']], + iwork=iwork, + rwork=rwork, + ) + + return result diff --git a/src/odrpack/result.py b/src/odrpack/result.py index a1f4675..03a8c42 100644 --- a/src/odrpack/result.py +++ b/src/odrpack/result.py @@ -53,7 +53,7 @@ class OdrResult(): Sum of squared differences between observed and fitted `y` values. iwork : NDArray[np.int32] Integer workspace array used internally by `odrpack`. - work : NDArray[np.float64] + rwork : NDArray[np.float64] Floating-point workspace array used internally by `odrpack`. """ beta: NDArray[np.float64] @@ -70,10 +70,24 @@ class OdrResult(): irank: int inv_condnum: float info: int - stopreason: str success: bool sum_square: float sum_square_delta: float sum_square_eps: float iwork: NDArray[np.int32] - work: NDArray[np.float64] + rwork: NDArray[np.float64] + + @property + def stopreason(self) -> str: + message = "" + if self.info == 1: + message = "Sum of squares convergence." + elif self.info == 2: + message = "Parameter convergence." + elif self.info == 3: + message = "Sum of squares and parameter convergence." + elif self.info == 4: + message = "Iteration limit reached." + elif self.info >= 5: + message = "Questionable results or fatal errors detected. See report and error message." + return message diff --git a/subprojects/nanobind.wrap b/subprojects/nanobind.wrap index 78e2e7c..96ae0e5 100644 --- a/subprojects/nanobind.wrap +++ b/subprojects/nanobind.wrap @@ -1,13 +1,13 @@ [wrap-file] -directory = nanobind-2.4.0 -source_url = https://github.com/wjakob/nanobind/archive/refs/tags/v2.4.0.tar.gz -source_filename = nanobind-2.4.0.tar.gz -source_hash = bb35deaed7efac5029ed1e33880a415638352f757d49207a8e6013fefb6c49a7 -patch_filename = nanobind_2.4.0-2_patch.zip -patch_url = https://wrapdb.mesonbuild.com/v2/nanobind_2.4.0-2/get_patch -patch_hash = cf493bda0b11ea4e8d9dd42229c3bbdd52af88cc4aedac75a1eccb102b86dd4a -source_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/nanobind_2.4.0-2/nanobind-2.4.0.tar.gz -wrapdb_version = 2.4.0-2 +directory = nanobind-2.7.0 +source_url = https://github.com/wjakob/nanobind/archive/refs/tags/v2.7.0.tar.gz +source_filename = nanobind-2.7.0.tar.gz +source_hash = 6c8c6bf0435b9d8da9312801686affcf34b6dbba142db60feec8d8e220830499 +patch_filename = nanobind_2.7.0-3_patch.zip +patch_url = https://wrapdb.mesonbuild.com/v2/nanobind_2.7.0-3/get_patch +patch_hash = 66d64f242b34e70bbe9edd80e9267e5bb5b6d631a67f1e6e4f7cedf4614dd0e8 +source_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/nanobind_2.7.0-3/nanobind-2.7.0.tar.gz +wrapdb_version = 2.7.0-3 [provide] nanobind = nanobind_dep diff --git a/subprojects/odrpack95.wrap b/subprojects/odrpack95.wrap index 805b1f9..aecd75a 100644 --- a/subprojects/odrpack95.wrap +++ b/subprojects/odrpack95.wrap @@ -1,7 +1,7 @@ [wrap-git] url = https://github.com/HugoMVale/odrpack95.git -revision = v1.0.1 -# revision = HEAD +# revision = v1.0.1 +revision = HEAD depth = 1 [provide] diff --git a/subprojects/robin-map.wrap b/subprojects/robin-map.wrap index 3da2993..0bea538 100644 --- a/subprojects/robin-map.wrap +++ b/subprojects/robin-map.wrap @@ -1,13 +1,13 @@ [wrap-file] -directory = robin-map-1.3.0 -source_url = https://github.com/Tessil/robin-map/archive/refs/tags/v1.3.0.tar.gz -source_filename = robin-map-1.3.0.tar.gz -source_hash = a8424ad3b0affd4c57ed26f0f3d8a29604f0e1f2ef2089f497f614b1c94c7236 -patch_filename = robin-map_1.3.0-1_patch.zip -patch_url = https://wrapdb.mesonbuild.com/v2/robin-map_1.3.0-1/get_patch -patch_hash = 6d090f988541ffb053512607e0942cbd0dbc2a4fa0563e44ff6a37e810b8c739 -source_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/robin-map_1.3.0-1/robin-map-1.3.0.tar.gz -wrapdb_version = 1.3.0-1 +directory = robin-map-1.4.0 +source_url = https://github.com/Tessil/robin-map/archive/refs/tags/v1.4.0.tar.gz +source_filename = robin-map-1.4.0.tar.gz +source_hash = 7930dbf9634acfc02686d87f615c0f4f33135948130b8922331c16d90a03250c +patch_filename = robin-map_1.4.0-1_patch.zip +patch_url = https://wrapdb.mesonbuild.com/v2/robin-map_1.4.0-1/get_patch +patch_hash = feb14b6752b7d439fb2f3ee968e595a9a3de00ef8cb029488af2ceb4f504b95d +source_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/robin-map_1.4.0-1/robin-map-1.4.0.tar.gz +wrapdb_version = 1.4.0-1 [provide] robin-map = robin_map_dep diff --git a/tests/example2.py b/tests/example2.py new file mode 100644 index 0000000..a4d79ad --- /dev/null +++ b/tests/example2.py @@ -0,0 +1,39 @@ +import numpy as np + +from odrpack import odr + +beta0 = np.array([-1.0, -3.0, 0.09, 0.02, 0.08]) +x = [[0.50, -0.12], + [1.20, -0.60], + [1.60, -1.00], + [1.86, -1.40], + [2.12, -2.54], + [2.36, -3.36], + [2.44, -4.00], + [2.36, -4.75], + [2.06, -5.25], + [1.74, -5.64], + [1.34, -5.97], + [0.90, -6.32], + [-0.28, -6.44], + [-0.78, -6.44], + [-1.36, -6.41], + [-1.90, -6.25], + [-2.50, -5.88], + [-2.88, -5.50], + [-3.18, -5.24], + [-3.44, -4.86]] +x = np.array(x).T +y = np.full(x.shape[-1], 4.0) + + +def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + v, h = x + return beta[2]*(v-beta[0])**2 + 2*beta[3]*(v-beta[0])*(h-beta[1]) \ + + beta[4]*(h-beta[1])**2 - 1 + + +sol = odr(f, beta0, y, x, iprint=1001) + +print("\n beta:", sol.beta) +print("\n delta:", sol.delta) diff --git a/tests/test_bindings.py b/tests/test_bindings.py new file mode 100644 index 0000000..bb25762 --- /dev/null +++ b/tests/test_bindings.py @@ -0,0 +1,89 @@ +import numpy as np + +from odrpack.__odrpack import loc_iwork, loc_rwork, odr, workspace_dimensions + + +def test_loc_iwork(): + res = loc_iwork(m=10, q=2, npar=5) + assert len(res) == 23 + assert all(idx >= 0 for idx in res.values()) + + +def test_loc_rwork(): + res = loc_rwork(n=10, m=2, q=2, npar=5, ldwe=1, ld2we=1, isodr=True) + assert len(res) == 52 + assert all(idx >= 0 for idx in res.values()) + + +def test_workspace_dimensions(): + n = 10 + q = 2 + m = 3 + npar = 5 + isodr = True + lrwork, liwork = workspace_dimensions(n, m, q, npar, isodr) + assert lrwork == 770 + assert liwork == 46 + + +def test_dimension_consistency(): + n = 11 + q = 2 + m = 3 + npar = 5 + for isodr in [True, False]: + dims = workspace_dimensions(n, m, q, npar, isodr) + iwork_idx = loc_iwork(m, q, npar) + rwork_idx = loc_rwork(n, m, q, npar, ldwe=1, ld2we=1, isodr=isodr) + assert dims[0] >= rwork_idx['lrwkmn'] + assert dims[1] >= iwork_idx['liwkmn'] + + +def test_odr(): + "example5 from odrpack" + + def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + return beta[0] * np.exp(beta[1]*x) + + def fjacb(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + jac = np.zeros((beta.size, x.size)) + jac[0, :] = np.exp(beta[1]*x) + jac[1, :] = beta[0]*x*np.exp(beta[1]*x) + return jac + + def fjacd(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + return beta[0] * beta[1] * np.exp(beta[1]*x) + + def fdummy(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + return np.array([42.0]) + + beta0 = np.array([2., 0.5]) + lower = np.array([0., 0.]) + upper = np.array([10., 0.9]) + x = np.array([0.982, 1.998, 4.978, 6.01]) + y = np.array([2.7, 7.4, 148.0, 403.0]) + n = x.size + m = 1 + q = 1 + npar = beta0.size + delta0 = np.zeros(x.shape) + + beta_ref = np.array([1.63337602E+00, 9.00000000E-01]) + + # solution without jacobian + beta = beta0.copy() + delta = delta0.copy() + info = odr(n, m, q, npar, 1, 1, 1, 1, 1, 1, 1, + f, fdummy, fdummy, beta, y, x, delta, + lower=lower, upper=upper, job=0) + assert info == 1 + np.allclose(beta, beta_ref) + + # solution with jacobian + beta = beta0.copy() + delta = delta0.copy() + info = odr(n, m, q, npar, 1, 1, 1, 1, 1, 1, 1, + f, fjacb, fjacd, beta, y, x, delta, + lower=lower, upper=upper, job=20) + assert info == 1 + np.allclose(beta, beta_ref) diff --git a/tests/test_odrpack_multiprocessing.py b/tests/test_multiprocessing.py similarity index 74% rename from tests/test_odrpack_multiprocessing.py rename to tests/test_multiprocessing.py index bf9c5f8..da87971 100644 --- a/tests/test_odrpack_multiprocessing.py +++ b/tests/test_multiprocessing.py @@ -3,7 +3,7 @@ import numpy as np -from odrpack import odr +from odrpack import odr_fit DELAY = 0.01 # s @@ -15,7 +15,7 @@ def f1(beta: np.ndarray, x: np.ndarray) -> np.ndarray: beta1 = np.array([1, -2., 0.1, -0.1]) -x1 = np.linspace(-10., 10., 21) +x1 = np.linspace(-10., 10., 21, dtype=np.float64) y1 = f1(beta1, x1) @@ -25,7 +25,7 @@ def f2(beta: np.ndarray, x: np.ndarray) -> np.ndarray: beta2 = np.array([2., 2.]) -x2 = np.linspace(-10., 10., 41) +x2 = np.linspace(-10., 10., 41, dtype=np.float64) x2 = np.vstack((x2, 10+x2/2)) y2 = f2(beta2, x2) @@ -39,26 +39,26 @@ def f3(beta: np.ndarray, x: np.ndarray) -> np.ndarray: beta3 = np.array([1., 2., 3.]) -x3 = np.linspace(-1., 1., 31) +x3 = np.linspace(-1., 1., 31, dtype=np.float64) x3 = np.vstack((x3, np.exp(x3), x3**2)) y3 = f3(beta3, x3) -case1 = (f1, np.ones_like(beta1), y1, x1) -case2 = (f2, np.ones_like(beta2), y2, x2) -case3 = (f3, np.ones_like(beta3), y3, x3) +case1 = (f1, x1, y1, np.ones_like(beta1)) +case2 = (f2, x2, y2, np.ones_like(beta2)) +case3 = (f3, x3, y3, np.ones_like(beta3)) def test_multiple_processes(): # ref solutions - sol1 = odr(*case1) - sol2 = odr(*case2) - sol3 = odr(*case3) + sol1 = odr_fit(*case1) + sol2 = odr_fit(*case2) + sol3 = odr_fit(*case3) # multiple processes pool = Pool() num_jobs = 10 - solutions = pool.starmap(odr, [case1, case2, case3]*num_jobs) + solutions = pool.starmap(odr_fit, [case1, case2, case3]*num_jobs) pool.close() pool.join() diff --git a/tests/test_odrpack_single.py b/tests/test_odr.py similarity index 86% rename from tests/test_odrpack_single.py rename to tests/test_odr.py index 5211be6..37da373 100644 --- a/tests/test_odrpack_single.py +++ b/tests/test_odr.py @@ -1,54 +1,18 @@ import os +from copy import deepcopy import numpy as np import pytest -from copy import deepcopy from odrpack import odr -from odrpack.__odrpack import diwinf, dwinf, workspace_dimensions +from odrpack.__odrpack import loc_rwork SEED = 1234567890 -def test_diwinf(): - res = diwinf(m=10, npar=5, nq=2) - assert len(res) == 23 - assert all(idx >= 0 for idx in res.values()) - - -def test_dwinf(): - res = dwinf(n=10, m=2, npar=5, nq=2, ldwe=1, ld2we=1, isodr=True) - assert len(res) == 52 - assert all(idx >= 0 for idx in res.values()) - - -def test_workspace_dimensions(): - n = 10 - nq = 2 - m = 3 - npar = 5 - isodr = True - dims = workspace_dimensions(n, m, npar, nq, isodr) - assert dims == (770, 46) - assert dims[1] == 20 + 2*npar + nq*(npar + m) - - -def test_dimension_consistency(): - n = 11 - nq = 2 - m = 3 - npar = 5 - for isodr in [True, False]: - dims = workspace_dimensions(n, m, npar, nq, isodr) - iworkidx = diwinf(m, npar, nq) - workidx = dwinf(n, m, npar, nq, ldwe=1, ld2we=1, isodr=isodr) - assert dims[0] >= workidx['lwkmn'] - assert dims[1] >= iworkidx['liwkmn'] - - @pytest.fixture def case1(): - # m=1, nq=1 + # m=1, q=1 def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: return beta[0] + beta[1] * x + beta[2] * x**2 + beta[3] * x**3 @@ -64,7 +28,7 @@ def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: @pytest.fixture def case2(): - # m=2, nq=1 + # m=2, q=1 def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: return (beta[0] * x[0, :])**3 + x[1, :]**beta[1] @@ -81,7 +45,7 @@ def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: @pytest.fixture def case3(): - # m=3, nq=2 + # m=3, q=2 def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: y = np.zeros((2, x.shape[-1])) y[0, :] = (beta[0] * x[0, :])**3 + x[1, :]**beta[1] + np.exp(x[2, :]/2) @@ -340,14 +304,14 @@ def test_we(case1, case3): sol = odr(**case1, we=1e10) assert np.allclose(sol.eps, np.zeros_like(sol.eps)) - # we (n,) and nq==1 + # we (n,) and q==1 we = np.ones_like(case1['y']) fix = (4, 7) we[fix,] = 1e10 sol = odr(**case1, we=we) assert np.allclose(sol.eps[fix,], np.zeros_like(sol.eps[fix,]), atol=ATOL) - # we (nq, n) + # we (q, n) we = np.ones_like(case3['y']) fix = (4, 13) we[:, fix,] = 1e10 @@ -356,13 +320,13 @@ def test_we(case1, case3): assert np.allclose(sol.eps[:, fix,], np.zeros((sol.eps.shape[0], len(fix))), atol=ATOL) - # we (nq, 1, n) + # we (q, 1, n) we = np.expand_dims(we, axis=1) sol = odr(**case3, we=we) assert np.allclose(sol.delta, sol1.delta) assert np.allclose(sol.eps, sol1.eps) - # we (nq,) + # we (q,) we = np.ones(case3['y'].shape[0]) fix = (1,) we[fix,] = 1e10 @@ -371,15 +335,15 @@ def test_we(case1, case3): assert np.allclose(sol.eps[fix, :], np.zeros((len(fix), sol.eps.shape[-1])), atol=ATOL) - # we (nq, 1, 1) + # we (q, 1, 1) we = we[..., np.newaxis, np.newaxis] sol = odr(**case3, we=we) assert np.allclose(sol.delta, sol1.delta) assert np.allclose(sol.eps, sol1.eps) - # we (nq, nq) - nq = case3['y'].shape[0] - we = np.zeros((nq, nq)) + # we (q, q) + q = case3['y'].shape[0] + we = np.zeros((q, q)) np.fill_diagonal(we, 1.) fix = (1,) we[fix, fix] = 1e10 @@ -388,13 +352,13 @@ def test_we(case1, case3): assert np.allclose(sol.eps[fix, :], np.zeros((len(fix), sol.eps.shape[-1])), atol=ATOL) - # we (nq, nq, 1) + # we (q, q, 1) we = we[..., np.newaxis] sol = odr(**case3, we=we) assert np.allclose(sol.delta, sol1.delta) assert np.allclose(sol.eps, sol1.eps) - # we (nq, nq, n) + # we (q, q, n) we = np.tile(we, (1, 1, case3['y'].shape[-1])) sol = odr(**case3, we=we) assert np.allclose(sol.delta, sol1.delta) @@ -420,20 +384,20 @@ def test_parameters(case1): sstol = 0.123 sol = odr(**case1, sstol=sstol) assert sol.info == 1 - idxwork = dwinf(case1['x'].size, 1, case1['beta0'].size, 1, 1, 1, True) - assert np.isclose(sol.work[idxwork['sstol']], sstol) + rwork_idx = loc_rwork(case1['x'].size, 1, 1, case1['beta0'].size, 1, 1, True) + assert np.isclose(sol.rwork[rwork_idx['sstol']], sstol) # partol partol = 0.456 sol = odr(**case1, partol=partol) assert sol.info == 2 - assert np.isclose(sol.work[idxwork['partl']], partol) + assert np.isclose(sol.rwork[rwork_idx['partl']], partol) # taufac taufac = 0.6969 sol = odr(**case1, taufac=taufac) assert sol.info == 1 - assert np.isclose(sol.work[idxwork['taufc']], taufac) + assert np.isclose(sol.rwork[rwork_idx['taufc']], taufac) def test_restart(case1): @@ -441,29 +405,31 @@ def test_restart(case1): # valid restart sol_ref = odr(**case1) sol1 = odr(**case1, maxit=2) - sol2 = odr(**case1, job=10_000, iwork=sol1.iwork, work=sol1.work) + sol2 = odr(**case1, job=10_000, iwork=sol1.iwork, rwork=sol1.rwork) assert sol2.info == 1 assert np.allclose(sol_ref.beta, sol2.beta) # invalid restarts with pytest.raises(ValueError): - _ = odr(**case1, iwork=sol1.iwork, work=sol1.work) + _ = odr(**case1, iwork=sol1.iwork, rwork=sol1.rwork) with pytest.raises(ValueError): _ = odr(**case1, job=10_000) with pytest.raises(ValueError): _ = odr(**case1, job=10_000, iwork=sol1.iwork) with pytest.raises(ValueError): - _ = odr(**case1, job=10_000, work=sol1.work) + _ = odr(**case1, job=10_000, rwork=sol1.rwork) with pytest.raises(ValueError): - _ = odr(**case1, job=10_000, iwork=sol1.iwork, work=np.ones(10000)) + _ = odr(**case1, job=10_000, iwork=sol1.iwork, rwork=np.ones(10000)) with pytest.raises(ValueError): - _ = odr(**case1, job=10_000, iwork=np.ones(10000, dtype=np.int32), work=sol1.work) + _ = odr(**case1, job=10_000, iwork=np.ones(10000, dtype=np.int32), rwork=sol1.rwork) def test_rptfile_and_errfile(case1): - # write to report file rptfile = 'rtptest.txt' + errfile = 'errtest.txt' + + # write to report file for iprint, rptsize in zip([0, 1001], [0, 2600]): if os.path.isfile(rptfile): os.remove(rptfile) @@ -472,11 +438,9 @@ def test_rptfile_and_errfile(case1): and abs(os.path.getsize(rptfile) - rptsize) < 200 # write to error file - errfile = 'errtest.txt' if os.path.isfile(errfile): os.remove(errfile) - _ = odr(**case1, iprint=1001, - errfile=errfile) + _ = odr(**case1, iprint=1001, errfile=errfile) assert os.path.isfile(errfile) # and os.path.getsize(errfile) > 0 # write to report and error file @@ -484,14 +448,18 @@ def test_rptfile_and_errfile(case1): os.remove(rptfile) if os.path.isfile(errfile): os.remove(errfile) - _ = odr(**case1, job=10, iprint=1001, - rptfile=rptfile, - errfile=errfile) + _ = odr(**case1, job=10, iprint=1001, rptfile=rptfile, errfile=errfile) assert os.path.isfile(rptfile) and os.path.getsize(rptfile) > 2500 assert os.path.isfile(errfile) # and os.path.getsize(errfile) > 0 # I can't get the error file to be written to.. + # Clean up + if os.path.isfile(rptfile): + os.remove(rptfile) + if os.path.isfile(errfile): + os.remove(errfile) + def test_jacobians(): @@ -556,6 +524,48 @@ def fjacd(beta: np.ndarray, x: np.ndarray) -> np.ndarray: _ = odr(f, beta0, y, x, job=0, fjacb=fjacb, fjacd=fjacd) +def test_implicit_model(): + + # model and data are from odrpack's example2 + beta0 = np.array([-1.0, -3.0, 0.09, 0.02, 0.08]) + x = [[0.50, -0.12], + [1.20, -0.60], + [1.60, -1.00], + [1.86, -1.40], + [2.12, -2.54], + [2.36, -3.36], + [2.44, -4.00], + [2.36, -4.75], + [2.06, -5.25], + [1.74, -5.64], + [1.34, -5.97], + [0.90, -6.32], + [-0.28, -6.44], + [-0.78, -6.44], + [-1.36, -6.41], + [-1.90, -6.25], + [-2.50, -5.88], + [-2.88, -5.50], + [-3.18, -5.24], + [-3.44, -4.86]] + x = np.array(x).T + y = np.full(x.shape[-1], 0.0) + + def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + v, h = x + return beta[2]*(v-beta[0])**2 + 2*beta[3]*(v-beta[0])*(h-beta[1]) \ + + beta[4]*(h-beta[1])**2 - 1 + + beta_ref = np.array([-9.99380462E-01, + -2.93104890E+00, + 8.75730642E-02, + 1.62299601E-02, + 7.97538109E-02]) + + sol = odr(f, beta0, y, x, job=1, wd=1) + assert np.allclose(sol.beta, beta_ref) + + def add_noise(array, noise, seed): """Adds random noise to an array.""" rng = np.random.default_rng(seed) diff --git a/tests/test_odr_fit.py b/tests/test_odr_fit.py new file mode 100644 index 0000000..17d5c66 --- /dev/null +++ b/tests/test_odr_fit.py @@ -0,0 +1,562 @@ +import os +from copy import deepcopy + +import numpy as np +import pytest + +from odrpack.__odrpack import loc_rwork +from odrpack.odr_scipy import odr_fit + +SEED = 1234567890 + + +def add_noise(array, noise, seed): + """Adds random noise to an array.""" + rng = np.random.default_rng(seed) + return array*(1 + noise*rng.uniform(-1, 1, size=array.shape)) + + +@pytest.fixture +def case1(): + # m=1, q=1 + def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + return beta[0] + beta[1] * x + beta[2] * x**2 + beta[3] * x**3 + + beta_star = np.array([1, -2., 0.1, -0.1]) + x = np.linspace(-10., 10., 21) + y = f(beta_star, x) + + x = add_noise(x, 5e-2, SEED) + y = add_noise(y, 10e-2, SEED) + + return {'xdata': x, 'ydata': y, 'f': f, 'beta0': np.zeros_like(beta_star)} + + +@pytest.fixture +def case2(): + # m=2, q=1 + def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + return (beta[0] * x[0, :])**3 + x[1, :]**beta[1] + + beta_star = np.array([2., 2.]) + x1 = np.linspace(-10., 10., 41) + x = np.vstack((x1, 10+x1/2)) + y = f(beta_star, x) + + x = add_noise(x, 5e-2, SEED) + y = add_noise(y, 10e-2, SEED) + + return {'xdata': x, 'ydata': y, 'f': f, 'beta0': np.array([1., 1.])} + + +@pytest.fixture +def case3(): + # m=3, q=2 + def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + y = np.zeros((2, x.shape[-1])) + y[0, :] = (beta[0] * x[0, :])**3 + x[1, :]**beta[1] + np.exp(x[2, :]/2) + y[1, :] = (beta[2] * x[0, :])**2 + x[1, :]**beta[1] + return y + + beta_star = np.array([1., 2., 3.]) + x1 = np.linspace(-1., 1., 31) + x = np.vstack((x1, np.exp(x1), x1**2)) + y = f(beta_star, x) + + x = add_noise(x, 5e-2, SEED) + y = add_noise(y, 10e-2, SEED) + + return {'xdata': x, 'ydata': y, 'f': f, 'beta0': np.array([5., 5., 5.])} + + +def test_base_cases(case1, case2, case3): + + # case 1 + sol1 = odr_fit(**case1) + assert sol1.success + assert sol1.info == 1 + + # case 2 + sol2 = odr_fit(**case2) + assert sol2.success + assert sol2.info == 1 + + # case 3 + sol3 = odr_fit(**case3) + assert sol3.success + assert sol3.info == 1 + + # invalid inputs: + with pytest.raises(ValueError): + # x and y don't have the same size + _ = odr_fit(f=case1['f'], + xdata=np.ones(case1['xdata'].size+1), + ydata=case1['ydata'], + beta0=case1['beta0']) + with pytest.raises(ValueError): + # x has invalid shape + _ = odr_fit(f=case2['f'], xdata=np.ones((1, 2, 10)), ydata=np.ones(10), + beta0=case2['beta0']) + with pytest.raises(ValueError): + # y has invalid shape + _ = odr_fit(f=case2['f'], xdata=np.ones((2, 10)), ydata=np.ones((1, 2, 10)), + beta0=case2['beta0']) + + +def test_beta0_related(case1): + + # reference + sol1 = odr_fit(**case1) + + # fix some parameters + sol = odr_fit(**case1, + fix_beta=np.array([True, False, False, True])) + assert np.isclose(sol.beta[0], 0) and np.isclose(sol.beta[-1], 0) + + # fix all parameters + sol = odr_fit(**case1, + fix_beta=np.array([True, True, True, True])) + assert np.allclose(sol.beta, [0]*4) + + # user-defined stpb + sol = odr_fit(**case1, + step_beta=np.full(4, 1e-5)) + assert sol.info == 1 + assert np.allclose(sol.beta, sol1.beta) + + # user-defined sclb + sol = odr_fit(**case1, + scale_beta=np.array([2, 2, 20, 20])) + assert sol.info == 1 + assert np.allclose(sol.beta, sol1.beta) + + # invalid inputs: + with pytest.raises(ValueError): + # lower > beta0 + lower = case1['beta0'].copy() + lower[1:] -= 1 + _ = odr_fit(**case1, bounds=(lower, None)) + with pytest.raises(ValueError): + # upper < beta0 + upper = case1['beta0'].copy() + upper[1:] += 1 + _ = odr_fit(**case1, bounds=(None, upper)) + with pytest.raises(ValueError): + # beta0 has invalid shape + _ = odr_fit(f=case1['f'], xdata=case1['xdata'], ydata=case1['ydata'], + beta0=np.zeros((4, 1))) + with pytest.raises(ValueError): + # beta0 has invalid shape + _ = odr_fit(f=case1['f'], xdata=case1['xdata'], ydata=case1['ydata'], + beta0=np.zeros((1, 4))) + with pytest.raises(ValueError): + # lower has invalid shape + _ = odr_fit(**case1, bounds=(np.zeros((1, 4)), None)) + with pytest.raises(ValueError): + # upper has invalid shape + _ = odr_fit(**case1, bounds=(None, np.zeros((1, 4)))) + with pytest.raises(ValueError): + # ifixb has invalid shape + _ = odr_fit(**case1, fix_beta=np.array([True, False, True])) + with pytest.raises(ValueError): + # stpb has invalid shape + _ = odr_fit(**case1, step_beta=np.array([1e-4, 1., 2.])) + with pytest.raises(ValueError): + # sclb has invalid shape + _ = odr_fit(**case1, scale_beta=np.array([1., 1., 1., 1., 1.])) + + +def test_delta0_related(case1, case3): + + # user-defined delta0 + sol = odr_fit(**case1, delta0=np.ones_like(case1['xdata'])) + assert sol.info == 1 + + # fix some x + fix_x = np.zeros_like(case1['xdata'], dtype=np.bool) + fix = (4, 8) + fix_x[fix,] = True + sol = odr_fit(**case1, fix_x=fix_x) + assert np.allclose(sol.delta[fix,], [0, 0]) + + # fix some x, broadcast (m,) + fix_x = np.zeros(case3['xdata'].shape[0], dtype=np.bool) + fix = (1) + fix_x[fix,] = True + sol = odr_fit(**case3, fix_x=fix_x) + assert np.allclose(sol.delta[fix, :], np.zeros(sol.delta.shape[1])) + + # fix some x, broadcast (n,) + fix_x = np.zeros(case3['xdata'].shape[-1], dtype=np.bool) + fix = (2, 7, 13) + fix_x[fix,] = True + sol = odr_fit(**case3, fix_x=fix_x) + assert np.allclose(sol.delta[:, fix], np.zeros((sol.delta.shape[0], len(fix)))) + + # fix all x (n,) + fix_x = np.ones_like(case1['xdata'], dtype=np.bool) + sol = odr_fit(**case1, fix_x=fix_x) + assert np.allclose(sol.delta, np.zeros_like(sol.delta)) + + # fix all x (m, n) + fix_x = np.ones_like(case3['xdata'], dtype=np.bool) + sol = odr_fit(**case3, fix_x=fix_x) + assert np.allclose(sol.delta, np.zeros_like(sol.delta)) + + # user stpd + sol3 = odr_fit(**case3) + for shape in [case3['xdata'].shape, + case3['xdata'].shape[0], + case3['xdata'].shape[-1]]: + step_delta = np.full(shape, 1e-5) + sol = odr_fit(**case3, step_delta=step_delta) + assert np.allclose(sol.delta, sol3.delta) + + # user scld + sol3 = odr_fit(**case3) + for shape in [case3['xdata'].shape, + case3['xdata'].shape[0], + case3['xdata'].shape[-1]]: + scale_delta = np.full(shape, 10.) + sol = odr_fit(**case3, scale_delta=scale_delta) + assert np.allclose(sol.delta, sol3.delta) + + # invalid inputs + with pytest.raises(ValueError): + # fix_x has invalid shape + _ = odr_fit(**case1, fix_x=np.array([True, False, True])) + with pytest.raises(ValueError): + # step_delta has invalid shape + _ = odr_fit(**case3, step_delta=np.array([1e-4, 1.])) + with pytest.raises(ValueError): + # sclale_delta has invalid shape + _ = odr_fit(**case3, scale_delta=np.array([1., 1., 1., 1.])) + with pytest.raises(ValueError): + # delta0 has invalid shape + _ = odr_fit(**case3, delta0=np.zeros_like(case1['ydata'])) + + +def test_weight_x(case1, case3): + + # weight_x scalar + sol = odr_fit(**case1, weight_x=1e10) + assert np.allclose(sol.delta, np.zeros_like(sol.delta)) + + # weight_x (n,) and m==1 + weight_x = np.ones_like(case1['xdata']) + fix = (4, 7) + weight_x[fix,] = 1e10 + sol = odr_fit(**case1, weight_x=weight_x) + assert np.allclose(sol.delta[fix,], np.zeros_like(sol.delta[fix,])) + + # weight_x (m, n) + weight_x = np.ones_like(case3['xdata']) + fix = (4, 13) + weight_x[:, fix,] = 1e10 + sol = odr_fit(**case3, weight_x=weight_x) + sol1 = deepcopy(sol) + assert np.allclose(sol.delta[:, fix,], np.zeros((sol.delta.shape[0], len(fix)))) + + # weight_x (m, 1, n) + weight_x = np.expand_dims(weight_x, axis=1) + sol = odr_fit(**case3, weight_x=weight_x) + assert np.allclose(sol.delta, sol1.delta) + assert np.allclose(sol.eps, sol1.eps) + + # weight_x (m,) + weight_x = np.ones(case3['xdata'].shape[0]) + fix = (1,) + weight_x[fix,] = 1e10 + sol = odr_fit(**case3, weight_x=weight_x) + sol1 = deepcopy(sol) + assert np.allclose(sol.delta[fix, :], np.zeros((len(fix), sol.delta.shape[-1]))) + + # weight_x (m, 1, 1) + weight_x = weight_x[..., np.newaxis, np.newaxis] + sol = odr_fit(**case3, weight_x=weight_x) + assert np.allclose(sol.delta, sol1.delta) + assert np.allclose(sol.eps, sol1.eps) + + # weight_x (m, m) + m = case3['xdata'].shape[0] + weight_x = np.zeros((m, m)) + np.fill_diagonal(weight_x, 1.) + fix = (1,) + weight_x[fix, fix] = 1e10 + sol = odr_fit(**case3, weight_x=weight_x) + sol1 = deepcopy(sol) + assert np.allclose(sol.delta[fix, :], np.zeros((len(fix), sol.delta.shape[-1]))) + + # weight_x (m, m, 1) + weight_x = weight_x[..., np.newaxis] + sol = odr_fit(**case3, weight_x=weight_x) + assert np.allclose(sol.delta, sol1.delta) + assert np.allclose(sol.eps, sol1.eps) + + # weight_x (m, m, n) + weight_x = np.tile(weight_x, (1, 1, case3['xdata'].shape[-1])) + sol = odr_fit(**case3, weight_x=weight_x) + assert np.allclose(sol.delta, sol1.delta) + assert np.allclose(sol.eps, sol1.eps) + + # weight_x has invalid shape + weight_x = np.ones((1, 1, 1)) + with pytest.raises(ValueError): + _ = odr_fit(**case3, weight_x=weight_x) + + # weight_x has invalid tye + with pytest.raises(TypeError): + _ = odr_fit(**case3, weight_x=[1., 1., 1.]) + + +def test_weight_y(case1, case3): + + ATOL = 1e-6 + + # weight_y scalar + sol = odr_fit(**case1, weight_y=1e10) + assert np.allclose(sol.eps, np.zeros_like(sol.eps)) + + # weight_y (n,) and q==1 + weight_y = np.ones_like(case1['ydata']) + fix = (4, 7) + weight_y[fix,] = 1e10 + sol = odr_fit(**case1, weight_y=weight_y) + assert np.allclose(sol.eps[fix,], np.zeros_like(sol.eps[fix,]), atol=ATOL) + + # weight_y (q, n) + weight_y = np.ones_like(case3['ydata']) + fix = (4, 13) + weight_y[:, fix,] = 1e10 + sol = odr_fit(**case3, weight_y=weight_y) + sol1 = deepcopy(sol) + assert np.allclose(sol.eps[:, fix,], np.zeros((sol.eps.shape[0], len(fix))), + atol=ATOL) + + # weight_y (q, 1, n) + weight_y = np.expand_dims(weight_y, axis=1) + sol = odr_fit(**case3, weight_y=weight_y) + assert np.allclose(sol.delta, sol1.delta) + assert np.allclose(sol.eps, sol1.eps) + + # weight_y (q,) + weight_y = np.ones(case3['ydata'].shape[0]) + fix = (1,) + weight_y[fix,] = 1e10 + sol = odr_fit(**case3, weight_y=weight_y) + sol1 = deepcopy(sol) + assert np.allclose(sol.eps[fix, :], np.zeros((len(fix), sol.eps.shape[-1])), + atol=ATOL) + + # weight_y (q, 1, 1) + weight_y = weight_y[..., np.newaxis, np.newaxis] + sol = odr_fit(**case3, weight_y=weight_y) + assert np.allclose(sol.delta, sol1.delta) + assert np.allclose(sol.eps, sol1.eps) + + # weight_y (q, q) + q = case3['ydata'].shape[0] + weight_y = np.zeros((q, q)) + np.fill_diagonal(weight_y, 1.) + fix = (1,) + weight_y[fix, fix] = 1e10 + sol = odr_fit(**case3, weight_y=weight_y) + sol1 = deepcopy(sol) + assert np.allclose(sol.eps[fix, :], np.zeros((len(fix), sol.eps.shape[-1])), + atol=ATOL) + + # weight_y (q, q, 1) + weight_y = weight_y[..., np.newaxis] + sol = odr_fit(**case3, weight_y=weight_y) + assert np.allclose(sol.delta, sol1.delta) + assert np.allclose(sol.eps, sol1.eps) + + # weight_y (q, q, n) + weight_y = np.tile(weight_y, (1, 1, case3['ydata'].shape[-1])) + sol = odr_fit(**case3, weight_y=weight_y) + assert np.allclose(sol.delta, sol1.delta) + assert np.allclose(sol.eps, sol1.eps) + + # weight_y has invalid shape + weight_y = np.ones((1, 1, 1)) + with pytest.raises(ValueError): + _ = odr_fit(**case3, weight_y=weight_y) + + # weight_y has invalid tye + with pytest.raises(TypeError): + _ = odr_fit(**case3, weight_y=[1., 1., 1.]) + + +def test_parameters(case1): + # maxit + sol = odr_fit(**case1, maxit=2) + assert sol.info == 4 + assert 'iteration limit' in sol.stopreason.lower() + + # sstol + sstol = 0.123 + sol = odr_fit(**case1, sstol=sstol) + assert sol.info == 1 + rwork_idx = loc_rwork(case1['xdata'].size, 1, 1, + case1['beta0'].size, 1, 1, True) + assert np.isclose(sol.rwork[rwork_idx['sstol']], sstol) + + # partol + partol = 0.456 + sol = odr_fit(**case1, partol=partol) + assert sol.info == 2 + assert np.isclose(sol.rwork[rwork_idx['partl']], partol) + + # taufac + taufac = 0.6969 + sol = odr_fit(**case1, taufac=taufac) + assert sol.info == 1 + assert np.isclose(sol.rwork[rwork_idx['taufc']], taufac) + + +def test_rptfile_and_errfile(case1): + + rptfile = 'rtptest.txt' + errfile = 'errtest.txt' + + # write to report file + for report, rptsize in zip(['none', 'short'], [0, 2600]): + if os.path.isfile(rptfile): + os.remove(rptfile) + _ = odr_fit(**case1, report=report, rptfile=rptfile) + assert os.path.isfile(rptfile) \ + and abs(os.path.getsize(rptfile) - rptsize) < 200 + + # write to error file + if os.path.isfile(errfile): + os.remove(errfile) + _ = odr_fit(**case1, report='short', errfile=errfile) + assert os.path.isfile(errfile) # and os.path.getsize(errfile) > 0 + + # write to report and error file + if os.path.isfile(rptfile): + os.remove(rptfile) + if os.path.isfile(errfile): + os.remove(errfile) + _ = odr_fit(**case1, diff_scheme='central', report='short', + rptfile=rptfile, errfile=errfile) + assert os.path.isfile(rptfile) and os.path.getsize(rptfile) > 2500 + assert os.path.isfile(errfile) # and os.path.getsize(errfile) > 0 + + # I can't get the error file to be written to.. + + # Clean up + if os.path.isfile(rptfile): + os.remove(rptfile) + if os.path.isfile(errfile): + os.remove(errfile) + + +def test_jacobians(): + + # model and data are from odrpack's example5 + xdata = np.array([0.982, 1.998, 4.978, 6.01]) + ydata = np.array([2.7, 7.4, 148.0, 403.0]) + beta0 = np.array([2., 0.5]) + bounds = (np.array([0., 0.]), np.array([10., 0.9])) + + def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + return beta[0] * np.exp(beta[1]*x) + + def jac_beta(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + jac = np.zeros((beta.size, x.size)) + jac[0, :] = np.exp(beta[1]*x) + jac[1, :] = beta[0]*x*np.exp(beta[1]*x) + return jac + + def jac_x(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + return beta[0] * beta[1] * np.exp(beta[1]*x) + + beta_ref = np.array([1.63337602, 0.9]) + delta_ref = np.array([-0.36886137, -0.31273038, 0.029287, 0.11031505]) + + # without jacobian + for diff_scheme in ['forward', 'central']: + sol = odr_fit(f, xdata, ydata, beta0, bounds=bounds, + diff_scheme=diff_scheme) + assert np.allclose(sol.beta, beta_ref, rtol=1e-4) + assert np.allclose(sol.delta, delta_ref, rtol=1e-3) + + # with jacobian + sol = odr_fit(f, xdata, ydata, beta0, bounds=bounds, + jac_beta=jac_beta, jac_x=jac_x) + assert np.allclose(sol.beta, beta_ref, rtol=1e-4) + assert np.allclose(sol.delta, delta_ref, rtol=1e-3) + + # invalid f shape + with pytest.raises(ValueError): + _ = odr_fit(lambda beta, x: np.reshape(f(beta, x), (-1, 1)), + xdata, ydata, beta0) + # invalid jac_beta shape + with pytest.raises(ValueError): + _ = odr_fit(f, xdata, ydata, beta0, + jac_beta=lambda beta, x: np.reshape(jac_beta(beta, x), (-1, 1)), + jac_x=jac_x) + # invalid jac_x shape + with pytest.raises(ValueError): + _ = odr_fit(f, xdata, ydata, beta0, + jac_beta=jac_beta, + jac_x=lambda beta, x: np.reshape(jac_x(beta, x), (-1, 1))) + # missing jac_beta + with pytest.raises(ValueError): + _ = odr_fit(f, xdata, ydata, beta0, + jac_x=jac_x) + # missing jac_x + with pytest.raises(ValueError): + _ = odr_fit(f, xdata, ydata, beta0, + jac_beta=jac_beta) + + +def test_implicit_model(): + + # model and data are from odrpack's example2 + beta0 = np.array([-1.0, -3.0, 0.09, 0.02, 0.08]) + x = [[0.50, -0.12], + [1.20, -0.60], + [1.60, -1.00], + [1.86, -1.40], + [2.12, -2.54], + [2.36, -3.36], + [2.44, -4.00], + [2.36, -4.75], + [2.06, -5.25], + [1.74, -5.64], + [1.34, -5.97], + [0.90, -6.32], + [-0.28, -6.44], + [-0.78, -6.44], + [-1.36, -6.41], + [-1.90, -6.25], + [-2.50, -5.88], + [-2.88, -5.50], + [-3.18, -5.24], + [-3.44, -4.86]] + xdata = np.array(x).T + ydata = np.full(xdata.shape[-1], 0.0) + + def f(beta: np.ndarray, x: np.ndarray) -> np.ndarray: + v, h = x + return beta[2]*(v-beta[0])**2 + 2*beta[3]*(v-beta[0])*(h-beta[1]) \ + + beta[4]*(h-beta[1])**2 - 1 + + beta_ref = np.array([-9.99380462E-01, + -2.93104890E+00, + 8.75730642E-02, + 1.62299601E-02, + 7.97538109E-02]) + + sol = odr_fit(f, xdata, ydata, beta0, task='implicit-ODR') + assert np.allclose(sol.beta, beta_ref) + + +def test_ols(case1): + + sol1 = odr_fit(**case1, task='OLS') + sol2 = odr_fit(**case1, weight_x=1e100) + assert np.allclose(sol1.beta, sol2.beta) + assert np.allclose(sol1.delta, np.zeros_like(sol1.delta))