diff --git a/.gitignore b/.gitignore index 2aef42d..7bf38a9 100644 --- a/.gitignore +++ b/.gitignore @@ -132,5 +132,8 @@ dmypy.json build_artifacts package/ -#ignore vs code specific settings +# poetry lock file +poetry.lock + +# ignore vs code specific settings .vscode/ diff --git a/docs/examples/WorkingWithGeospatialData.ipynb b/docs/examples/WorkingWithGeospatialData.ipynb index 7ab4627..604191d 100644 --- a/docs/examples/WorkingWithGeospatialData.ipynb +++ b/docs/examples/WorkingWithGeospatialData.ipynb @@ -1682,6 +1682,62 @@ "\n", "This example is simply demonstrative and will not work as-is because this group label doesn't exist for this dataset." ] + }, + { + "cell_type": "markdown", + "id": "94a148c1", + "metadata": {}, + "source": [ + "# Boyce Index\n", + "\n", + "The Boyce index measures the reliability of habitat suitability models by evaluating how predicted suitability values correlate with observed presence data. It focuses on whether areas predicted to have high suitability also contain a proportionately higher number of observed presences, making it useful for evaluating model performance over specific suitability ranges.\n", + "\n", + "AUC (Area Under the Curve) alone is not sufficient in MaxEnt modeling because it evaluates the overall discrimination between presence and absence but does not assess how well the model predicts relative suitability across different areas. It can be misleading, especially when the model performs well in distinguishing presence and background but poorly in identifying areas of high suitability. The Boyce index complements AUC by focusing on this aspect. Based on work of [Hirzel et al. 2006](https://www.whoi.edu/cms/files/hirzel_etal_2006_53457.pdf) (Eq.4)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "f2b5b7f3", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHFCAYAAAAOmtghAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABYtklEQVR4nO3dd1hTd/sG8DvsHWQPEXAPBFGqFauto1pU3KO11m1r1VqraPW1dVfqqlat2uF421r91W2ti2pV0LoFB25QQEEFZEsgyfn9Qc0rgppgwgnJ/bmuXBf55pzkzpGQx+9zhkQQBAFEREREBsJE7ABERERE2sTihoiIiAwKixsiIiIyKCxuiIiIyKCwuCEiIiKDwuKGiIiIDAqLGyIiIjIoLG6IiIjIoLC4ISIiIoPC4oaIVNavXw+JRFLq5urqirfeegu7d+8WO94re/L+zpw5o7XnvH37NiQSCdavX6+15ySiV8PihojKWLduHf755x8cP34cP/zwA0xNTREeHo4//vhD7GhERC9lJnYAItI/AQEBCAkJUd1/5513UK1aNWzcuBHh4eEiJiMiejnO3BDRS1lZWcHCwgLm5ualxjMzMzF69Gh4e3vDwsICNWvWxLRp0yCTyVTLtG/fHvXr18ez1+gVBAG1a9dGly5dVGMymQyzZ89GgwYNYGVlBWdnZ7Rt2xbHjx8vtd7KlSvRpEkTWFtbo1q1aujTpw8SEhIq9N6GDBkCOzs73Lx5E507d4adnR18fHwwceLEUu8DAO7du4d+/frB3t4eUqkU/fv3R1paWrnPe+bMGXTr1g1OTk6wsrJCcHAwfv/9d9Xj6enp8PHxQWhoKIqLi1Xj8fHxsLW1xQcffFCh90NELG6IqBwKhQJyuRzFxcVISUnB+PHjkZ+fjwEDBqiWKSwsRNu2bfHzzz9jwoQJ+PPPPzFw4EAsWLAAvXr1Ui336aef4tq1azh48GCp19i7dy9u3bqFMWPGAADkcjnCwsIwZ84cdO3aFdu3b8f69esRGhqKpKQk1XofffQRxo8fjw4dOmDHjh1YuXIlLl++jNDQUNy/f79C77e4uBjdunVD+/btsXPnTgwbNgxLlizB/PnzVcs8fvwYHTp0wIEDBxAZGYnNmzfDw8MD/fv3L/N8f//9N1q1aoWsrCysXr0aO3fuRJMmTdC/f3/VvjkuLi7YtGkTTp8+jc8//xwAUFBQgL59+6JGjRpYvXp1hd4LEQEQiIj+tW7dOgFAmZulpaWwcuXKUsuuXr1aACD8/vvvpcbnz58vABAOHDggCIIgKBQKoWbNmkL37t1LLRcWFibUqlVLUCqVgiAIws8//ywAEH788cfn5vvnn38EAMLixYtLjScnJwvW1tbC5MmT1Xp/p0+fVo0NHjy43PfRuXNnoV69eqr7q1atEgAIO3fuLLXcyJEjBQDCunXrVGP169cXgoODheLi4lLLdu3aVfD09BQUCoVq7Mn22r59uzB48GDB2tpauHDhwgvfBxG9GGduiKiMn3/+GadPn8bp06exd+9eDB48GGPGjMGKFStUyxw6dAi2trbo06dPqXWHDBkCAKqZGhMTE4wdOxa7d+9WzcDcunUL+/btw+jRoyGRSACUzORYWVlh2LBhz821e/duSCQSDBw4EHK5XHXz8PBAUFAQDh8+XKH3K5FIyuxLFBgYiDt37qju//3337C3t0e3bt1KLff0bBYA3Lx5E1evXsX7778PAKVydu7cGampqbh27Zpq+UmTJqFLly5477338N///hfLly9H48aNK/Q+iKgEixsiKqNBgwYICQlBSEgI3nnnHXz//ffo2LEjJk+ejKysLABARkYGPDw8VMXJE25ubjAzM0NGRoZqbNiwYbC2tla1Wr777jtYW1uXKmQePnwILy8vmJg8/8/S/fv3IQgC3N3dYW5uXup24sQJpKenV+j92tjYwMrKqtSYpaUlCgsLVfczMjLg7u5eZl0PD48yGQEgIiKiTMbRo0cDQKmcEokEQ4YMQWFhITw8PLivDZEW8GgpIlJLYGAg9u/fj+vXr6N58+ZwdnbGyZMnIQhCqQLnwYMHkMvlcHFxUY1JpVIMHjwYP/30EyIiIrBu3ToMGDAAjo6OqmVcXV0RExMDpVL53ALHxcUFEokE0dHRsLS0LPN4eWPa4uzsjFOnTpUZf3aH4ifve+rUqaX2PXpavXr1VD+npqZizJgxaNKkCS5fvoyIiAgsW7ZMi8mJjA9nbohILbGxsQBKihCg5CiovLw87Nixo9RyP//8s+rxp40bNw7p6eno06cPsrKyMHbs2FKPh4WFobCw8IUnw+vatSsEQcDdu3dVM0tP33TZzmnbti1yc3Oxa9euUuO//fZbqfv16tVDnTp1EBcXV27GkJAQ2NvbAyjZcfu9996DRCLB3r17ERkZieXLl2Pbtm06ex9ExoAzN0RUxqVLlyCXywGUtGO2bduGqKgo9OzZE/7+/gCAQYMG4bvvvsPgwYNx+/ZtNG7cGDExMZg3bx46d+6MDh06lHrOunXr4p133sHevXvxxhtvICgoqNTj7733HtatW4dRo0bh2rVraNu2LZRKJU6ePIkGDRrg3XffRatWrfDhhx9i6NChOHPmDNq0aQNbW1ukpqYiJiYGjRs3xscff6yTbTJo0CAsWbIEgwYNwldffYU6depgz5492L9/f5llv//+e4SFhaFTp04YMmQIvL29kZmZiStXruDcuXPYvHkzAGDGjBmIjo7GgQMH4OHhgYkTJ+LIkSMYPnw4goODVduaiDQk8g7NRKRHyjtaSiqVCk2aNBG++eYbobCwsNTyGRkZwqhRowRPT0/BzMxM8PX1FaZOnVpmuSfWr18vABA2bdpU7uOPHz8Wpk+fLtSpU0ewsLAQnJ2dhXbt2gnHjx8vtdzatWuFFi1aCLa2toK1tbVQq1YtYdCgQcKZM2fUen/PHi1la2tbZtkZM2YIz/6JTElJEXr37i3Y2dkJ9vb2Qu/evYXjx4+XOVpKEAQhLi5O6Nevn+Dm5iaYm5sLHh4eQrt27YTVq1cLgiAIBw4cEExMTIQZM2aUWi8jI0OoUaOG8NprrwkymeyF74eIyicRhGfOrEVEpCO9e/fGiRMncPv27TInBCQi0ha2pYhIp2QyGc6dO4dTp05h+/bt+Oabb1jYEJFOceaGiHTq9u3b8Pf3h4ODAwYMGIAVK1bA1NRU7FhEZMBY3BAREZFB4aHgREREZFBY3BAREZFBYXFDREREBsXojpZSKpW4d+8e7O3ty1wTh4iIiPSTIAjIzc196TXoACMsbu7duwcfHx+xYxAREVEFJCcno3r16i9cxuiKmyfXdElOToaDg4PIaYiIiEgdOTk58PHxUX2Pv4jRFTdPWlEODg4sboiIiKoYdXYp4Q7FREREZFBY3BAREZFBYXFDREREBoXFDRERERkUFjdERERkUFjcEBERkUFhcUNEREQGhcUNERERGRQWN0RERGRQjO4MxURERKQbCqWAU4mZeJBbCDd7KzT3d4KpSeVfpJrFDREREb2yfZdSMeuPeKRmF6rGPKVWmBHeEO8EeFZqFraliIiI6JXsu5SKj389V6qwAYC07EJ8/Os57LuUWql5WNwQERFRhSmUAmb9EQ+hnMeejM36Ix4KZXlL6AaLGyIiIqqwU4mZZWZsniYASM0uxKnEzErLxOKGiIiIKuxB7vMLm4ospw0sboiIiKjC3OyttLqcNrC4ISIiogpr7u8EFzuL5z4uQclRU839nSotE4sbIiIiqjBBEGBrWf6ZZZ6c4WZGeMNKPd8NixsiIiKqsB+iE3AnowDW5iZws7cs9ZiH1AqrBjat9PPc8CR+REREVCHX0nKxNOoGAOCrno3RvYk3z1BMREREVVOxQomIzXEoUijRoYEbegZ7QyKRoGUtZ7GjsS1FREREmlt9+BYu3s2G1Noc83o2hkRS+TM0z8PihoiIiDRyJTUHyw6VtKNmd28EN4fKO8xbHSxuiIiISG3FCiUm/h6HYoWATo3c0S3IS+xIZbC4ISIiIrV99/dNxKfmoJqNOeb20K921BMsboiIiEgtl+5mY8WhmwCA2d0D4PrMod/6gsUNERERvVSRvOToKLlSQOfGHugaWLnnrtEEixsiIiJ6qRWHbuBqWi6cbS0wp3uAXrajnmBxQ0RERC90MSUb3x2+BQCY0yMAznb62Y56gsUNERERPZdMrsDEzbFQKAV0DfRE58b62456QtTi5ujRowgPD4eXlxckEgl27Njx0nU2bNiAoKAg2NjYwNPTE0OHDkVGRobuwxIRERmhb/+6gev38+BiZ4HZ3QPEjqMWUYub/Px8BAUFYcWKFWotHxMTg0GDBmH48OG4fPkyNm/ejNOnT2PEiBE6TkpERGR84pKzsPpISTtqbo/GcLK1EDmRekS9tlRYWBjCwsLUXv7EiRPw8/PDuHHjAAD+/v746KOPsGDBAl1FJCIiMkqFxQpM3BwHpQD0aOKFdwI8xI6ktiq1z01oaChSUlKwZ88eCIKA+/fvY8uWLejSpctz15HJZMjJySl1IyIiohdb8td13HyQB1d7S8zs1kjsOBqpcsXNhg0b0L9/f1hYWMDDwwOOjo5Yvnz5c9eJjIyEVCpV3Xx8fCoxMRERUdVzLukRfjyaAACY17MxHG2qRjvqiSpV3MTHx2PcuHGYPn06zp49i3379iExMRGjRo167jpTp05Fdna26pacnFyJiYmIiKqWwmIFIv5tR/Vq6o23G7qLHUljou5zo6nIyEi0atUKkyZNAgAEBgbC1tYWrVu3xty5c+HpWfbwNEtLS1ha6vfx+ERERPpi8YFrSHiYD3cHS8zoWrXaUU9UqZmbgoICmJiUjmxqagoAEARBjEhEREQG48ztTPwUkwgA+LpXIKQ25iInqhhRi5u8vDzExsYiNjYWAJCYmIjY2FgkJSUBKGkpDRo0SLV8eHg4tm3bhlWrViEhIQHHjh3DuHHj0Lx5c3h56d8l14mIiKqKx0Ul7ShBAPo2q4629d3EjlRhoralzpw5g7Zt26ruT5gwAQAwePBgrF+/HqmpqapCBwCGDBmC3NxcrFixAhMnToSjoyPatWuH+fPnV3p2IiIiQ7Jg/1XcziiAp9QKX3RtKHacVyIRjKyfk5OTA6lUiuzsbDg4OIgdh4iISHQnEzLw7o8nIAjAf4c1x5t1XcWOVIYm399Vap8bIiIi0q6CIjkmbbkAQQDefc1HLwsbTbG4ISIiMmLz915FUmYBvKRWmNalgdhxtILFDRERkZE6fisd//3nDgBgQZ8g2FtVzaOjnsXihoiIyAjly+SYvOUCAOD9FjXwRh0XkRNpD4sbIiIiIxS59wpSHj1G9WrWmNrZMNpRT7C4ISIiMjIxN9Lx64mSU60s6BMIO8sqdcGCl2JxQ0REZERyC4vx+daSdtSglr4IrWU47agnWNwQEREZkXl7ruBu1mPUcLLB5+/UFzuOTrC4ISIiMhJHrz/ExlPJAICFfQJha2DtqCdY3BARERmBnKfaUUNC/dCiprPIiXSHxQ0REZERmLs7HqnZhfBztsHkd+qJHUenWNwQEREZuL+vPsDvZ1IgkQAL+wbBxsIw21FPsLghIiIyYNkFxZiyraQdNbyVP17zcxI5ke6xuCEiIjJgs3fH436ODDVdbBHRybDbUU+wuCEiIjJQf8Xfx9ZzKTD5tx1lZW4qdqRKweKGiIjIAGUVFGHq9osAgJGta6KZbzWRE1UeFjdEREQGaOauy3iYK0MtV1t89nZdseNUKhY3REREBmb/5TTsiL0HEwmwuF8To2lHPcHihoiIyIBk5hdh2r/tqI/erIUmPo7iBhIBixsiIiIDMmPXZaTnFaGuux3Gd6gjdhxRsLghIiIyEHsvpuKPuHswNZFgUd8gWJoZVzvqCRY3REREBiAjT4YvdlwCAIx+qxYCqzuKG0hELG6IiIgMwPSdl5GRX4T6Hvb4pJ1xtqOeYHFDRERUxe2+cA9/XkyF2b/tKAsz4/56N+53T0REVMU9zJXhy3/bUWPa1kaAt1TkROJjcUNERFRFCYKAL3ZcxKOCYjT0dMCYtrXFjqQXWNwQERFVUbvi7mH/5fswN2U76mncCkRERFXQg5xCTN95GQDwSbs6aOjlIHIi/cHihoiIqIoRBAH/2X4R2Y+LEeDtgI/fqiV2JL3C4oaIiKiK2X7+Lv668gDmphIs7tsE5qb8On8atwYREVEVcj+nEDN3lbSjxneoi3oe9iIn0j8sboiIiKoIQRAwddtF5BTKEVRdio/a1BQ7kl5icUNERFRFbDmbgkNXH8DC1ASL+gbBjO2ocnGrEBERVQGp2Y8x+494AMCEjnVRx53tqOdhcUNERKTnBEHA51svIlcmR3ANR4xszXbUi7C4ISIi0nO/n0nG0esPYWFmgoV9gmBqIhE7kl5jcUNERKTH7mY9xpzdVwAAkzrWQ203O5ET6T8WN0RERHpKEAR8vuUC8mRyNPOthmFv+IsdqUpgcUNERKSnfjuVhJib6bAyN8HCPoFsR6mJxQ0REZEeSs4swLw/S9pRkzvVR01XtqPUxeKGiIhIzyiVAj7fegH5RQo093PCkFA/sSNVKSxuiIiI9MyGk3dw/FYGrM1NsaBPIEzYjtIIixsiIiI9kpRRgHl7rgIApoTVh5+LrciJqh5Ri5ujR48iPDwcXl5ekEgk2LFjx0vXkclkmDZtGnx9fWFpaYlatWph7dq1ug9LRESkY0qlgIgtcXhcrMDrNZ3wweu+YkeqkszEfPH8/HwEBQVh6NCh6N27t1rr9OvXD/fv38eaNWtQu3ZtPHjwAHK5XMdJiYiIdO/nf27jVGImbCxMsbBPENtRFSRqcRMWFoawsDC1l9+3bx+OHDmChIQEODk5AQD8/Px0lI6IiKjy3E7Px9f7StpRUzs3gI+TjciJqq4qtc/Nrl27EBISggULFsDb2xt169ZFREQEHj9+/Nx1ZDIZcnJySt2IiIj0iUIpIGJzHAqLlWhV2xnvN68hdqQqTdSZG00lJCQgJiYGVlZW2L59O9LT0zF69GhkZmY+d7+byMhIzJo1q5KTEhERqW/dsUScufMIthammN+bR0e9qio1c6NUKiGRSLBhwwY0b94cnTt3xjfffIP169c/d/Zm6tSpyM7OVt2Sk5MrOTUREdHzJTzMw8L91wAA07o0RPVqbEe9qio1c+Pp6Qlvb29IpVLVWIMGDSAIAlJSUlCnTp0y61haWsLS0rIyYxIREanlSTtKJleidR0XvNfcR+xIBqFKzdy0atUK9+7dQ15enmrs+vXrMDExQfXq1UVMRkREpLk1MQk4l5QFe0szzO8dCImE7ShtELW4ycvLQ2xsLGJjYwEAiYmJiI2NRVJSEoCSltKgQYNUyw8YMADOzs4YOnQo4uPjcfToUUyaNAnDhg2DtbW1GG+BiIioQm4+yMWiA9cBAF92bQgvR36PaYuoxc2ZM2cQHByM4OBgAMCECRMQHByM6dOnAwBSU1NVhQ4A2NnZISoqCllZWQgJCcH777+P8PBwLFu2TJT8REREFSFXKDFx8wUUyZV4q54r+oaw+6BNEkEQBLFDVKacnBxIpVJkZ2fDwcFB7DhERGSEVh2+hfn7rsLeygxRn70JD6mV2JH0nibf31VqnxsiIqKq7vr9XCyJKmlHzQhvxMJGB1jcEBERVZJihRITf49DkUKJ9vXd0Lupt9iRDBKLGyIiokry/ZFbuHg3G1Jrc8zr1ZhHR+lIhc5zc+vWLSxduhRXrlyBRCJBgwYN8Omnn6JWrVrazkdERGQQrqbl4NuDNwAAM7s1hLsD21G6ovHMzf79+9GwYUOcOnUKgYGBCAgIwMmTJ9GoUSNERUXpIiMREVGV9qQdVawQ8HZDd/RownaULml8tFRwcDA6deqEr7/+utT4lClTcODAAZw7d06rAbWNR0sREVFl+/avG1jy13U42pjjwGdt4GbPWRtN6fRoqStXrmD48OFlxocNG4b4+HhNn46IiMigXb6XjeWHStpRs7sHsLCpBBoXN66urqozCj8tNjYWbm5u2shERERkEIrkJe0ouVJAWIAHwgM9xY5kFDTeoXjkyJH48MMPkZCQgNDQUEgkEsTExGD+/PmYOHGiLjISERFVSSv+vomrablwsrXAnB4BPDqqkmhc3Hz55Zewt7fH4sWLMXXqVACAl5cXZs6ciXHjxmk9IBERUVV06W42vvv7JgBgTvcAuNhZipzIeLzS5Rdyc3MBAPb29loLpGvcoZiIiHRNJleg2/JjuHY/F10CPfHdgKZiR6ryNPn+rtB5bp6oSkUNERFRZVl28Aau3c+Fi50F5nQPEDuO0VGruGnatCkOHjyIatWqITg4+IU9Q30/FJyIiEiX4pKzsPpIAgBgbo8AONlaiJzI+KhV3HTv3h2Wlpaqn7lDFBERUVmFxQpEbI6DQimgW5AX3gng0VFieKV9bqoi7nNDRES68vXeq1h95BZc7CwR9VkbVOOsjdbo9CR+NWvWREZGRpnxrKws1KxZU9OnIyIiMgjnkh7hh6O3AADzegawsBGRxsXN7du3oVAoyozLZDKkpKRoJRQREVFV8qQdpRSAnsHe6NjIQ+xIRk3to6V27dql+nn//v2QSqWq+wqFAgcPHoS/v7920xEREVUB30RdR8LDfLjZW2JGeEOx4xg9tYubHj16AAAkEgkGDx5c6jFzc3P4+flh8eLFWg1HRESk787eycSP0SVHR0X2agxHG7ajxKZ2caNUKgEA/v7+OH36NFxcXHQWioiIqCp4XKRAxOYLEASgT7PqaN/AXexIhAqcxC8xMVEXOYiIiKqchfuvITE9Hx4OVviyK9tR+qJCZyjOz8/HkSNHkJSUhKKiolKP8fpSRERkDE4lZmLd8ZL/8Ef2bgyptbnIiegJjYub8+fPo3PnzigoKEB+fj6cnJyQnp4OGxsbuLm5sbghIiKDV1Akx6QtcRAEoH+ID9rWcxM7Ej1F40PBP/vsM4SHhyMzMxPW1tY4ceIE7ty5g2bNmmHRokW6yEhERKRXFuy7hjsZBfCSWmFa1wZix6FnaFzcxMbGYuLEiTA1NYWpqSlkMhl8fHywYMEC/Oc//9FFRiIiIr3xz60MrD9+GwAwv08gHKzYjtI3Ghc35ubmqmtLubu7IykpCQAglUpVPxMRERmifJkck7fGAQDea14Dreu4ipyIyqPxPjfBwcE4c+YM6tati7Zt22L69OlIT0/HL7/8gsaNG+siIxERkV74eu9VJGc+hrejNaZ1YTtKX2k8czNv3jx4epZc5XTOnDlwdnbGxx9/jAcPHuCHH37QekAiIiJ9cOxmOn45cQcAsKBPIOwsK3TAMVUCjf5lBEGAq6srGjVqBABwdXXFnj17dBKMiIhIX+QWFmPylgsAgA9e90Wr2jyRrT7TaOZGEATUqVOHF8gkIiKjMm/PVdzNegwfJ2tMCasvdhx6CY2KGxMTE9SpUwcZGRm6ykNERKRXjl5/iI2nSg6YWdA7CLZsR+k9jfe5WbBgASZNmoRLly7pIg8REZHeyCksxpStJe2oIaF+aFnLWeREpA6Ny8+BAweioKAAQUFBsLCwgLW1danHMzMztRaOiIhITF/tvoJ72YXwdbbB5HfqiR2H1KRxcbN06VIdxCAiItIvf197gP87kwyJBFjYJwg2FmxHVRUa/0sNHjxYFzmIiIj0RvbjYkzdehEAMDTUH839nURORJrQeJ8bIiIiQzdndzzScgrh72KLSZ3YjqpqWNwQERE95eCV+9hyNgUSCbCobyCsLUzFjkQaYnFDRET0r6yCIkzdVtKOGtm6Jpr5sh1VFbG4ISIi+tesP+LxIFeGmq62mPB2XbHjUAWxuCEiIgJw4HIatp+/CxMJsKhvEKzM2Y6qqtQ6WqpXr15qP+G2bdsqHIaIiEgMj/KL8J/tJSen/bBNLTStUU3kRPQq1Jq5kUqlqpuDgwMOHjyIM2fOqB4/e/YsDh48CKlUqrOgREREujJj12Wk58lQx80O4zvUETsOvSK1ipt169apbu7u7ujXrx8SExOxbds2bNu2DQkJCXj33Xfh4qLZVVKPHj2K8PBweHl5QSKRYMeOHWqve+zYMZiZmaFJkyYavSYREdHT9l1Kxa64ezA1kbAdZSA03udm7dq1iIiIgKnp//7xTU1NMWHCBKxdu1aj58rPz0dQUBBWrFih0XrZ2dkYNGgQ2rdvr9F6RERET8vIk2Hav+2oUW/WRJCPo7iBSCs0PkOxXC7HlStXUK9e6ZMaXblyBUqlUqPnCgsLQ1hYmKYR8NFHH2HAgAEwNTXVaLaHiIjoadN3XUZGfhHqudtjXHu2owyFxsXN0KFDMWzYMNy8eROvv/46AODEiRP4+uuvMXToUK0HfNa6detw69Yt/Prrr5g7d67OX4+IiAzT7gv38OeFVJiaSLC4XxAszdiOMhQaFzeLFi2Ch4cHlixZgtTUVACAp6cnJk+ejIkTJ2o94NNu3LiBKVOmIDo6GmZm6kWXyWSQyWSq+zk5ObqKR0REVcTDXBm+3FHSjhrTtjYCvHlAjCHRuLgxMTHB5MmTMXnyZFWh4ODgoPVgz1IoFBgwYABmzZqFunXVP7FSZGQkZs2apcNkRERUlQiCgC93XMKjgmI08HTA2La1xY5EWiYRBEHQdCW5XI7Dhw/j1q1bGDBgAOzt7XHv3j04ODjAzs6uYkEkEmzfvh09evQo9/GsrCxUq1at1I7MSqUSgiDA1NQUBw4cQLt27cqsV97MjY+PD7KzsyulKCMiIv2yK+4exm08DzMTCXaObYVGXpy1qQpycnIglUrV+v7WeObmzp07eOedd5CUlASZTIa3334b9vb2WLBgAQoLC7F69eoKB38RBwcHXLx4sdTYypUrcejQIWzZsgX+/v7lrmdpaQlLS0udZCIioqrlQW4hpu8saUd90q4OCxsDpXFx8+mnnyIkJARxcXFwdnZWjffs2RMjRozQ6Lny8vJw8+ZN1f3ExETExsbCyckJNWrUwNSpU3H37l38/PPPMDExQUBAQKn13dzcYGVlVWaciIjoWYIgYNr2S8gqKEYjLweMbltL7EikIxoXNzExMTh27BgsLCxKjfv6+uLu3bsaPdeZM2fQtm1b1f0JEyYAAAYPHoz169cjNTUVSUlJmkYkIiIqY0fsXUTF34e5acnRUeamvLyiodK4uFEqlVAoFGXGU1JSYG9vr9FzvfXWW3jRLj/r169/4fozZ87EzJkzNXpNIiIyPvdzCjFzVzwA4NP2dVDfg/tcGjKNy9a3334bS5cuVd2XSCTIy8vDjBkz0LlzZ21mIyIiemWCIOA/2y4i+3ExGntLMepNtqMMncYzN0uWLEHbtm3RsGFDFBYWYsCAAbhx4wZcXFywceNGXWQkIiKqsK3n7uLg1QewMDXB4n5BMGM7yuBpXNx4eXkhNjYWmzZtwtmzZ6FUKjF8+HC8//77sLa21kVGIiKiCknNfoxZf1wGAHz2dl3Uddds9wmqmjQ+z83Ro0cRGhpa5gzBcrkcx48fR5s2bbQaUNs0OU6eiIiqLkEQMGTdaRy5/hBBPo7YOqolZ22qME2+vzX+V27bti0yMzPLjGdnZ5c68omIiEhMm8+k4Mj1h7AwM8HivoEsbIyIxv/SgiBAIpGUGc/IyICtra1WQhEREb2Ku1mPMWd3ydFRER3rorYb21HGRO19bnr16gWg5OioIUOGlDrrr0KhwIULFxAaGqr9hERERBoQBAFTtl5ArkyOpjUcMfyNmmJHokqmdnEjlZacoloQBNjb25faedjCwgKvv/46Ro4cqf2EREREGth4KhnRN9JhaWaCRX2DYGpStttAhk3t4mbdunUAAD8/P0yaNAk2NjY6C0VERFQRKY8K8NWfJe2oSZ3qoaZrxS7mTFWbxvvcDBo0qNzLLNy4cQO3b9/WRiYiIiKNKZUCJm+5gPwiBV7zq4ahrcq/oDIZPo2LmyFDhuD48eNlxk+ePIkhQ4ZoIxMREZHGNpxKwvFbGbAyN8HCPmxHGTONi5vz58+jVatWZcZff/11xMbGaiMTERGRRpIyChC55woA4PN36sPPhUfvGjONixuJRILc3Nwy49nZ2eVeUJOIiEiXlEoBk7bEoaBIgeb+Thjc0k/sSCQyjYub1q1bIzIyslQho1AoEBkZiTfeeEOr4YiIiF7mlxN3cDIxEzYWpljUJwgmbEcZPY2vLbVgwQK0adMG9erVQ+vWrQEA0dHRyMnJwaFDh7QekIiI6Hlup+fj671XAQBTw+qjhjOP5KUKzNw0bNgQFy5cQL9+/fDgwQPk5uZi0KBBuHr1KgICAnSRkYiIqIwn7ajHxQqE1nLG+y18xY5EekLjmRug5Mrg8+bN03YWIiIita07fhunbz+CrYUp5vcOZDuKVCp0FbHo6GgMHDgQoaGhqnPe/PLLL4iJidFqOCIiovIkPMzDwv0l7aj/dGkAHye2o+h/NC5utm7dik6dOsHa2hrnzp2DTCYDAOTm5nI2h4iIdE6hFDBpywUUFivxRm0XDGheQ+xIpGc0Lm7mzp2L1atX48cff4S5ublqPDQ0FOfOndNqOCIiometjUnE2TuPYGdphvl9AiGRsB1FpWlc3Fy7dg1t2rQpM+7g4ICsrCxtZCIiIirXzQd5WHjgGgDgiy4N4O1o/ZI1yBhpXNx4enri5s2bZcZjYmJQsyYvK09ERLqhUAqI2ByHIrkSbeq6ov9rPmJHIj2lcXHz0Ucf4dNPP8XJkychkUhw7949bNiwARERERg9erQuMhIREeHH6ATEJmfB3soM83s3ZjuKnkvjQ8EnT56M7OxstG3bFoWFhWjTpg0sLS0RERGBsWPH6iIjEREZuRv3c/HNgesAgOldG8JTynYUPZ9EEAShIisWFBQgPj4eSqUSDRs2hJ2dnbaz6UROTg6kUimys7Ph4OAgdhwiInoJuUKJXquO40JKNtrVd8OawSGctTFCmnx/V+gkfgBgY2MDd3d3SCSSKlPYEBFR1fP90QRcSMmGg5UZ5vVkO4peTuN9buRyOb788ktIpVL4+fnB19cXUqkUX3zxBYqLi3WRkYiIjNTVtBws/aukHTWzWyN4SK1ETkRVgcYzN2PHjsX27duxYMECtGzZEgDwzz//YObMmUhPT8fq1au1HpKIiIxPsUKJiM1xKFYI6NDAHT2DvcWORFWExsXNxo0bsWnTJoSFhanGAgMDUaNGDbz77rssboiISCtWHb6FS3dz4Ghjjnm9AtiOIrVp3JaysrKCn59fmXE/Pz9YWFhoIxMRERm5y/eysezgDQDArG6N4GbPdhSpT+PiZsyYMZgzZ47qmlIAIJPJ8NVXX/FQcCIiemVFciUiNl+AXCmgUyN3dAvyEjsSVTEat6XOnz+PgwcPonr16ggKCgIAxMXFoaioCO3bt0evXr1Uy27btk17SYmIyCh89/dNXEnNQTUbc8ztwaOjSHMaFzeOjo7o3bt3qTEfH54Cm4iIXt2lu9n47u+SS/zM6REAV3tLkRNRVaRxcbNu3Tpd5CAiIiMnkysQsTkOcqWAzo090DWQ7SiqGI33ubl8+fJzH9u3b98rhSEiIuO1/OBNXE3LhbOtBeZ0DxA7DlVhGhc3ISEhWL58eakxmUyGsWPHomfPnloLRkRExuNCShZWHbkFAJjbIwDOdmxHUcVpXNxs2LABs2bNQlhYGNLS0hAbG4vg4GAcOnQIx44d00VGIiIyYDK5AhN/j4NCKSA8yAthjT3FjkRVnMbFTa9evXDhwgXI5XIEBASgZcuWeOutt3D27Fk0bdpUFxmJiMiALf3rBm48yIOLnSVmd2skdhwyABoXNwCgUChQVFQEhUIBhUIBDw8PWFpyCpGIiDRzPukRvv+3HfVVzwBUs+XJYOnVaVzcbNq0CYGBgZBKpbh+/Tr+/PNP/PDDD2jdujUSEhJ0kZGIiAxQYXHJ0VFKAejRxAudGnmIHYkMhMbFzfDhwzFv3jzs2rULrq6uePvtt3Hx4kV4e3ujSZMmOohIRESGaEnUddx6mA9Xe0vMZDuKtEjj89ycO3cO9erVKzVWrVo1/P777/jll1+0FoyIiAzX2TuZ+CG6ZLY/smdjONqwHUXao/HMzbOFzdM++OCDVwpDRESG73GRAhGbL0AQgF5NvdGhobvYkcjAqF3cNGzYEJmZmar7H374IR4+fKi6/+DBA9jY2Gj04kePHkV4eDi8vLwgkUiwY8eOFy6/bds2vP3223B1dYWDgwNatmyJ/fv3a/SaREQkrkUHriExPR/uDpaY0ZXtKNI+tYubq1evQi6Xq+5v2rQJubm5qvuCIKCwsFCjF8/Pz0dQUBBWrFih1vJHjx7F22+/jT179uDs2bNo27YtwsPDcf78eY1el4iIxHH6dibWHksEAHzdKxBSG3ORE5Eh0nifmycEQSgzpumVW8PCwhAWFqb28kuXLi11f968edi5cyf++OMPBAcHa/TaRERUuQqK5Ji0OQ6CAPQLqY629d3EjkQGqsLFjT5QKpXIzc2Fk5PTc5eRyWSQyWSq+zk5OZURjYiInrFg3zXcziiAp9QKX3RtKHYcMmBqt6UkEkmZmRlNZ2q0bfHixcjPz0e/fv2eu0xkZCSkUqnq5uPjU4kJiYgIAE4kZGD98dsAgK97B8LBiu0o0h21Z24EQUD79u1hZlayyuPHjxEeHg4Li5LD957eH6cybNy4ETNnzsTOnTvh5vb8qc2pU6diwoQJqvs5OTkscIiIKlG+TI5JW+IAAO8198GbdV1FTkSGTu3iZvr06aVmarp3715mmd69e2sn1Uv83//9H4YPH47NmzejQ4cOL1zW0tKSl4YgIhLR/H1XkZz5GN6O1vhP5wZixyEjoHZxExERATs7O11mUcvGjRsxbNgwbNy4EV26dBE7DhERvcDxm+n4+Z87AID5vQNhz3YUVQK197lxcXFBWFgYVq1ahXv37mnlxfPy8hAbG4vY2FgAQGJiImJjY5GUlASgpKU0aNAg1fIbN27EoEGDsHjxYrz++utIS0tDWloasrOztZKHiIi0J08mx6QtFwAA77eogTfquIiciIyF2sXNtWvX0LlzZ2zduhX+/v547bXXMGfOHFy4cKHCL37mzBkEBwerDuOeMGECgoODMX36dABAamqqqtABgO+//x5yuRxjxoyBp6en6vbpp59WOAMREenGvD1XcDfrMapXs8ZUtqOoEkmE8k5Y8xLZ2dnYs2cPdu7ciX379qFatWro1q0bunfvjjfffBOmpqa6yKoVOTk5kEqlyM7OhoODg9hxiIgMUvSNh/hgzSkAwG8jWyC0Fmdt6NVo8v2t8bWlAEAqleK9997Dpk2bkJ6eju+//x5KpRJDhw6Fq6srNmzYUKHgRERU9eUWFuPzf9tRg1v6srChSlehmZsXOX/+PORyOV577TVtPq3WcOaGiEi3pmy9gE2nk1HDyQb7xreGjUWVPl8s6QlNvr/V+o3TZL8aXgaBiMh4Hb72AJtOJwMAFvYJZGFDolDrt65JkyaQSCQQBOGlZyVWKBRaCUZERFVL9uNiTNl6EQAwtJUfWtR0FjkRGSu19rlJTExEQkICEhMTVUdLrVy5EufPn8f58+excuVK1KpVC1u3btV1XiIi0lNzd8cjLacQfs42mNypvthxyIipNXPj6+ur+rlv375YtmwZOnfurBoLDAyEj48PvvzyS/To0UPrIYmISL8dunofm8+mQCIBFvUNgrWF/h41S4ZP46OlLl68CH9//zLj/v7+iI+P10ooIiKqOrIL/teOGt7KHyF+TiInImOncXHToEEDzJ07F4WFhaoxmUyGuXPnokEDnqSJiMjYzPrjMh7kylDTxRYRneqJHYdI/WtLPbF69WqEh4fDx8cHQUFBAIC4uDhIJBLs3r1b6wGJiEh/RcXfx7bzd2EiARb1C4KVOdtRJD6Ni5vmzZsjMTERv/76K65evQpBENC/f38MGDAAtra2ushIRER66FF+Ef6zvaQdNbJNTTStUU3kREQlKnQCAhsbG3z44YfazkJERFXIzD8u42GuDLXd7PBZh7pixyFSqdDlF3755Re88cYb8PLywp07JZeyX7JkCXbu3KnVcEREpJ/2XUrDzth7Je2ovmxHkX7RuLhZtWoVJkyYgLCwMDx69Eh10r5q1aph6dKl2s5HRER6JjO/CF/sKGlHjXqzFpr4OIobiOgZGhc3y5cvx48//ohp06bBzOx/Xa2QkBBcvHhRq+GIiEj/TN95Cel5RajrbodPO9QROw5RGRoXN4mJieVeP8rS0hL5+flaCUVERPrpzwup2H0hFaYmEizu2wSWZmxHkf7RuLjx9/dHbGxsmfG9e/eiYcOG2shERER6KD1Phi93XgIAjH6rFhpXl4qciKh8Gh8tNWnSJIwZMwaFhYUQBAGnTp3Cxo0bERkZiZ9++kkXGYmISGSCIODLHZeQmV+E+h72+KQd21GkvzQuboYOHQq5XI7JkyejoKAAAwYMgLe3N7799lu8++67ushIREQi230hFXsvpcHMRIJFfYNgYVahg22JKoVEEAShoiunp6dDqVTCzc1Nm5l0KicnB1KpFNnZ2XBwcBA7DhGR3nuQW4iOS44iq6AYn7avg8/e5jltqPJp8v2tcendrl07ZGVlAQBcXFxUhU1OTg7atWuneVoiItJbgiBg2vZLyCooRkNPB4xpW1vsSEQvpXFxc/jwYRQVFZUZLywsRHR0tFZCERGRftgZew9R8fdhbsp2FFUdau9zc+HCBdXP8fHxSEtLU91XKBTYt28fvL29tZuOiIhE8yCnEDN2XQYAjGtXBw292MqnqkHt4qZJkyaQSCSQSCTltp+sra2xfPlyrYYjIiJxCIKA/2y/iOzHxWjsLcWot2qJHYlIbWoXN4mJiRAEATVr1sSpU6fg6uqqeszCwgJubm4wNeXJnIiIDMG2c3fx15UHsDA1waK+QTA3ZTuKqg61ixtfX18AgFKp1FkYIiISX1p2IWb+UdKO+rRDHdTzsBc5EZFmNC7FIyMjsXbt2jLja9euxfz587USioiIxCEIAqZsu4DcQjmCqkvxUZuaYkci0pjGxc3333+P+vXrlxlv1KgRVq9erZVQREQkjs1nU3D42kNYmJW0o8zYjqIqSOPf2rS0NHh6epYZd3V1RWpqqlZCERFR5buX9Rhz/ogHAEx4uy7quLMdRVWTxsWNj48Pjh07Vmb82LFj8PLy0kooIiKqXIIg4POtF5ArkyO4hiNGtmY7iqouja8tNWLECIwfPx7FxcWqQ8IPHjyIyZMnY+LEiVoPSEREurfpdDKib6TD8t92lKmJROxIRBWmcXEzefJkZGZmYvTo0aozFVtZWeHzzz/H1KlTtR6QiIh0K+VRAb768woAYFKneqjlaidyIqJXU+ELZ+bl5eHKlSuwtrZGnTp1YGlpqe1sOsELZxIR/Y8gCBi45iSO3cxAiG81/N9HLTlrQ3pJk+9vjWdunrCzs8Nrr71W0dWJiEgPbDiZhGM3M2BlboKFbEeRgVCruOnVqxfWr18PBwcH9OrV64XLbtu2TSvBiIhIt5IzCzBvT0k7anKn+vB3sRU5EZF2qFXcSKVSSCQS1c9ERFS1KZUCJm2JQ0GRAs39nDAk1E/sSERaU+F9bqoq7nNDRAT8/M9tTN95Gdbmptg3vjV8nTlrQ/pNk+9vnnqSiMjI3MnIR+SeqwCAKWH1WdiQwVGrLRUcHKxqS73MuXPnXikQERHpjlIpYNLmC3hcrMDrNZ3wweu+Ykci0jq1ipsePXqofi4sLMTKlSvRsGFDtGzZEgBw4sQJXL58GaNHj9ZJSCIi0o71x2/j1O1M2FiYYmGfIJjw6CgyQGoVNzNmzFD9PGLECIwbNw5z5swps0xycrJ20xERkdYkpudjwf6SdtR/OjeAj5ONyImIdEPjfW42b96MQYMGlRkfOHAgtm7dqpVQRESkXQqlgEmb41BYrESr2s54v0UNsSMR6YzGxY21tTViYmLKjMfExMDKykoroYiISLvWHUvEmTuPYGthivm9A9Xej5KoKtL4DMXjx4/Hxx9/jLNnz+L1118HULLPzdq1azF9+nStByQioldz80EeFu6/BgD4omtDVK/GdhQZNo1nbqZMmYKff/4Z58+fx7hx4zBu3DicP38e69evx5QpUzR6rqNHjyI8PBxeXl6QSCTYsWPHS9c5cuQImjVrBisrK9SsWROrV6/W9C0QERkNxb8n65PJlWhdxwXvvuYjdiQinavQtaX69euHfv36vfKL5+fnIygoCEOHDkXv3r1funxiYiI6d+6MkSNH4tdff8WxY8cwevRouLq6qrU+EZGx+Sk6AeeTsmBvacZ2FBmNChU3WVlZ2LJlCxISEhAREQEnJyecO3cO7u7u8Pb2Vvt5wsLCEBYWpvbyq1evRo0aNbB06VIAQIMGDXDmzBksWrSIxQ0R0TNu3M/F4qjrAIAvuzaEl6O1yImIKofGxc2FCxfQoUMHSKVS3L59GyNGjICTkxO2b9+OO3fu4Oeff9ZFTgDAP//8g44dO5Ya69SpE9asWYPi4mKYm5uXWUcmk0Emk6nu5+Tk6CwfEZG+kCuUiNgchyK5Em/Vc0XfkOpiRyKqNBrvczNhwgQMGTIEN27cKHV0VFhYGI4eParVcM9KS0uDu7t7qTF3d3fI5XKkp6eXu05kZCSkUqnq5uPDfjMRGb7vjyYgLiUb9lZm+LoX21FkXDQubk6fPo2PPvqozLi3tzfS0tK0EupFnv2APrnu5/M+uFOnTkV2drbqxhMNEpGhu5aWi2//ugEAmBneCB5SnqaDjIvGbSkrK6tyWzvXrl2Dq6urVkI9j4eHR5kC6sGDBzAzM4Ozs3O561haWsLS0lKnuYiI9EXxk3aUQon29d3Qq6n6+0ESGQqNZ266d++O2bNno7i4GEDJjElSUhKmTJmi8516W7ZsiaioqFJjBw4cQEhISLn72xARGZvVh2/h4t1sSK3NMa9XY7ajyChpXNwsWrQIDx8+hJubGx4/fow333wTtWvXhr29Pb766iuNnisvLw+xsbGIjY0FUHKod2xsLJKSkgCUtJSevtTDqFGjcOfOHUyYMAFXrlzB2rVrsWbNGkRERGj6NoiIDE78vRwsO1TSjprVrRHcHdiOIuOkcVvKwcEBMTExOHToEM6dOwelUommTZuiQ4cOGr/4mTNn0LZtW9X9CRMmAAAGDx6M9evXIzU1VVXoAIC/vz/27NmDzz77DN999x28vLywbNkyHgZOREbvSTuqWCGgY0N3dG/iJXYkItFIhCd75KpBLpfDysoKsbGxCAgI0GUuncnJyYFUKkV2djYcHBzEjkNEpBVL/7qOpX/dgKONOQ581gZu9py1IcOiyfe3Rm0pMzMz+Pr6QqFQvFJAIiLSnkt3s7Hi0E0AwOzuASxsyOhpvM/NF198galTpyIzM1MXeYiISANF8pJ2lFwpICzAA+GBnmJHIhKdxvvcLFu2DDdv3oSXlxd8fX1ha2tb6vFz585pLRwREb3Y8kM3cDUtF062FpjTI4BHRxGhAsVN9+7d+eEhItIDF1OysfLwLQDAnO4BcLHjOb2IgAoUNzNnztRBDCIi0oRMrsDEzbFQKAV0CfREF7ajiFTU3uemoKAAY8aMgbe3N9zc3DBgwIDnXs+JiIh069u/buD6/Ty42FlgTveqefQqka6oXdzMmDED69evR5cuXfDuu+8iKioKH3/8sS6zERFROWKTs7D6SEk7am6PxnCytRA5EZF+UbsttW3bNqxZswbvvvsuAGDgwIFo1aoVFAoFTE1NdRaQiIj+p7BYgYm/x0IpAN2beOGdAA+xIxHpHbVnbpKTk9G6dWvV/ebNm8PMzAz37t3TSTAiIipryV/XcethPlzsLDEzvJHYcYj0ktrFjUKhgIVF6alPMzMzyOVyrYciIqKyzt55hB+PJgAA5vUMQDW2o4jKpXZbShAEDBkyBJaW/zvUsLCwEKNGjSp1rptt27ZpNyEREaGwWIFJm+OgFIBewd7o2IjtKKLnUbu4GTx4cJmxgQMHajUMERGVb9H+a0hIz4ebvSVmsB1F9EJqFzfr1q3TZQ4iInqOM7czseZYIgDg696NIbUxFzkRkX7T+NpSRERUeR4XKRCxOQ6CAPRpVh3t6ruLHYlI77G4ISLSYwv2X8XtjAJ4OFjhy64NxY5DVCWwuCEi0lMnEjKw7thtAP+2o6zZjiJSB4sbIiI9VFAkx+QtFwAA777mg7fquYmciKjqYHFDRKSH5u+9iqTMAnhJrTCtSwOx4xBVKSxuiIj0zPFb6fjvP3cAAPP7BMLeiu0oIk2wuCEi0iN5sv+1owa0qIHWdVxFTkRU9bC4ISLSI5F7riDl0WN4O1rjP53ZjiKqCBY3RER6IuZGOjacTAIALOwTCDtLtc+zSkRPYXFDRKQHcguL8fnWknbUB6/7IrS2i8iJiKouFjdERHpg3p4ruJv1GD5O1pgSVl/sOERVGosbIiKRHbn+EBtPJQMAFvYJgi3bUUSvhMUNEZGIcgqLMeXfdtSQUD+8XtNZ5EREVR+LGyIiEc3dHY/U7EL4Ottg8jv1xI5DZBBY3BARieTvqw/w+5kUSCQl7SgbC7ajiLSBxQ0RkQiyC4oxZVtJO2pYK38093cSORGR4WBxQ0Qkglm7L+N+jgw1XWwR0ZHtKCJtYnFDRFTJ/oq/j23n7pa0o/oGwtrCVOxIRAaFxQ0RUSXKKijC1O0XAQAjW9dEM1+2o4i0jcUNEVElmrnrMh7mylDL1RYT3q4rdhwig8Tihoiokuy/nIYdsfdgIgEW9Q2ClTnbUUS6wOKGiKgSZOYXYdq/7agP29RCcI1qIiciMlwsboiIKsGMXZeRnleEOm52GN+hjthxiAwaixsiIh3bczEVf8Tdg6mJhO0ookrA4oaISIfS82T4YsclAMDHb9ZCkI+juIGIjACLGyIiHZq+8xIy84tQ38Men7SvLXYcIqPA4oaISEd2X7iHPRfTVO0oSzO2o4gqA4sbIiIdeJgrw5f/tqPGtK2NAG+pyImIjAeLGyIiLRMEAV/suIhHBcVo4OmAsW3ZjiKqTCxuiIi0bFfcPey/fB9mJhIs7hsECzP+qSWqTKJ/4lauXAl/f39YWVmhWbNmiI6OfuHyGzZsQFBQEGxsbODp6YmhQ4ciIyOjktISEb3Yg5xCTN95GQDwSbs6aOjlIHIiIuMjanHzf//3fxg/fjymTZuG8+fPo3Xr1ggLC0NSUlK5y8fExGDQoEEYPnw4Ll++jM2bN+P06dMYMWJEJScnIipLEAT8Z/tFZD8uRiMvB4xuW0vsSERGSdTi5ptvvsHw4cMxYsQINGjQAEuXLoWPjw9WrVpV7vInTpyAn58fxo0bB39/f7zxxhv46KOPcObMmUpOTkRU1vbzd/HXlQcwN5Vgcb8gmJuKPjlOZJRE++QVFRXh7Nmz6NixY6nxjh074vjx4+WuExoaipSUFOzZsweCIOD+/fvYsmULunTp8tzXkclkyMnJKXUjItK2tOxCzNxV0o4a36Eu6nuwHUUkFtGKm/T0dCgUCri7u5cad3d3R1paWrnrhIaGYsOGDejfvz8sLCzg4eEBR0dHLF++/LmvExkZCalUqrr5+Pho9X0QEQmCgKnbLiCnUI7A6lJ81Kam2JGIjJroc6YSiaTUfUEQyow9ER8fj3HjxmH69Ok4e/Ys9u3bh8TERIwaNeq5zz916lRkZ2erbsnJyVrNT0S05WwK/r72EBamJljUNwhmbEcRicpMrBd2cXGBqalpmVmaBw8elJnNeSIyMhKtWrXCpEmTAACBgYGwtbVF69atMXfuXHh6epZZx9LSEpaWltp/A0REAFKzH2P2H/EAgM/erou67vYiJyIi0f57YWFhgWbNmiEqKqrUeFRUFEJDQ8tdp6CgACYmpSObmpaczlwQBN0EJSJ6DkEQ8PnWi8iVydHExxEjW/uLHYmIIHJbasKECfjpp5+wdu1aXLlyBZ999hmSkpJUbaapU6di0KBBquXDw8Oxbds2rFq1CgkJCTh27BjGjRuH5s2bw8vLS6y3QURG6v9OJ+Po9YewMGM7ikifiNaWAoD+/fsjIyMDs2fPRmpqKgICArBnzx74+voCAFJTU0ud82bIkCHIzc3FihUrMHHiRDg6OqJdu3aYP3++WG+BiIzU3azHmPvnFQBARMe6qO1mJ3IiInpCIhhZPycnJwdSqRTZ2dlwcOChmkSkOUEQ8MGaU4i5mY6mNRyxeVQoTE3KPxCCiLRDk+9vzqESEWnot1NJiLmZDst/21EsbIj0C4sbIiINJGcW4Kt/21GT36mPmq5sRxHpGxY3RERqUioFTN5yAQVFCrzmVw1DQ/3EjkRE5WBxQ0Skpg0n7+CfhAxYmZtgYZ8gmLAdRaSXWNwQEakhKaMA8/ZcBQBMeac+/FxsRU5ERM/D4oaI6CWUSgERW+LwuFiBFv5OGNTST+xIRPQCLG6IiF7iv//cxqnETNhYmLIdRVQFsLghInqB2+n5mL+vpB01Naw+ajjbiJyIiF6GxQ0R0XMolAIiNsehsFiJ0FrOeL+Fr9iRiEgNLG6IiJ5j3bFEnLnzCLYWppjfO5DtKKIqgsUNEVE5bj3Mw8L91wAA07o0hI8T21FEVQWLGyKiZyiUAiZtjoNMrkTrOi54r7mP2JGISAMsboiInrEmJgHnkrJgZ2mGr3sHQiJhO4qoKmFxQ0T0lJsPcrHowHUAwJddG8Db0VrkRESkKRY3RET/kiuUmLj5AorkSrxZ1xX9QtiOIqqKWNwQEf3rh+gExCVnwd7KDF/3bsx2FFEVxeKGiAjA9fu5WBp1AwAwvWtDeErZjiKqqljcEJHRK1YoMfH3OBQplGhX3w19mlUXOxIRvQIWN0Rk9L4/cgsX72bDwcoMkb3YjiKq6ljcEJFRu5Kag28PlrSjZnVvBHcHK5ETEdGrYnFDREarWKFExOY4FCsEdGjgjh5NvMWORERawOKGiIzWyr9v4fK9HDjamGNerwC2o4gMBIsbIjJKl+9lY/mhf9tR3RrBzZ7tKCJDweKGiIxOkbzk6Ci5UsA7jTzQLchL7EhEpEUsbojI6Kw4dANX03JRzcYcc3qwHUVkaFjcEJFRuXQ3G98dvgUAmNMjAK72liInIiJtY3FDREZDJldg4u9xUCgFdGnsia6BbEcRGSIWN0RkNJYdvIFr93PhbGuB2d0biR2HiHSExQ0RGYW45Cys+rcdNbdHAJzt2I4iMlQsbojI4BUWKzBxcxyUAhAe5IWwxp5iRyIiHWJxQ0QGb+lfN3DzQR5c7CwxuxvbUUSGjsUNERm0c0mP8MPRknbUvJ4BqGZrIXIiItI1FjdEZLAKixWI+Lcd1TPYGx0beYgdiYgqAYsbIjJYiw9cQ8LDfLjaW2JGeEOx4xBRJWFxQ0QG6eydTPwUkwgAiOzZGI42bEcRGQszsQMYCoVSwKnETDzILYSbvRWa+zvB1ISndCeqLE9/Bh2tzTFj12UIAtC7aXV0aOgudjwiqkQsbrRg36VUzPojHqnZhaoxT6kVZoQ3xDsBPOSUSNfK+wwCgNTaDNPZjiIyOmxLvaJ9l1Lx8a/nyvxRTcsuxMe/nsO+S6kiJSMyDs/7DAJA9mM5/rmVLkIqIhITZ25egUIpYNYf8RDKeezJ2PSdl9HA04EtKiIdUCgFfLnzcrmfQQCQAJj1RzzebujBzyCREWFx8wpOJWaW+7/Fpz3IleHNhYcrJxARlSIASM0uxKnETLSs5Sx2HCKqJCxuXsGD3BcXNk+YmUj4v0YiHVAoBciVz5u3+R91P6tEZBhY3LwCN3srtZb7ZXgL/q+RSAf+uZWB93488dLl1P2sEpFhEH2H4pUrV8Lf3x9WVlZo1qwZoqOjX7i8TCbDtGnT4OvrC0tLS9SqVQtr166tpLSlNfd3gqfUCs+bk5Gg5Kip5v5OlRmLyGjwM0hE5RG1uPm///s/jB8/HtOmTcP58+fRunVrhIWFISkp6bnr9OvXDwcPHsSaNWtw7do1bNy4EfXr16/E1P9jaiJRnfX02T+uT+7PCG/IlhSRjvAzSETlkQiC8PKGtY60aNECTZs2xapVq1RjDRo0QI8ePRAZGVlm+X379uHdd99FQkICnJwq9j+xnJwcSKVSZGdnw8HBocLZS+XieW6IRMXPIJHh0+T7W7TipqioCDY2Nti8eTN69uypGv/0008RGxuLI0eOlFln9OjRuH79OkJCQvDLL7/A1tYW3bp1w5w5c2Btba3W6+qiuAF4hmIisfEzSGTYNPn+Fm2H4vT0dCgUCri7lz4turu7O9LS0spdJyEhATExMbCyssL27duRnp6O0aNHIzMz87n73chkMshkMtX9nJwc7b2Jp5iaSLjTMJGI+BkkoidE36FYIin9PytBEMqMPaFUKiGRSLBhwwY0b94cnTt3xjfffIP169fj8ePH5a4TGRkJqVSquvn4+Gj9PRAREZH+EK24cXFxgampaZlZmgcPHpSZzXnC09MT3t7ekEqlqrEGDRpAEASkpKSUu87UqVORnZ2tuiUnJ2vvTRAREZHeEa24sbCwQLNmzRAVFVVqPCoqCqGhoeWu06pVK9y7dw95eXmqsevXr8PExATVq1cvdx1LS0s4ODiUuhEREZHhErUtNWHCBPz0009Yu3Ytrly5gs8++wxJSUkYNWoUgJJZl0GDBqmWHzBgAJydnTF06FDEx8fj6NGjmDRpEoYNG6b2DsVERERk2EQ9Q3H//v2RkZGB2bNnIzU1FQEBAdizZw98fX0BAKmpqaXOeWNnZ4eoqCh88sknCAkJgbOzM/r164e5c+eK9RaIiIhIz4h6nhsx6OpQcCIiItIdTb6/RT9aioiIiEibWNwQERGRQWFxQ0RERAZF1B2KxfBkFyNdnamYiIiItO/J97Y6uwobXXGTm5sLADxTMRERURWUm5tb6mS+5TG6o6WUSiXu3bsHe3v7517mwdDk5OTAx8cHycnJPELsJbit1MdtpT5uK/VxW6nP2LaVIAjIzc2Fl5cXTExevFeN0c3cvOhsxoaOZ2hWH7eV+rit1MdtpT5uK/UZ07Z62YzNE9yhmIiIiAwKixsiIiIyKCxujIClpSVmzJgBS0tLsaPoPW4r9XFbqY/bSn3cVurjtno+o9uhmIiIiAwbZ26IiIjIoLC4ISIiIoPC4oaIiIgMCosbIiIiMigsbgzAypUr4e/vDysrKzRr1gzR0dHPXTYmJgatWrWCs7MzrK2tUb9+fSxZsqQS04pLk231tGPHjsHMzAxNmjTRbUA9o8n2Onz4MCQSSZnb1atXKzGxeDT93ZLJZJg2bRp8fX1haWmJWrVqYe3atZWUVlyabKshQ4aU+3vVqFGjSkwsHk1/rzZs2ICgoCDY2NjA09MTQ4cORUZGRiWl1SMCVWmbNm0SzM3NhR9//FGIj48XPv30U8HW1la4c+dOucufO3dO+O2334RLly4JiYmJwi+//CLY2NgI33//fSUnr3yabqsnsrKyhJo1awodO3YUgoKCKiesHtB0e/39998CAOHatWtCamqq6iaXyys5eeWryO9Wt27dhBYtWghRUVFCYmKicPLkSeHYsWOVmFocmm6rrKysUr9PycnJgpOTkzBjxozKDS4CTbdVdHS0YGJiInz77bdCQkKCEB0dLTRq1Ejo0aNHJScXH4ubKq558+bCqFGjSo3Vr19fmDJlitrP0bNnT2HgwIHajqZ3Krqt+vfvL3zxxRfCjBkzjKq40XR7PSluHj16VAnp9Ium22rv3r2CVCoVMjIyKiOeXnnVv1nbt28XJBKJcPv2bV3E0yuabquFCxcKNWvWLDW2bNkyoXr16jrLqK/YlqrCioqKcPbsWXTs2LHUeMeOHXH8+HG1nuP8+fM4fvw43nzzTV1E1BsV3Vbr1q3DrVu3MGPGDF1H1Cuv8rsVHBwMT09PtG/fHn///bcuY+qFimyrXbt2ISQkBAsWLIC3tzfq1q2LiIgIPH78uDIii0Ybf7PWrFmDDh06wNfXVxcR9UZFtlVoaChSUlKwZ88eCIKA+/fvY8uWLejSpUtlRNYrRnfhTEOSnp4OhUIBd3f3UuPu7u5IS0t74brVq1fHw4cPIZfLMXPmTIwYMUKXUUVXkW1148YNTJkyBdHR0TAzM66PSkW2l6enJ3744Qc0a9YMMpkMv/zyC9q3b4/Dhw+jTZs2lRFbFBXZVgkJCYiJiYGVlRW2b9+O9PR0jB49GpmZmQa9382r/M0CgNTUVOzduxe//fabriLqjYpsq9DQUGzYsAH9+/dHYWEh5HI5unXrhuXLl1dGZL1iXH+xDZREIil1XxCEMmPPio6ORl5eHk6cOIEpU6agdu3aeO+993QZUy+ou60UCgUGDBiAWbNmoW7dupUVT+9o8rtVr1491KtXT3W/ZcuWSE5OxqJFiwy6uHlCk22lVCohkUiwYcMG1VWOv/nmG/Tp0wffffcdrK2tdZ5XTBX5mwUA69evh6OjI3r06KGjZPpHk20VHx+PcePGYfr06ejUqRNSU1MxadIkjBo1CmvWrKmMuHqDxU0V5uLiAlNT0zJV/IMHD8pU+8/y9/cHADRu3Bj379/HzJkzDbq40XRb5ebm4syZMzh//jzGjh0LoOQLSRAEmJmZ4cCBA2jXrl2lZBfDq/xuPe3111/Hr7/+qu14eqUi28rT0xPe3t6qwgYAGjRoAEEQkJKSgjp16ug0s1he5fdKEASsXbsWH3zwASwsLHQZUy9UZFtFRkaiVatWmDRpEgAgMDAQtra2aN26NebOnQtPT0+d59YX3OemCrOwsECzZs0QFRVVajwqKgqhoaFqP48gCJDJZNqOp1c03VYODg64ePEiYmNjVbdRo0ahXr16iI2NRYsWLSoruii09bt1/vx5g/+DWpFt1apVK9y7dw95eXmqsevXr8PExATVq1fXaV4xvcrv1ZEjR3Dz5k0MHz5clxH1RkW2VUFBAUxMSn+tm5qaAij5O29UxNmPmbTlyaGCa9asEeLj44Xx48cLtra2qiMJpkyZInzwwQeq5VesWCHs2rVLuH79unD9+nVh7dq1goODgzBt2jSx3kKl0XRbPcvYjpbSdHstWbJE2L59u3D9+nXh0qVLwpQpUwQAwtatW8V6C5VG022Vm5srVK9eXejTp49w+fJl4ciRI0KdOnWEESNGiPUWKk1FP4cDBw4UWrRoUdlxRaXptlq3bp1gZmYmrFy5Urh165YQExMjhISECM2bNxfrLYiGxY0B+O677wRfX1/BwsJCaNq0qXDkyBHVY4MHDxbefPNN1f1ly5YJjRo1EmxsbAQHBwchODhYWLlypaBQKERIXvk02VbPMrbiRhA0217z588XatWqJVhZWQnVqlUT3njjDeHPP/8UIbU4NP3dunLlitChQwfB2tpaqF69ujBhwgShoKCgklOLQ9NtlZWVJVhbWws//PBDJScVn6bbatmyZULDhg0Fa2trwdPTU3j//feFlJSUSk4tPokgGNtcFRERERky7nNDREREBoXFDRERERkUFjdERERkUFjcEBERkUFhcUNEREQGhcUNERERGRQWN0RERGRQWNwQkVreeustjB8//oXL+Pn5YenSpZWS51Xdvn0bEokEsbGxOn+ew4cPQyKRICsrC8D/LgD5xMyZM9GkSZNXykFE/8PihsiADRkypNwrKD/7Zastp0+fxocffqi6L5FIsGPHDo2fpzKKJB8fH6SmpiIgIACA7rYJAISGhiI1NbXUhTKfFhERgYMHD6ruP+/fjYjUw6uCE5HWuLq6ih1BbaampvDw8KiU17KwsHjha9nZ2cHOzq5SshAZA87cEBEyMjLw3nvvoXr16rCxsUHjxo2xcePGMsvJ5XKMHTsWjo6OcHZ2xhdffFHqasNPz7j4+fkBAHr27AmJRKK6f+vWLXTv3h3u7u6ws7PDa6+9hr/++kv1HG+99Rbu3LmDzz77DBKJBBKJ5Lm5Z86ciRo1asDS0hJeXl4YN26c6rHyZo0cHR2xfv16AKXbSbdv30bbtm0BANWqVYNEIsGQIUMAAPv27cMbb7yhes9du3bFrVu3ymS5evUqQkNDYWVlhUaNGuHw4cOqx142K/R0W2rmzJn473//i507d6re/+HDh9GuXTuMHTu21HoZGRmwtLTEoUOHnruNiIwRixsiQmFhIZo1a4bdu3fj0qVL+PDDD/HBBx/g5MmTpZb773//CzMzM5w8eRLLli3DkiVL8NNPP5X7nKdPnwYArFu3Dqmpqar7eXl56Ny5M/766y+cP38enTp1Qnh4OJKSkgAA27ZtQ/Xq1TF79mykpqYiNTW13OffsmULlixZgu+//x43btzAjh070Lhx4wq9fx8fH2zduhUAcO3aNaSmpuLbb78FAOTn52PChAk4ffo0Dh48CBMTE/Ts2RNKpbLUc0yaNAkTJ07E+fPnERoaim7duiEjI0PjLBEREejXrx/eeecd1fsPDQ3FiBEj8Ntvv0Emk6mW3bBhA7y8vFSFGRGVYFuKyMDt3r27TMtDoVCUuu/t7Y2IiAjV/U8++QT79u3D5s2b0aJFC9W4j48PlixZAolEgnr16uHixYtYsmQJRo4cWeZ1n7SoHB0dS7VkgoKCEBQUpLo/d+5cbN++Hbt27cLYsWPh5OQEU1NT2Nvbv7CVk5SUBA8PD3To0AHm5uaoUaMGmjdvruZWKc3U1BROTk4AADc3t1I7+/bu3bvUsmvWrIGbmxvi4+NV++sAwNixY1XLrlq1Cvv27cOaNWswefJkjbLY2dnB2toaMpms1Pvv3bs3PvnkE+zcuRP9+vUDUFI4Dhky5IWzW0TGiDM3RAaubdu2iI2NLXV7drZFoVDgq6++QmBgIJydnWFnZ4cDBw6oZlOeeP3110t9kbZs2RI3btwoUyy9SH5+PiZPnoyGDRvC0dERdnZ2uHr1apnXepm+ffvi8ePHqFmzJkaOHInt27dDLpdr9BzquHXrFgYMGICaNWvCwcEB/v7+AFAmb8uWLVU/m5mZISQkBFeuXNFaDktLSwwcOBBr164FAMTGxiIuLk7VPiOi/+HMDZGBs7W1Re3atUuNpaSklLq/ePFiLFmyBEuXLkXjxo1ha2uL8ePHo6ioSOt5Jk2ahP3792PRokWoXbs2rK2t0adPH41fy8fHB9euXUNUVBT++usvjB49GgsXLsSRI0dgbm4OiURSan8gACguLtY4b3h4OHx8fPDjjz/Cy8sLSqUSAQEBauXV9ozKiBEj0KRJE6SkpGDt2rVo3749fH19tfoaRIaAMzdEhOjoaHTv3h0DBw5EUFAQatasiRs3bpRZ7sSJE2Xu16lTB6ampuU+r7m5eZlZnejoaAwZMgQ9e/ZE48aN4eHhgdu3b5daxsLCQq3ZIGtra3Tr1g3Lli3D4cOH8c8//+DixYsAStpiT++vc+PGDRQUFDz3uSwsLACUbtllZGTgypUr+OKLL9C+fXs0aNAAjx49Knf9p7eNXC7H2bNnUb9+/Ze+h+dlKe/9N27cGCEhIfjxxx/x22+/YdiwYRV6fiJDx+KGiFC7dm1ERUXh+PHjuHLlCj766COkpaWVWS45ORkTJkzAtWvXsHHjRixfvhyffvrpc5/Xz88PBw8eRFpamqooqF27NrZt26ZqqwwYMKDMzrl+fn44evQo7t69i/T09HKfe/369VizZg0uXbqEhIQE/PLLL7C2tlbNZLRr1w4rVqzAuXPncObMGYwaNQrm5ubPzerr6wuJRILdu3fj4cOHyMvLQ7Vq1eDs7IwffvgBN2/exKFDhzBhwoRy1//uu++wfft2XL16FWPGjMGjR48qXHz4+fnhwoULuHbtGtLT00vNOI0YMQJff/01FAoFevbsWaHnJzJ0LG6ICF9++SWaNm2KTp064a233oKHh0e5J5EbNGgQHj9+jObNm2PMmDH45JNPSp2071mLFy9GVFQUfHx8EBwcDABYsmQJqlWrhtDQUISHh6NTp05o2rRpqfVmz56N27dvo1atWs89d46joyN+/PFHtGrVCoGBgTh48CD++OMPODs7q17bx8cHbdq0wYABAxAREQEbG5vnZvX29sasWbMwZcoUuLu7Y+zYsTAxMcGmTZtw9uxZBAQE4LPPPsPChQvLXf/rr7/G/PnzERQUhOjoaOzcuRMuLi7Pfb0XGTlyJOrVq4eQkBC4urri2LFjqsfee+89mJmZYcCAAbCysqrQ8xMZOonwbFOaiIj0VnJyMvz8/HD69OkyRSERlWBxQ0RUBRQXFyM1NRVTpkzBnTt3Ss3mEFFpbEsREVUBx44dg6+vL86ePYvVq1eLHYdIr3HmhoiIiAwKZ26IiIjIoLC4ISIiIoPC4oaIiIgMCosbIiIiMigsboiIiMigsLghIiIig8LihoiIiAwKixsiIiIyKCxuiIiIyKD8P2HsKp397NwkAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "{'F.ratio': array([0.625, 0.625, 1.875]),\n", + " 'Spearman.cor': 0.8660254037844387,\n", + " 'HS': array([[0.1, 0.4],\n", + " [0.4, 0.7],\n", + " [0.7, 1. ]])}" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# presence suitability scores (e.g., predictions at presences)\n", + "suitability_presence = np.array([0.3, 0.7, 0.8, 0.9])\n", + "\n", + "# background suitability scores (e.g., predictions at backgrounds)\n", + "suitability_background = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])\n", + "\n", + "# Call the continuous_boyce_index function to calculate the Boyce index and Spearman correlation\n", + "results = ela.continuous_boyce_index(yobs=suitability_presence, ypred=suitability_background, nbins=3, to_plot=True)\n", + "\n", + "results" + ] } ], "metadata": { @@ -1700,7 +1756,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.16" + "version": "3.11.10" } }, "nbformat": 4, diff --git a/elapid/__init__.py b/elapid/__init__.py index cf4968f..8e00fd1 100644 --- a/elapid/__init__.py +++ b/elapid/__init__.py @@ -1,5 +1,6 @@ """User entrypoint to elapid""" +from elapid.evaluate import boyce_index, continuous_boyce_index, get_intervals, plot_PE_curve from elapid.features import ( CategoricalTransformer, HingeTransformer, diff --git a/elapid/evaluate.py b/elapid/evaluate.py new file mode 100644 index 0000000..7d67326 --- /dev/null +++ b/elapid/evaluate.py @@ -0,0 +1,176 @@ +import math +import warnings +from typing import Dict, List, Tuple, Union + +import geopandas as gpd +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from scipy.stats import spearmanr + + +def plot_PE_curve(x: Union[np.ndarray, List[float]], y: Union[np.ndarray, List[float]]) -> None: + """ + Plots the Predicted/Expected (P/E) curve for the Boyce Index. + + Args: + x (array-like): Habitat suitability values (e.g., interval midpoints). + y (array-like): Predicted/Expected ratios corresponding to the habitat suitability intervals. + + Returns: + None + """ + plt.figure() + plt.plot(x, y, marker="o") + plt.xlabel("Habitat suitability") + plt.ylabel("Predicted/Expected ratio") + plt.title("Boyce Index") + plt.show() + + +def get_intervals( + bins: Union[int, float, List[float], np.ndarray, str] = "default", + value_range: Tuple[float, float] = (0, 1), +) -> np.ndarray: + """ + Generates habitat suitability intervals based on a 'bins' parameter. + Args: + bins (int | float | list | np.ndarray | str, optional): Defines the binning strategy. + - float: Bin width. Number of bins will be computed as ceil(range / width). + - int: Number of bins. + - list or ndarray: Custom bin edges. Must contain at least 2 values. + - "default": Uses 10 equally spaced bins. Default is "default". + value_range (tuple): The (min, max) range over which to create intervals. Default is (0, 1). + + Returns: + np.ndarray: An (N, 2) array of bin intervals, where each row is (lower_bound, upper_bound). N = total number of bins. + """ + mini, maxi = value_range + + if isinstance(bins, float): + div_result = (maxi - mini) / bins + nbins = int(np.ceil(div_result)) + if not np.isclose(div_result, nbins): + warnings.warn( + f"Bin size results in a non-integer number of bins. Using ceil: {div_result} -> {nbins}", + UserWarning, + ) + boundaries = np.linspace(mini, maxi, nbins + 1) + intervals = np.column_stack((boundaries[:-1], boundaries[1:])) + elif isinstance(bins, (list, np.ndarray)): + bins = np.sort(np.array(bins)) + if len(bins) < 2: + raise ValueError("bins list must have at least two elements.") + intervals = np.column_stack((bins[:-1], bins[1:])) + elif isinstance(bins, int) and bins > 1: + boundaries = np.linspace(mini, maxi, bins + 1) + intervals = np.column_stack((boundaries[:-1], boundaries[1:])) + elif bins == "default": + boundaries = np.linspace(mini, maxi, 11) + intervals = np.column_stack((boundaries[:-1], boundaries[1:])) + else: + raise ValueError("Invalid `bins` value. Must be a float (width), int (count), or list of edges.") + + return intervals + + +def boyce_index( + ypred_observed: np.ndarray, ypred_background: np.ndarray, interval: Union[Tuple[float, float], List[float]] +) -> float: + """ + Calculates the Boyce index for a given interval. + + Uses the convention as defined in Hirzel et al. 2006 to compute the ratio of observed to expected frequencies. + + Args: + ypred_observed (np.ndarray): Suitability values at observed locations (e.g., predictions at presence points). + ypred_background (np.ndarray): Suitability values at random locations (e.g., predictions at background points). + interval (tuple or list): Two elements representing the lower and upper bounds of the interval (i.e., habitat suitability). + + Returns: + float: The ratio of observed to expected frequencies for the given interval. + """ + lower, upper = interval + yobs_bin = (ypred_observed >= lower) & (ypred_observed < upper) + ypred_bin = (ypred_background >= lower) & (ypred_background < upper) + + # Include upper edge for last interval + if np.isclose(upper, np.max(ypred_background)): + yobs_bin |= ypred_observed == upper + ypred_bin |= ypred_background == upper + + pi = np.sum(yobs_bin) / len(ypred_observed) + ei = np.sum(ypred_bin) / len(ypred_background) + + return np.nan if ei == 0 else pi / ei + + +def continuous_boyce_index( + ypred_observed: Union[np.ndarray, pd.Series, gpd.GeoSeries], + ypred_background: Union[np.ndarray, pd.Series, gpd.GeoSeries], + bins: Union[int, float, List[float], np.ndarray, str] = "default", + to_plot: bool = False, +) -> Tuple[np.ndarray, float, np.ndarray]: + """ + Compute the continuous Boyce index to evaluate habitat suitability models. + + Uses the convention as defined in Hirzel et al. 2006 to compute the ratio of observed to expected frequencies. + + Args: + ypred_observed (numpy.ndarray | pd.Series | gpd.GeoSeries): Suitability values at observed locations (i.e., presence points). + ypred_background (numpy.ndarray | pd.Series | gpd.GeoSeries): Suitability values at random/background locations. + bins (int | float | list | np.ndarray | str, optional): Defines the binning strategy: + - int: number of bins + - float: bin width + - list/ndarray: custom bin edges + - 'default': 10 equally spaced bins over the prediction range + to_plot (bool, optional): Whether to plot the predicted-to-expected (P/E) curve. Defaults to False. + + Returns: + Tuple: + - f_scores (numpy.ndarray): The Boyce index scores for each interval. + - corr (float): Spearman correlation coefficient between the P/E ratios and the midpoints of the intervals. + - intervals (numpy.ndarray): The intervals used for the Boyce index calculation. + """ + if not isinstance(ypred_background, (np.ndarray, pd.Series, gpd.GeoSeries)): + raise TypeError( + "The 'ypred_background' parameter must be a NumPy array, Pandas Series, or GeoPandas GeoSeries." + ) + if not isinstance(ypred_observed, (np.ndarray, pd.Series, gpd.GeoSeries)): + raise TypeError("The 'ypred_observed' parameter must be a NumPy array, Pandas Series, or GeoPandas GeoSeries.") + + # Remove NaNs + if np.isnan(ypred_background).any(): + warnings.warn("'ypred_background' contains NaN values, which will be ignored.", UserWarning) + ypred_background = ypred_background[~np.isnan(ypred_background)] + if np.isnan(ypred_observed).any(): + warnings.warn("'ypred_observed' contains NaN values, which will be ignored.", UserWarning) + ypred_observed = ypred_observed[~np.isnan(ypred_observed)] + + ypred_background = np.asarray(ypred_background) + ypred_observed = np.asarray(ypred_observed) + + if ypred_background.ndim != 1 or len(ypred_background) == 0: + raise ValueError("'ypred_background' must be a non-empty one-dimensional array.") + if ypred_observed.ndim != 1 or len(ypred_observed) == 0: + raise ValueError("'ypred_observed' must be a non-empty one-dimensional array.") + + mini, maxi = np.min(ypred_background), np.max(ypred_background) + intervals = get_intervals(bins, value_range=(mini, maxi)) + + f_scores = np.array([boyce_index(ypred_observed, ypred_background, interval) for interval in intervals]) + + valid = ~np.isnan(f_scores) + f_valid = f_scores[valid] + intervals_mid = np.mean(intervals[valid], axis=1) + + if np.sum(valid) <= 2: + warnings.warn("Not enough valid intervals to compute Spearman correlation.", UserWarning) + corr = np.nan + else: + corr, _ = spearmanr(f_valid, intervals_mid) + + if to_plot: + plot_PE_curve(x=intervals_mid, y=f_valid) + + return f_scores, corr, intervals diff --git a/tests/test_evaluate.py b/tests/test_evaluate.py new file mode 100644 index 0000000..dd94b37 --- /dev/null +++ b/tests/test_evaluate.py @@ -0,0 +1,139 @@ +import geopandas as gpd +import numpy as np +import pandas as pd +import pytest +from shapely.geometry import Point + +from elapid.evaluate import boyce_index, continuous_boyce_index, get_intervals + + +# Test Case 1: Normal case with random data +def test_normal_case(): + np.random.seed(0) + ypred_background = np.random.rand(1000) + ypred_observed = np.random.choice(ypred_background, size=100, replace=False) + f_scores, corr, intervals = continuous_boyce_index(ypred_observed, ypred_background, bins=10, to_plot=False) + # Check outputs + assert isinstance(f_scores, np.ndarray) + assert isinstance(corr, float) + assert isinstance(intervals, np.ndarray) + assert len(f_scores) == 10 + assert not np.any(np.isnan(f_scores)) + assert -1 <= corr <= 1 + + +# Test Case 2: Empty background array +def test_empty_background(): + ypred_background = np.array([]) + ypred_observed = np.array([0.5, 0.6, 0.7]) + with pytest.raises(ValueError) as exc_info: + continuous_boyce_index(ypred_observed, ypred_background, bins=10) + assert "'ypred_background' must be a non-empty one-dimensional array." in str(exc_info.value) + + +# Test Case 3: Empty observed array +def test_empty_observed(): + ypred_background = np.random.rand(1000) + ypred_observed = np.array([]) + with pytest.raises(ValueError) as exc_info: + continuous_boyce_index(ypred_observed, ypred_background, bins=10) + assert "'ypred_observed' must be a non-empty one-dimensional array." in str(exc_info.value) + + +# Test Case 4: Observed containing NaNs +def test_observed_with_nans(recwarn): + ypred_background = np.random.rand(1000) + ypred_observed = np.random.choice(ypred_background, size=100, replace=False) + ypred_observed[::10] = np.nan + f_scores, corr, intervals = continuous_boyce_index(ypred_observed, ypred_background, bins=10) + # warning should be issued + w = recwarn.pop(UserWarning) + assert "'ypred_observed' contains NaN values, which will be ignored." in str(w.message) + # outputs still valid + assert len(f_scores) == 10 + assert np.isnan(corr) or -1 <= corr <= 1 + + +# Test Case 5: Invalid bins value +def test_invalid_bins(): + ypred_background = np.random.rand(1000) + ypred_observed = np.random.choice(ypred_background, size=100, replace=False) + with pytest.raises(ValueError): + continuous_boyce_index(ypred_observed, ypred_background, bins=1) + + +# Test Case 6: Custom float bin width +def test_custom_bin_width(): + ypred_background = np.random.rand(1000) + ypred_observed = np.random.choice(ypred_background, size=100, replace=False) + f_scores, corr, intervals = continuous_boyce_index(ypred_observed, ypred_background, bins=0.1) + # number of bins should be ceil(range/0.1) + expected_nbins = int(np.ceil((ypred_background.max() - ypred_background.min()) / 0.1)) + assert len(f_scores) == expected_nbins + assert -1 <= corr <= 1 or np.isnan(corr) + + +# Test Case 7: Background containing NaNs +def test_background_with_nans(recwarn): + ypred_background = np.random.rand(1000) + ypred_background[::50] = np.nan + ypred_observed = np.random.choice(ypred_background[~np.isnan(ypred_background)], size=100, replace=False) + f_scores, corr, intervals = continuous_boyce_index(ypred_observed, ypred_background, bins=10) + w = recwarn.pop(UserWarning) + assert "'ypred_background' contains NaN values, which will be ignored." in str(w.message) + assert len(f_scores) == 10 + + +# Test Case 8: Observed outside background range +def test_observed_outside_range(): + ypred_background = np.random.rand(1000) + ypred_observed = np.array([ypred_background.max() + 0.5, ypred_background.max() + 1.0]) + f_scores, corr, intervals = continuous_boyce_index(ypred_observed, ypred_background, bins=10) + assert len(f_scores) == 10 + # all f_scores should be zero or nan + assert all((np.isnan(f) or f == 0) for f in f_scores) + assert np.isnan(corr) or -1 <= corr <= 1 + + +# Test Case 9: Large dataset performance +def test_large_dataset(): + ypred_background = np.random.rand(1000000) + ypred_observed = np.random.choice(ypred_background, size=10000, replace=False) + f_scores, corr, intervals = continuous_boyce_index(ypred_observed, ypred_background, bins=20) + assert len(f_scores) == 20 + assert -1 <= corr <= 1 + + +# Test Case 10: Pandas Series inputs +def test_with_pandas_series(): + np.random.seed(0) + ypred_background = pd.Series(np.random.rand(1000)) + ypred_observed = ypred_background.sample(n=100) + f_scores, corr, intervals = continuous_boyce_index(ypred_observed, ypred_background, bins=10) + assert len(f_scores) == 10 + assert -1 <= corr <= 1 + + +# Test Case 11: GeoPandas GeoSeries inputs +def test_with_geopandas_geoseries(): + np.random.seed(0) + num_points = 500 + coords = np.random.rand(num_points, 2) + suitability = np.random.rand(num_points) + geometry = [Point(xy) for xy in coords] + gdf = gpd.GeoDataFrame({"suitability": suitability}, geometry=geometry) + ypred_background = gdf["suitability"] + ypred_observed = ypred_background.sample(n=50) + f_scores, corr, intervals = continuous_boyce_index(ypred_observed, ypred_background, bins=10) + assert len(f_scores) == 10 + assert -1 <= corr <= 1 + + +# Test Case 12: Both with NaNs leading to empty after removal +def test_empty_after_nan_removal(): + ypred_background = np.array([np.nan, np.nan]) + ypred_observed = np.array([np.nan, np.nan]) + with pytest.raises(ValueError) as exc_info: + continuous_boyce_index(ypred_observed, ypred_background, bins=10) + # after dropping NaNs, arrays empty + assert "must be a non-empty one-dimensional array" in str(exc_info.value)